TABLE OF CONTENTS
ABINIT/make_efg_ion [ Functions ]
NAME
make_efg_ion
FUNCTION
compute the electric field gradient due to ionic cores
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group (JJ) This file is distributed under the terms of the GNU General Public License, see ~ABINIT/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
natom, number of atoms in the unit cell nsym=number of symmetries in space group ntypat, the number of types of atoms in the unit cell rprimd(3,3), the matrix giving the transformation from crystal to cartesian coordinates symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations tnons(3,nsym) = nonsymmorphic translations typat(natom), the type of each atom in the unit cell ucvol, the volume of the unit cell in atomic units xred(3,natom) the location of each atom in the cell in crystallographic coordinates zion(ntypat) the net charge on each type of atom
OUTPUT
efg(3,3,natom), the 3x3 efg tensors at each atomic site
SIDE EFFECTS
NOTES
This routine computes the electric field gradient, specifically the components
$\partial^2 V/\partial x_\alpha \partial x_\beta$ of the potential generated by the ionic cores,
at each atomic site in the unit cell.
Key references:
Profeta, Mauri, and Pickard, ``Accurate first principles prediction of $^{17}$O NMR parameters in
SiO$_2$: Assignment of the zeolite ferrierite spectrum'', J. Am. Chem. Soc. 125, 541--548 (2003);
A. Honma, ``Dipolar lattice-sums with applications to the exciton bands of anthracene crystal and
the crystal field due to point charges'', J. Phys. Soc. Jpn. 42, 1129--1135 (1977);
and Kresse and Joubert, ``From ultrasoft pseudopotentials to the projector augmented wave method'',
Phys. Rev. B. 59, 1758--1775 (1999). In Kresse and Joubert's notation, the ionic cores are $n_{Zc}$;
these charges are given by the net core charges on the pseudoatoms. Due to otherwise slow convergence,
the sum over atoms is carried out by an Ewald method as detailed in the Honma reference, specifically
his Eq. 4.8.
PARENTS
calc_efg
CHILDREN
matpointsym,matr3inv,xred2xcart
SOURCE
55 #if defined HAVE_CONFIG_H 56 #include "config.h" 57 #endif 58 59 #include "abi_common.h" 60 61 subroutine make_efg_ion(efg,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,zion) 62 63 use defs_basis 64 use m_profiling_abi 65 66 use m_geometry, only : xred2xcart 67 use m_special_funcs, only : abi_derfc 68 69 !This section has been created automatically by the script Abilint (TD). 70 !Do not modify the following lines by hand. 71 #undef ABI_FUNC 72 #define ABI_FUNC 'make_efg_ion' 73 use interfaces_32_util 74 use interfaces_45_geomoptim 75 !End of the abilint section 76 77 implicit none 78 79 !Arguments ------------------------------------ 80 !scalars 81 integer,intent(in) :: natom,nsym,ntypat 82 real(dp) :: ucvol 83 !arrays 84 integer,intent(in) :: symrel(3,3,nsym),typat(natom) 85 real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym) 86 real(dp),intent(in) :: zion(ntypat) 87 real(dp),intent(inout) :: xred(3,natom) 88 real(dp),intent(out) :: efg(3,3,natom) 89 !Local variables------------------------------- 90 !scalars 91 integer :: iatom,ishell,ii,jatom,jj,nshell,sx,sy,sz 92 real(dp) :: cph,dampfac,derfc_karg,derivs,gsq,karg 93 real(dp) :: lenrho,phase,qk,rlkcut,trace,xi0 94 real(dp) :: glkcut 95 !arrays 96 real(dp) :: cvec(3),gvec(3),gpl(3),gprimd(3,3) 97 real(dp) :: rhok(3),rhored(3),rpl(3) 98 real(dp),allocatable :: efg_g(:,:,:),efg_r(:,:,:) 99 real(dp),allocatable :: xcart(:,:) 100 101 ! ************************************************************************ 102 103 !DEBUG 104 !write(std_out,*)' make_efg_ion : enter' 105 !ENDDEBUG 106 107 ABI_ALLOCATE(efg_g,(3,3,natom)) 108 ABI_ALLOCATE(efg_r,(3,3,natom)) 109 ABI_ALLOCATE(xcart,(3,natom)) 110 efg(:,:,:) = zero ! final efg tensor 111 efg_g(:,:,:) = zero ! part of tensor accumulated in G space 112 efg_r(:,:,:) = zero ! part of tensor accumulated in R space 113 114 call xred2xcart(natom,rprimd,xcart,xred) ! get atomic locations in cartesian coords 115 116 do ii = 1, 3 ! generate the lengths of the unit cell edges in atomic units 117 rpl(ii) = sqrt(rprimd(1,ii)**2+rprimd(2,ii)**2+rprimd(3,ii)**2) 118 end do 119 xi0 = sqrt(pi/(maxval(rpl)*minval(rpl))) ! this estimate for xi0 is from Honma's paper 120 121 call matr3inv(rprimd,gprimd) ! gprimd holds the inverse transpose of rprimd 122 !remember ordering: rprimd( (x_comp,y_comp,z_comp), (edge 1, edge 2, edge 3) ) 123 !while gprimd( (edge 1, edge 2, edge 3),(x_comp, y_comp, z_comp) ) 124 do ii = 1, 3 ! generate the lengths of the reciprocal cell edges 125 gpl(ii) = sqrt(gprimd(ii,1)**2+gprimd(ii,2)**2+gprimd(ii,3)**2) 126 end do 127 128 !go out enough shells such that g**2/4*xi0**2 is of order 30 129 nshell = int(anint(sqrt(30.0)*xi0/(pi*minval(gpl)))) 130 glkcut = (0.95*nshell*two*pi*minval(gpl))**2 131 132 do ishell = 0, nshell ! loop over shells 133 do sx = -ishell, ishell 134 do sy = -ishell, ishell 135 do sz = -ishell, ishell 136 if ( .not. (sx==0 .and. sy==0 .and. sz==0) ) then ! avoid origin 137 ! constrain to be on shell surface, not interior 138 if ( abs(sx)==ishell .or. abs(sy)==ishell .or. abs(sz)==ishell ) then 139 cvec(1)=sx;cvec(2)=sy;cvec(3)=sz 140 ! make the g vector in cartesian coords 141 gvec(:) = zero 142 do ii = 1, 3 143 do jj = 1, 3 144 gvec(ii) = gvec(ii) + gprimd(ii,jj)*cvec(jj)*two*pi 145 end do 146 end do 147 gsq = dot_product(gvec,gvec) 148 if(gsq < glkcut) then 149 dampfac = exp(-gsq/(4.0*xi0*xi0)) ! see Honma eq. 4.8 150 do iatom = 1, natom 151 do jatom = 1, natom 152 qk = zion(typat(jatom)) ! charge on neighbor atom 153 rhok = xcart(:,jatom)-xcart(:,iatom) 154 phase = dot_product(gvec,rhok) 155 cph = cos(phase) 156 do ii = 1, 3 157 do jj = 1, 3 158 derivs = -3.0*gvec(ii)*gvec(jj)/gsq 159 if (ii == jj) derivs = 1.0 + derivs 160 efg_g(ii,jj,iatom) = efg_g(ii,jj,iatom) + & 161 & qk*cph*derivs*dampfac 162 end do ! end loop over jj 163 end do ! end loop over ii 164 end do ! end loop over jatom 165 end do ! end loop over iatom 166 end if ! constrain to gsq < glkcut 167 end if ! end selection on shell edge 168 end if ! end avoidance of origin 169 end do ! end loop over sz 170 end do ! end loop over sy 171 end do ! end loop over sx 172 end do ! end loop over ishell 173 174 !sum in real space begins here 175 176 !go out enough shells such that (r*xi0)**2 is of order 30 177 nshell = int(anint(sqrt(30.)/(minval(rpl)*xi0))) 178 rlkcut = nshell*minval(rpl)*0.95 179 ! 180 !go out enough shells so that rlkcut is of order 30 bohr 181 !nshell=int(anint(30.0/minval(rpl))) 182 !rlkcut = 0.95*nshell*minval(rpl) 183 184 do ishell = 0, nshell ! total set of cells to loop over 185 do sx = -ishell, ishell ! loop over all cells in each dimension 186 do sy = -ishell, ishell 187 do sz = -ishell, ishell 188 ! constrain to shell surface, not interior 189 if ( abs(sx)==ishell .or. abs(sy)==ishell .or. abs(sz)==ishell ) then 190 do jatom = 1, natom ! loop over atoms in shell cell 191 do iatom = 1, natom ! loop over atoms in central unit cell 192 if (.NOT. (jatom == iatom .AND. sx == 0 .AND. sy == 0 .AND. sz == 0)) then ! avoid self term 193 qk = zion(typat(jatom)) ! charge on each neighbor atom 194 ! ! rhored is the vector in crystal coords from neighbor to target 195 rhored(1) = xred(1,jatom) + sx - xred(1,iatom) 196 rhored(2) = xred(2,jatom) + sy - xred(2,iatom) 197 rhored(3) = xred(3,jatom) + sz - xred(3,iatom) 198 ! ! rhok is rhored in cartesian coords 199 rhok(1) = rprimd(1,1)*rhored(1)+rprimd(1,2)*rhored(2)+rprimd(1,3)*rhored(3) 200 rhok(2) = rprimd(2,1)*rhored(1)+rprimd(2,2)*rhored(2)+rprimd(2,3)*rhored(3) 201 rhok(3) = rprimd(3,1)*rhored(1)+rprimd(3,2)*rhored(2)+rprimd(3,3)*rhored(3) 202 trace = dot_product(rhok,rhok) 203 lenrho = sqrt(trace) 204 if (lenrho < rlkcut) then ! this restriction is critical as it ensures 205 ! ! that we sum over a sphere of atoms in real space 206 ! ! no matter what shape the unit cell has 207 karg = xi0*lenrho 208 derfc_karg = abi_derfc(karg) 209 ! see Honma eq. 2.10 for derivation of the following damping factor 210 dampfac = (1.0+3.0/(2.0*karg*karg))*exp(-karg*karg)+3.0*sqrt(pi)*derfc_karg/(4.0*karg**3) 211 do ii = 1, 3 ! loop over tensor elements 212 do jj = 1, 3 ! loop over tensor elements 213 derivs = -3.0*rhok(ii)*rhok(jj)/trace 214 if(ii == jj) derivs = derivs + 1.0 ! see Honma eq 4.8 re: sign 215 ! accumulate real space tensor element, 216 ! weighted by charge of neighbor and Ewald damping factor 217 efg_r(ii,jj,iatom) = efg_r(ii,jj,iatom) + qk*derivs*dampfac 218 end do ! end loop over jj in efg(ii,jj,iatom) 219 end do ! end loop over ii in efg(ii,jj,iatom) 220 end if ! end if statement restricting to a sphere of radius rlkcut 221 end if ! end if statement avoiding the self atom term 222 end do ! end loop over i atoms in cell 223 end do ! end loop over j atoms in cell 224 end if ! end selection on outer shell of cells only 225 end do ! end loop over sz cells 226 end do ! end loop over sy cells 227 end do ! end loop over sx cells 228 end do ! end loop over shells 229 230 !now combine the g-space and r-space parts, properly weighted (see Honma) 231 do iatom = 1, natom 232 do ii = 1, 3 233 do jj = 1, 3 234 efg(ii,jj,iatom) = four_pi*efg_g(ii,jj,iatom)/(three*ucvol)-& 235 & four*xi0**3*efg_r(ii,jj,iatom)/(three*sqrt(pi)) 236 ! note extra factor of two: compare Honma eq. 4.6 237 end do 238 end do 239 end do 240 241 ! symmetrize tensor at each atomic site using point symmetry operations 242 do iatom = 1, natom 243 call matpointsym(iatom,efg(:,:,iatom),natom,nsym,rprimd,symrel,tnons,xred) 244 end do 245 246 ABI_DEALLOCATE(efg_g) 247 ABI_DEALLOCATE(efg_r) 248 ABI_DEALLOCATE(xcart) 249 250 !DEBUG 251 !write(std_out,*)' make_efg_ion : exit ' 252 !stop 253 !ENDDEBUG 254 255 end subroutine make_efg_ion