TABLE OF CONTENTS
ABINIT/calc_efg [ Functions ]
NAME
calc_efg
FUNCTION
calculation and output of electric field gradient tensor at each atomic site
INPUTS
mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor natom=number of atoms in cell. nfft=number of points on fft grid ngfft(18)=details of fft nhat(nfft,nspden)=compensation charge density nspden=number of spin densities nsym=number of symmetries in space group ntypat=number of atom types nucefg=1 to print summary output, 2 for detailed output ptcharge(ntypat)=user input charges on atoms to make simple point charge calc paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh pawang <type(pawang_type)>=paw angular mesh and related data pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data quadmom(ntypat)=quadrupole moments in barns of different atomic nuclei rhor(nfft,nspden)=electron density on grid (strictly $\tilde{n}+\hat{n}$) rprimd(3,3)=matrix relating cartesian coordinates to crystal coordinates symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations tnons(3,nsym)=nonsymmorphic translations typat(natom)=type (integer) for each atom ucvol=unit cell volume in Bohr^3 usepaw=1 if we are using PAW formalism, 0 else xred(3,natom)=vectors locating each atom in the unit cell, in crystal coords zion(ntypat)=net core charge on each type of atom
OUTPUT
(only writing, printing)
SOURCE
100 subroutine calc_efg(mpi_enreg,my_natom,natom,nfft,ngfft,nhat,nspden,nsym,nucefg,ntypat,& 101 paw_an,pawang,pawrad,pawrhoij,pawtab,& 102 ptcharge,quadmom,rhor,rprimd,symrel,tnons,typat,ucvol,usepaw,xred,zion,& 103 mpi_atmtab,comm_atom) ! optional arguments (parallelism) 104 105 !Arguments ------------------------------------ 106 !scalars 107 integer,intent(in) :: my_natom,natom,nfft,nspden,nsym,nucefg,ntypat,usepaw 108 integer,optional,intent(in) :: comm_atom 109 real(dp),intent(in) :: ucvol 110 type(MPI_type),intent(in) :: mpi_enreg 111 type(pawang_type),intent(in) :: pawang 112 !arrays 113 integer,intent(in) :: ngfft(18),symrel(3,3,nsym),typat(natom) 114 integer,optional,target,intent(in) :: mpi_atmtab(:) 115 real(dp),intent(in) :: nhat(nfft,nspden),ptcharge(ntypat) 116 real(dp),intent(in) :: quadmom(ntypat),rhor(nfft,nspden),rprimd(3,3) 117 real(dp),intent(in) :: tnons(3,nsym),zion(ntypat) 118 real(dp),intent(inout) :: xred(3,natom) 119 type(paw_an_type),intent(in) :: paw_an(my_natom) 120 type(pawrad_type),intent(in) :: pawrad(ntypat) 121 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 122 type(pawtab_type),intent(in) :: pawtab(ntypat) 123 124 !Local variables------------------------------- 125 !scalars 126 integer :: ii,INFO,LDA,LWORK,N,iatom,my_comm_atom 127 logical :: my_atmtab_allocated,paral_atom 128 real(dp),parameter :: efg_si = 9.7173624292E21 ! convert EFG in au to SI units V/m2 129 real(dp) :: cq,eta,vxx,vyy,vzz 130 character(len=500) :: message 131 !arrays 132 integer,pointer :: my_atmtab(:) 133 real(dp) :: eigval(3),matr(3,3),work(8) 134 real(dp),allocatable :: efg(:,:,:),efg_el(:,:,:),efg_ion(:,:,:),efg_paw(:,:,:) 135 real(dp),allocatable :: efg_point_charge(:,:,:) 136 137 ! ************************************************************************ 138 139 !Compatibility tests 140 if (usepaw /= 1) then 141 message = ' usepaw /= 1 but EFG calculation requires PAW ' 142 ABI_ERROR(message) 143 end if 144 145 !Set up parallelism over atoms 146 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 147 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 148 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 149 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 150 151 ABI_MALLOC(efg,(3,3,natom)) 152 ABI_MALLOC(efg_el,(3,3,natom)) 153 ABI_MALLOC(efg_ion,(3,3,natom)) 154 ABI_MALLOC(efg_paw,(3,3,natom)) 155 ABI_MALLOC(efg_point_charge,(3,3,natom)) 156 efg_el(:,:,:) = zero 157 efg_ion(:,:,:) = zero 158 efg_paw(:,:,:) = zero 159 efg_point_charge(:,:,:) = zero 160 161 call make_efg_el(efg_el,mpi_enreg,natom,nfft,ngfft,nhat,nspden,nsym,rhor,rprimd,symrel,tnons,xred) 162 163 call make_efg_ion(efg_ion,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,zion) 164 165 if (paral_atom) then 166 call make_efg_onsite(efg_paw,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab,& 167 & rprimd,symrel,tnons,xred,comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 168 else 169 call make_efg_onsite(efg_paw,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab,& 170 & rprimd,symrel,tnons,xred) 171 end if 172 173 !calculate efg due to pure point charges, as input in variable ptcharge(ntypat) 174 !note here all atoms of the same type will have the same valence; in the future this 175 !could be made more flexible by having ptcharge(natom) but that will require a slightly 176 !different version than the existing make_efg_ion routine 177 if(nucefg > 2) then 178 call make_efg_ion(efg_point_charge,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,ptcharge) 179 end if 180 181 efg(:,:,:) = efg_el(:,:,:) + efg_ion(:,:,:) + efg_paw(:,:,:) 182 183 write(message,'(a,a,a)' ) ch10,' Electric Field Gradient Calculation ',ch10 184 call wrtout(ab_out,message,'COLL') 185 186 LDA=3; LWORK=8;N=3 ! these parameters are needed for the LAPACK dsyev routine 187 do iatom = 1, natom 188 matr(:,:) = efg(:,:,iatom) 189 call dsyev('V','U',N,matr,LDA,eigval,work,LWORK,INFO) ! get eigenvalues and eigenvectors 190 if (eigval(3) > abs(eigval(1)) ) then ! In NMR, the convention is that whatever component is 191 ! largest in magnitude is called Vzz, next comes Vxx, then Vyy 192 vzz = eigval(3) 193 vxx = eigval(1) 194 vyy = eigval(2) 195 else 196 vzz = eigval(1) 197 vxx = eigval(3) 198 vyy = eigval(2) 199 end if 200 ! Cq = (eq)*(eQ)/h where eq is the efg vzz; Q is the nuclear quad moment in barns (10E-28 m2) 201 ! multiply vzz (in au) by efg_si to get volts/m^2 202 ! multiply quadmom by e_Cb (the electric charge unit) and 1D-28 to get eQ in SI units 203 ! divide by Plancks Constant to get freq 204 ! multiply by 1D-6 to get MHz 205 cq = 1.0D-6*vzz*efg_si*quadmom(typat(iatom))*e_Cb*1.0D-28/6.62607015D-34 206 if (abs(cq) > tol8) then 207 eta = abs(vxx-vyy)/abs(vzz) 208 else 209 cq = zero 210 eta = zero ! if Cq is small then eta is meaningless 211 end if 212 213 write(message,'(a,a,i4,a,i4)')ch10,' atom : ',iatom,' typat : ',typat(iatom) 214 call wrtout(ab_out,message,'COLL') 215 if (nucefg > 1) then 216 write(message,'(2a,f9.4,a,f9.4,a,f9.4)') ch10,' Nuclear quad. mom. (barns) : ',quadmom(typat(iatom)),& 217 & ' Cq (MHz) : ',cq,' eta : ',eta 218 call wrtout(ab_out,message,'COLL') 219 end if 220 221 ! for printing and test portability, it's better to simply set very small eigvals to zero 222 do ii=1,3 223 if (abs(eigval(ii))<tol8) eigval(ii)=zero 224 end do 225 write(message,'(2a,f13.6,a,es16.8,a,a,3f13.6)')ch10,' efg eigval (au) : ',eigval(1),' ; (V/m^2) : ',eigval(1)*efg_si,ch10,& 226 & '- eigvec : ',matr(1,1),matr(2,1),matr(3,1) 227 call wrtout(ab_out,message,'COLL') 228 write(message,'(2a,f13.6,a,es16.8,a,a,3f13.6)')ch10,' efg eigval (au) : ',eigval(2),' ; (V/m^2) : ',eigval(2)*efg_si,ch10,& 229 & '- eigvec : ',matr(1,2),matr(2,2),matr(3,2) 230 call wrtout(ab_out,message,'COLL') 231 write(message,'(2a,f13.6,a,es16.8,a,a,3f13.6)')ch10,' efg eigval (au) : ',eigval(3),' ; (V/m^2) : ',eigval(3)*efg_si,ch10,& 232 & '- eigvec : ',matr(1,3),matr(2,3),matr(3,3) 233 call wrtout(ab_out,message,'COLL') 234 write(message,'(a,a,3f13.6)')ch10,' total efg : ',efg(1,1,iatom),efg(1,2,iatom),efg(1,3,iatom) 235 call wrtout(ab_out,message,'COLL') 236 write(message,'(a,3f13.6)')' total efg : ',efg(2,1,iatom),efg(2,2,iatom),efg(2,3,iatom) 237 call wrtout(ab_out,message,'COLL') 238 write(message,'(a,3f13.6,a)')' total efg : ',efg(3,1,iatom),efg(3,2,iatom),efg(3,3,iatom),ch10 239 call wrtout(ab_out,message,'COLL') 240 write(message,'(a,a,3f13.6)')ch10,' efg_el : ',efg_el(1,1,iatom),efg_el(1,2,iatom),efg_el(1,3,iatom) 241 call wrtout(ab_out,message,'COLL') 242 write(message,'(a,3f13.6)')' efg_el : ',efg_el(2,1,iatom),efg_el(2,2,iatom),efg_el(2,3,iatom) 243 call wrtout(ab_out,message,'COLL') 244 write(message,'(a,3f13.6,a)')' efg_el : ',efg_el(3,1,iatom),efg_el(3,2,iatom),efg_el(3,3,iatom),ch10 245 call wrtout(ab_out,message,'COLL') 246 write(message,'(a,3f13.6)')' efg_ion : ',efg_ion(1,1,iatom),efg_ion(1,2,iatom),efg_ion(1,3,iatom) 247 call wrtout(ab_out,message,'COLL') 248 write(message,'(a,3f13.6)')' efg_ion : ',efg_ion(2,1,iatom),efg_ion(2,2,iatom),efg_ion(2,3,iatom) 249 call wrtout(ab_out,message,'COLL') 250 write(message,'(a,3f13.6,a)')' efg_ion : ',efg_ion(3,1,iatom),efg_ion(3,2,iatom),efg_ion(3,3,iatom),ch10 251 call wrtout(ab_out,message,'COLL') 252 write(message,'(a,3f13.6)')' efg_paw : ',efg_paw(1,1,iatom),efg_paw(1,2,iatom),efg_paw(1,3,iatom) 253 call wrtout(ab_out,message,'COLL') 254 write(message,'(a,3f13.6)')' efg_paw : ',efg_paw(2,1,iatom),efg_paw(2,2,iatom),efg_paw(2,3,iatom) 255 call wrtout(ab_out,message,'COLL') 256 write(message,'(a,3f13.6,a)')' efg_paw : ',efg_paw(3,1,iatom),efg_paw(3,2,iatom),efg_paw(3,3,iatom),ch10 257 call wrtout(ab_out,message,'COLL') 258 if (nucefg > 2) then ! write output of pure pointcharge calculation 259 matr(:,:) = efg_point_charge(:,:,iatom) 260 call dsyev('V','U',N,matr,LDA,eigval,work,LWORK,INFO) ! get eigenvalues and eigenvectors 261 if (eigval(3) > abs(eigval(1)) ) then ! In NMR, the convention is that whatever component is 262 ! largest in magnitude is called Vzz, next comes Vxx, then Vyy 263 vzz = eigval(3) 264 vxx = eigval(1) 265 vyy = eigval(2) 266 else 267 vzz = eigval(1) 268 vxx = eigval(3) 269 vyy = eigval(2) 270 end if 271 cq = 1.0D-6*vzz*efg_si*quadmom(typat(iatom))*e_Cb*1.0D-28/6.62607015D-34 272 if (abs(cq) > tol8) then 273 eta = abs(vxx-vyy)/abs(vzz) 274 else 275 eta = zero ! if Cq is small then eta is meaningless 276 end if 277 write(message,'(a,f9.4,a,f9.4)') ' Point charge Cq = ',cq,' MHz eta = ',eta 278 call wrtout(ab_out,message,'COLL') 279 write(message,'(2a,f13.6,a,es16.8,a,a,3f13.6)')ch10,' point charge eigval (au) : ',eigval(1),' ; (V/m^2) : ',eigval(1)*efg_si,ch10,& 280 & '- eigvec : ',matr(1,1),matr(2,1),matr(3,1) 281 call wrtout(ab_out,message,'COLL') 282 write(message,'(2a,f13.6,a,es16.8,a,a,3f13.6)')ch10,' point charge eigval (au) : ',eigval(2),' ; (V/m^2) : ',eigval(2)*efg_si,ch10,& 283 & '- eigvec : ',matr(1,2),matr(2,2),matr(3,2) 284 call wrtout(ab_out,message,'COLL') 285 write(message,'(2a,f13.6,a,es16.8,a,a,3f13.6)')ch10,' point charge eigval (au) : ',eigval(3),' ; (V/m^2) : ',eigval(3)*efg_si,ch10,& 286 & '- eigvec : ',matr(1,3),matr(2,3),matr(3,3) 287 call wrtout(ab_out,message,'COLL') 288 write(message,'(a,a,3f13.6)')ch10,' point charge efg : ',efg_point_charge(1,1,iatom),& 289 & efg_point_charge(1,2,iatom),efg_point_charge(1,3,iatom) 290 call wrtout(ab_out,message,'COLL') 291 write(message,'(a,3f13.6)')' point charge efg : ',efg_point_charge(2,1,iatom),& 292 & efg_point_charge(2,2,iatom),efg_point_charge(2,3,iatom) 293 call wrtout(ab_out,message,'COLL') 294 write(message,'(a,3f13.6,a)')' point charge efg : ',efg_point_charge(3,1,iatom),& 295 & efg_point_charge(3,2,iatom),efg_point_charge(3,3,iatom),ch10 296 call wrtout(ab_out,message,'COLL') 297 end if 298 end do 299 write(message,'(3a)')ch10,ch10,ch10 300 call wrtout(ab_out,message,'COLL') 301 302 ABI_FREE(efg) 303 ABI_FREE(efg_el) 304 ABI_FREE(efg_ion) 305 ABI_FREE(efg_paw) 306 ABI_FREE(efg_point_charge) 307 308 !Destroy atom table used for parallelism 309 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 310 311 !DEBUG 312 !write(std_out,*)' calc_efg : exit ' 313 !stop 314 !ENDDEBUG 315 316 end subroutine calc_efg
ABINIT/calc_fc [ Functions ]
NAME
calc_fc
FUNCTION
calculation and output of Fermi-contact term at each atomic site
INPUTS
mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms my_natom=number of atoms treated by current processor natom=number of atoms in cell. nspden=number of spin density components ntypat=number of atom types pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data typat(natom)=type (integer) for each atom usepaw=1 if PAW is activated
OUTPUT
(only writing, printing)
SIDE EFFECTS
NOTES
SOURCE
350 subroutine calc_fc(my_natom,natom,nspden,ntypat,pawrad,pawrhoij,pawtab,typat,usepaw,& 351 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 352 353 !Arguments ------------------------------------ 354 !scalars 355 integer,intent(in) :: my_natom,natom,nspden,ntypat,usepaw 356 integer,optional,intent(in) :: comm_atom 357 !arrays 358 integer,intent(in) :: typat(natom) 359 integer,optional,target,intent(in) :: mpi_atmtab(:) 360 type(pawrad_type),intent(in) :: pawrad(ntypat) 361 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 362 type(pawtab_type),intent(in) :: pawtab(ntypat) 363 364 !Local variables------------------------------- 365 !scalars 366 integer :: iatom,my_comm_atom 367 logical :: my_atmtab_allocated,paral_atom 368 character(len=500) :: message 369 !arrays 370 integer,pointer :: my_atmtab(:) 371 real(dp),allocatable :: fc(:,:) 372 373 !*********************************************************************** 374 375 !Compatibility tests 376 if (usepaw /= 1) then 377 message = ' usepaw /= 1 but Fermi-contact calculation requires PAW ' 378 ABI_ERROR(message) 379 end if 380 381 !Set up parallelism over atoms 382 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 383 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 384 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 385 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 386 387 !Initialization 388 ABI_MALLOC(fc,(nspden,natom)) 389 390 !Computation 391 if (paral_atom) then 392 call make_fc_paw(fc,my_natom,natom,nspden,ntypat,pawrhoij,pawrad,pawtab,& 393 & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 394 else 395 call make_fc_paw(fc,my_natom,natom,nspden,ntypat,pawrhoij,pawrad,pawtab) 396 end if 397 398 !Printing 399 write(message,'(a,a,a)' ) ch10,' Fermi-contact Term Calculation ',ch10 400 call wrtout(ab_out,message,'COLL') 401 402 do iatom = 1, natom 403 if (nspden == 2) then 404 write(message,'(a,i3,a,i3,a,f12.4)') ' Atom ',iatom,', typat ',typat(iatom),': FC total = ',& 405 & fc(1,iatom)+fc(2,iatom) 406 call wrtout(ab_out,message,'COLL') 407 write(message,'(a,i3,a,i3,a,f12.4)') ' Atom ',iatom,', typat ',typat(iatom),': FC up - down = ',& 408 & fc(1,iatom)-fc(2,iatom) 409 call wrtout(ab_out,message,'COLL') 410 else 411 write(message,'(a,i3,a,i3,a,f12.4)') ' Atom ',iatom,', typat ',typat(iatom),': FC = ',& 412 & fc(1,iatom) 413 call wrtout(ab_out,message,'COLL') 414 end if 415 end do 416 417 write(message,'(3a)')ch10,ch10,ch10 418 call wrtout(ab_out,message,'COLL') 419 420 !Memory deallocation 421 ABI_FREE(fc) 422 423 !Destroy atom table used for parallelism 424 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 425 426 end subroutine calc_fc
ABINIT/m_nucprop [ Modules ]
NAME
m_nucprop
FUNCTION
routines used to compute properties at the nuclear sites, including electric field gradient and Fermi contact
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (MT, JWZ) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 module m_nucprop 24 25 use defs_basis 26 use m_abicore 27 use m_errors 28 29 use defs_abitypes, only : MPI_type 30 use m_mpinfo, only : ptabs_fourdp 31 use m_xmpi, only : xmpi_comm_self, xmpi_sum 32 use m_geometry, only : xred2xcart 33 use m_linalg_interfaces, only: dsyev 34 use m_paw_an, only : paw_an_type 35 use m_pawang, only : pawang_type 36 use m_pawrad, only : pawrad_type 37 use m_pawtab, only : pawtab_type 38 use m_pawrhoij, only : pawrhoij_type 39 use m_paw_nmr, only : make_efg_onsite,make_fc_paw 40 use m_paral_atom, only : get_my_atmtab,free_my_atmtab 41 use m_special_funcs, only : abi_derfc 42 use m_symtk, only : matr3inv, matpointsym 43 use m_fft, only : fourdp 44 45 implicit none 46 47 private
ABINIT/make_efg_el [ Functions ]
NAME
make_efg_el
FUNCTION
compute the electric field gradient due to electron density
INPUTS
mpi_enreg=information about MPI parallelization natom, number of atoms in unit cell nfft,ngfft(18), number of FFT points and details of FFT nhat(nfft,nspden) compensation charge density nspden, number of spin components nsym=number of symmetries in space group rhor(nfft,nspden), valence electron density, here $\tilde{n} + \hat{n}$ rprimd(3,3), conversion from crystal coordinates to cartesian coordinates symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations tnons(3,nsym) = nonsymmorphic translations xred(3,natom), location of atoms in crystal coordinates.
OUTPUT
efg(3,3,natom), the 3x3 efg tensor at each atomic site due to rhor
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 valence
electrons, at each atomic site in the unit cell. Key references: Kresse and Joubert, ``From
ultrasoft pseudopotentials to the projector augmented wave method'', Phys. Rev. B. 59, 1758--1775 (1999) [[cite:Kresse1999]],
and 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) [[cite:Profeta2003]]. This
routine computes the second derivatives of the potential generated by $\tilde{n}$ (see Kresse and Joubert
for notation, Fourier-transforming the density, doing the sum in G space, and then transforming back at
each atomic site. The final formula is
\begin{displaymath}
\frac{\partial^2 V}{\partial x_\alpha\partial x_\beta} = -4\pi^2\sum_G (G_\alpha G_\beta - \delta_{\alpha,\beta}G^2/3)
\left(\frac{\tilde{n}(G)}{\pi G^2}\right)e^{2\pi i G\cdot R}
\end{displaymath}
SOURCE
692 subroutine make_efg_el(efg,mpi_enreg,natom,nfft,ngfft,nhat,nspden,nsym,rhor,rprimd,symrel,tnons,xred) 693 694 !Arguments ------------------------------------ 695 !scalars 696 integer,intent(in) :: natom,nfft,nspden,nsym 697 type(MPI_type),intent(in) :: mpi_enreg 698 !arrays 699 integer,intent(in) :: ngfft(18),symrel(3,3,nsym) 700 real(dp),intent(in) :: nhat(nfft,nspden),rhor(nfft,nspden),rprimd(3,3),tnons(3,nsym),xred(3,natom) 701 real(dp),intent(out) :: efg(3,3,natom) 702 703 !Local variables------------------------------- 704 !scalars 705 integer :: cplex,fftdir,fofg_index,iatom,i1,i2,i2_local,i23,i3,id1,id2,id3 706 integer :: ierr,ig,ig2,ig3,ii,ii1,ing,jj 707 integer :: me_fft,n1,n2,n3,nproc_fft,tim_fourdp 708 real(dp) :: cph,derivs,phase,sph,trace 709 ! type(MPI_type) :: mpi_enreg_seq 710 !arrays 711 integer :: id(3) 712 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 713 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 714 real(dp) :: gprimd(3,3),gqred(3),gvec(3),ratom(3) 715 real(dp),allocatable :: fofg(:,:),fofr(:),gq(:,:),xcart(:,:) 716 717 ! ************************************************************************ 718 719 !DEBUG 720 !write(std_out,*)' make_efg_el : enter' 721 !ENDDEBUG 722 723 ABI_MALLOC(fofg,(2,nfft)) 724 ABI_MALLOC(fofr,(nfft)) 725 ABI_MALLOC(xcart,(3,natom)) 726 727 efg(:,:,:) = zero 728 call xred2xcart(natom,rprimd,xcart,xred) ! get atomic locations in cartesian coords 729 call matr3inv(rprimd,gprimd) 730 731 tim_fourdp = 0 ! timing code, not using 732 fftdir = -1 ! FT from R to G 733 cplex = 1 ! fofr is real 734 !here we are only interested in the valence pseudo charge density, which is rhor(:,1)-nhat(:,1) 735 !regardless of the value of nspden. This may change in the future depending on 736 !developments with noncollinear magnetization and so forth. Such a change will 737 !require an additional loop over nspden. 738 !Multiply by -1 to convert the electron particle density to the charge density 739 fofr(:) = -(rhor(:,1)-nhat(:,1)) 740 741 ! Get the distrib associated with this fft_grid See hartre.F90 for another example where 742 ! this is done 743 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 744 nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft 745 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 746 747 call fourdp(cplex,fofg,fofr,fftdir,mpi_enreg,nfft,1,ngfft,tim_fourdp) ! construct charge density in G space 748 749 ! the following loops over G vectors has been copied from hartre.F90 in order to be compatible with 750 ! possible FFT parallelism 751 752 ! In order to speed the routine, precompute the components of g 753 ! Also check if the booked space was large enough... 754 ABI_MALLOC(gq,(3,max(n1,n2,n3))) 755 do ii=1,3 756 id(ii)=ngfft(ii)/2+2 757 do ing=1,ngfft(ii) 758 ig=ing-(ing/id(ii))*ngfft(ii)-1 759 gq(ii,ing)=ig 760 end do 761 end do 762 id1=n1/2+2;id2=n2/2+2;id3=n3/2+2 763 764 ! Triple loop on each dimension 765 do i3=1,n3 766 ig3=i3-(i3/id3)*n3-1 767 gqred(3) = gq(3,i3) 768 769 do i2=1,n2 770 ig2=i2-(i2/id2)*n2-1 771 if (fftn2_distrib(i2) == me_fft) then 772 773 gqred(2) = gq(2,i2) 774 i2_local = ffti2_local(i2) 775 i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1)) 776 ! Do the test that eliminates the Gamma point outside of the inner loop 777 ii1=1 778 if(i23==0 .and. ig2==0 .and. ig3==0) ii1=2 779 780 ! Final inner loop on the first dimension (note the lower limit) 781 do i1=ii1,n1 782 ! gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 783 gqred(1) = gq(1,i1) 784 gvec(1:3) = MATMUL(gprimd,gqred) 785 fofg_index=i1+i23 786 trace = dot_product(gvec,gvec) 787 do ii = 1, 3 ! sum over components of efg tensor 788 do jj = 1, 3 ! sum over components of efg tensor 789 derivs = gvec(ii)*gvec(jj) ! This term is $G_\alpha G_\beta$ 790 if (ii == jj) derivs = derivs - trace/three 791 do iatom = 1, natom ! sum over atoms in unit cell 792 ratom(:) = xcart(:,iatom) ! extract location of atom iatom 793 phase = two_pi*dot_product(gvec,ratom) ! argument of $e^{2\pi i G\cdot R}$ 794 cph = cos(phase) 795 sph = sin(phase) 796 efg(ii,jj,iatom) = efg(ii,jj,iatom) - & 797 & four_pi*derivs*(fofg(1,fofg_index)*cph-fofg(2,fofg_index)*sph)/trace ! real part of efg tensor 798 end do ! end loop over atoms in cell 799 end do ! end loop over jj in V_ij 800 end do ! end loop over ii in V_ij 801 end do ! End loop on i1 802 end if 803 end do ! End loop on i2 804 end do ! End loop on i3 805 806 call xmpi_sum(efg,mpi_enreg%comm_fft,ierr) 807 808 ! symmetrize tensor at each atomic site using point symmetry operations 809 do iatom = 1, natom 810 call matpointsym(iatom,efg(:,:,iatom),natom,nsym,rprimd,symrel,tnons,xred) 811 end do 812 813 ABI_FREE(fofg) 814 ABI_FREE(fofr) 815 ABI_FREE(xcart) 816 ABI_FREE(gq) 817 818 !DEBUG 819 !write(std_out,*)' make_efg_el : exit ' 820 !stop 821 !ENDDEBUG 822 823 end subroutine make_efg_el
ABINIT/make_efg_ion [ Functions ]
NAME
make_efg_ion
FUNCTION
compute the electric field gradient due to ionic cores
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) [[cite:Profeta2003]];
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) [[cite:Honma1977]];
and Kresse and Joubert, ``From ultrasoft pseudopotentials to the projector augmented wave method'',
Phys. Rev. B. 59, 1758--1775 (1999) [[cite:Kresse1999]]. 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.
SOURCE
470 subroutine make_efg_ion(efg,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,zion) 471 472 !Arguments ------------------------------------ 473 !scalars 474 integer,intent(in) :: natom,nsym,ntypat 475 real(dp) :: ucvol 476 !arrays 477 integer,intent(in) :: symrel(3,3,nsym),typat(natom) 478 real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym) 479 real(dp),intent(in) :: zion(ntypat) 480 real(dp),intent(inout) :: xred(3,natom) 481 real(dp),intent(out) :: efg(3,3,natom) 482 !Local variables------------------------------- 483 !scalars 484 integer :: iatom,ishell,ii,jatom,jj,nshell,sx,sy,sz 485 real(dp) :: cph,dampfac,derfc_karg,derivs,gsq,karg 486 real(dp) :: lenrho,phase,qk,rlkcut,trace,xi0 487 real(dp) :: glkcut 488 !arrays 489 real(dp) :: cvec(3),gvec(3),gpl(3),gprimd(3,3) 490 real(dp) :: rhok(3),rhored(3),rpl(3) 491 real(dp),allocatable :: efg_g(:,:,:),efg_r(:,:,:) 492 real(dp),allocatable :: xcart(:,:) 493 494 ! ************************************************************************ 495 496 !DEBUG 497 !write(std_out,*)' make_efg_ion : enter' 498 !ENDDEBUG 499 500 ABI_MALLOC(efg_g,(3,3,natom)) 501 ABI_MALLOC(efg_r,(3,3,natom)) 502 ABI_MALLOC(xcart,(3,natom)) 503 efg(:,:,:) = zero ! final efg tensor 504 efg_g(:,:,:) = zero ! part of tensor accumulated in G space 505 efg_r(:,:,:) = zero ! part of tensor accumulated in R space 506 507 call xred2xcart(natom,rprimd,xcart,xred) ! get atomic locations in cartesian coords 508 509 do ii = 1, 3 ! generate the lengths of the unit cell edges in atomic units 510 rpl(ii) = sqrt(rprimd(1,ii)**2+rprimd(2,ii)**2+rprimd(3,ii)**2) 511 end do 512 xi0 = sqrt(pi/(maxval(rpl)*minval(rpl))) ! this estimate for xi0 is from Honma's paper 513 514 call matr3inv(rprimd,gprimd) ! gprimd holds the inverse transpose of rprimd 515 !remember ordering: rprimd( (x_comp,y_comp,z_comp), (edge 1, edge 2, edge 3) ) 516 !while gprimd( (edge 1, edge 2, edge 3),(x_comp, y_comp, z_comp) ) 517 do ii = 1, 3 ! generate the lengths of the reciprocal cell edges 518 gpl(ii) = sqrt(gprimd(ii,1)**2+gprimd(ii,2)**2+gprimd(ii,3)**2) 519 end do 520 521 !go out enough shells such that g**2/4*xi0**2 is of order 30 522 nshell = int(anint(sqrt(30.0)*xi0/(pi*minval(gpl)))) 523 glkcut = (0.95*nshell*two*pi*minval(gpl))**2 524 525 do ishell = 0, nshell ! loop over shells 526 do sx = -ishell, ishell 527 do sy = -ishell, ishell 528 do sz = -ishell, ishell 529 if ( .not. (sx==0 .and. sy==0 .and. sz==0) ) then ! avoid origin 530 ! constrain to be on shell surface, not interior 531 if ( abs(sx)==ishell .or. abs(sy)==ishell .or. abs(sz)==ishell ) then 532 cvec(1)=sx;cvec(2)=sy;cvec(3)=sz 533 ! make the g vector in cartesian coords 534 gvec(:) = zero 535 do ii = 1, 3 536 do jj = 1, 3 537 gvec(ii) = gvec(ii) + gprimd(ii,jj)*cvec(jj)*two*pi 538 end do 539 end do 540 gsq = dot_product(gvec,gvec) 541 if(gsq < glkcut) then 542 dampfac = exp(-gsq/(4.0*xi0*xi0)) ! see Honma eq. 4.8 543 do iatom = 1, natom 544 do jatom = 1, natom 545 qk = zion(typat(jatom)) ! charge on neighbor atom 546 rhok = xcart(:,jatom)-xcart(:,iatom) 547 phase = dot_product(gvec,rhok) 548 cph = cos(phase) 549 do ii = 1, 3 550 do jj = 1, 3 551 derivs = -3.0*gvec(ii)*gvec(jj)/gsq 552 if (ii == jj) derivs = 1.0 + derivs 553 efg_g(ii,jj,iatom) = efg_g(ii,jj,iatom) + & 554 & qk*cph*derivs*dampfac 555 end do ! end loop over jj 556 end do ! end loop over ii 557 end do ! end loop over jatom 558 end do ! end loop over iatom 559 end if ! constrain to gsq < glkcut 560 end if ! end selection on shell edge 561 end if ! end avoidance of origin 562 end do ! end loop over sz 563 end do ! end loop over sy 564 end do ! end loop over sx 565 end do ! end loop over ishell 566 567 !sum in real space begins here 568 569 !go out enough shells such that (r*xi0)**2 is of order 30 570 nshell = int(anint(sqrt(30.)/(minval(rpl)*xi0))) 571 rlkcut = nshell*minval(rpl)*0.95 572 ! 573 !go out enough shells so that rlkcut is of order 30 bohr 574 !nshell=int(anint(30.0/minval(rpl))) 575 !rlkcut = 0.95*nshell*minval(rpl) 576 577 do ishell = 0, nshell ! total set of cells to loop over 578 do sx = -ishell, ishell ! loop over all cells in each dimension 579 do sy = -ishell, ishell 580 do sz = -ishell, ishell 581 ! constrain to shell surface, not interior 582 if ( abs(sx)==ishell .or. abs(sy)==ishell .or. abs(sz)==ishell ) then 583 do jatom = 1, natom ! loop over atoms in shell cell 584 do iatom = 1, natom ! loop over atoms in central unit cell 585 if (.NOT. (jatom == iatom .AND. sx == 0 .AND. sy == 0 .AND. sz == 0)) then ! avoid self term 586 qk = zion(typat(jatom)) ! charge on each neighbor atom 587 ! ! rhored is the vector in crystal coords from neighbor to target 588 rhored(1) = xred(1,jatom) + sx - xred(1,iatom) 589 rhored(2) = xred(2,jatom) + sy - xred(2,iatom) 590 rhored(3) = xred(3,jatom) + sz - xred(3,iatom) 591 ! ! rhok is rhored in cartesian coords 592 rhok(1) = rprimd(1,1)*rhored(1)+rprimd(1,2)*rhored(2)+rprimd(1,3)*rhored(3) 593 rhok(2) = rprimd(2,1)*rhored(1)+rprimd(2,2)*rhored(2)+rprimd(2,3)*rhored(3) 594 rhok(3) = rprimd(3,1)*rhored(1)+rprimd(3,2)*rhored(2)+rprimd(3,3)*rhored(3) 595 trace = dot_product(rhok,rhok) 596 lenrho = sqrt(trace) 597 if (lenrho < rlkcut) then ! this restriction is critical as it ensures 598 ! ! that we sum over a sphere of atoms in real space 599 ! ! no matter what shape the unit cell has 600 karg = xi0*lenrho 601 derfc_karg = abi_derfc(karg) 602 ! see Honma eq. 2.10 for derivation of the following damping factor 603 dampfac = (1.0+3.0/(2.0*karg*karg))*exp(-karg*karg)+3.0*sqrt(pi)*derfc_karg/(4.0*karg**3) 604 do ii = 1, 3 ! loop over tensor elements 605 do jj = 1, 3 ! loop over tensor elements 606 derivs = -3.0*rhok(ii)*rhok(jj)/trace 607 if(ii == jj) derivs = derivs + 1.0 ! see Honma eq 4.8 re: sign 608 ! accumulate real space tensor element, 609 ! weighted by charge of neighbor and Ewald damping factor 610 efg_r(ii,jj,iatom) = efg_r(ii,jj,iatom) + qk*derivs*dampfac 611 end do ! end loop over jj in efg(ii,jj,iatom) 612 end do ! end loop over ii in efg(ii,jj,iatom) 613 end if ! end if statement restricting to a sphere of radius rlkcut 614 end if ! end if statement avoiding the self atom term 615 end do ! end loop over i atoms in cell 616 end do ! end loop over j atoms in cell 617 end if ! end selection on outer shell of cells only 618 end do ! end loop over sz cells 619 end do ! end loop over sy cells 620 end do ! end loop over sx cells 621 end do ! end loop over shells 622 623 !now combine the g-space and r-space parts, properly weighted (see Honma) 624 do iatom = 1, natom 625 do ii = 1, 3 626 do jj = 1, 3 627 efg(ii,jj,iatom) = four_pi*efg_g(ii,jj,iatom)/(three*ucvol)-& 628 & four*xi0**3*efg_r(ii,jj,iatom)/(three*sqrt(pi)) 629 ! note extra factor of two: compare Honma eq. 4.6 630 end do 631 end do 632 end do 633 634 ! symmetrize tensor at each atomic site using point symmetry operations 635 do iatom = 1, natom 636 call matpointsym(iatom,efg(:,:,iatom),natom,nsym,rprimd,symrel,tnons,xred) 637 end do 638 639 ABI_FREE(efg_g) 640 ABI_FREE(efg_r) 641 ABI_FREE(xcart) 642 643 !DEBUG 644 !write(std_out,*)' make_efg_ion : exit ' 645 !stop 646 !ENDDEBUG 647 648 end subroutine make_efg_ion