|
Chapters:
|
MuSR /
EwaldPython
from ase import Atom
from numpy import *
def vsum(lattice,musite,N):
"""This routine calculates Coulomb point charge potential
in the periodic lattice, by the Ewald trick
(see
http://www.fis.unipr.it/~derenzi/dispense/pmwiki.php?n=MuSR.ASEStart#potential
). The result is the inverse Fourier transform of V(Q), calculated
on reciprocal lattice sites. Since the factor phi(Q) is propto
1/Q**2 this sum converges faster than direct sums. It can be limited
within a reciprocal Lorentz sphere of radius G(Angstroems-1)
with G = b, 2b, 3b,..., nb and b=max(b1,b2,b3)
For each value of G a V(r) value is added.
"""
qe=1.60217646e-19 # electron charge in C (abs val)
A2m=1e-10 # Angstroem to m
# to account for distances in (1A=1e-10 m)
# and charges in e=1.60217646e-19 C
epsilon_0=8.85418781e-12 # in SI
reciprocal_cell=lattice.get_reciprocal_cell() # ditto
b1=reciprocal_cell[0,:]
b2=reciprocal_cell[1,:]
b3=reciprocal_cell[2,:]
b=sqrt(dot(reciprocal_cell,reciprocal_cell).diagonal()).max()
# print b # for debugging
# max b vector length
# shorter b vectors requires more cell to reach the same radius
n1=ceil(N*(sqrt(dot(b1,b1))/b)).__int__() # number of cells in b1 direction
n2=ceil(N*(sqrt(dot(b2,b2))/b)).__int__() # number of cells in b2 direction
n3=ceil(N*(sqrt(dot(b3,b3))/b)).__int__() # number of cells in b3 direction
# whose index is number_sites +1 for mu, -1 since indices start from 0
v=zeros((N))
r=lattice.get_positions() # ion positions in the unit cell
rmu=musite.get_position() # supposed muon position
number_sites=lattice.get_number_of_atoms()
rr=r-outer(ones((number_sites,1)),rmu) # in order to compute exp i Q dot rr
q=lattice.get_charges() # charges of each ion in the unit cell
v_cell=lattice.get_volume()
for n in range(-n1,n1):
for m in range(-n2,n2):
for l in range(-n3,n3):
Q=n*b1+m*b2+l*b3 # reciprocal lattice vector
QQ=sqrt(dot(Q,Q)) # its modulus
j=int(floor_divide(QQ,b)) # index of the Lorentz shell
if j < N and j > 0:
rho_c=0.0
for k in range(number_sites):
rho_c=rho_c+q[k]*qe*exp(1j*dot(Q,rr[k,:])) # FT of cell charge density
v[j]=v[j]+2*pi**4/v_cell/epsilon_0*rho_c*(A2m/QQ)**2
return v
|