TABLE OF CONTENTS
ABINIT/calc_efg [ Functions ]
NAME
calc_efg
FUNCTION
calculation and output of electric field gradient tensor at each atomic site
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group (JZ,MT) 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
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 nspden=number of spin densities nsym=number of symmetries in space group ntypat=number of atom types paral_kgb 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 prtefg=1 to print summary output, 2 for detailed output 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)
SIDE EFFECTS
NOTES
PARENTS
outscfcv
CHILDREN
dsyev,free_my_atmtab,get_my_atmtab,make_efg_el,make_efg_ion make_efg_onsite,wrtout
SOURCE
62 #if defined HAVE_CONFIG_H 63 #include "config.h" 64 #endif 65 66 #include "abi_common.h" 67 68 subroutine calc_efg(mpi_enreg,my_natom,natom,nfft,ngfft,nspden,nsym,ntypat,paral_kgb,& 69 & paw_an,pawang,pawrad,pawrhoij,pawtab,& 70 & ptcharge,prtefg,quadmom,rhor,rprimd,symrel,tnons,typat,ucvol,usepaw,xred,zion,& 71 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 72 73 use defs_basis 74 use defs_abitypes, only : MPI_type 75 use m_profiling_abi 76 use m_errors 77 use m_xmpi 78 79 use m_linalg_interfaces 80 81 use m_pawang, only : pawang_type 82 use m_pawrad, only : pawrad_type 83 use m_pawtab, only : pawtab_type 84 use m_paw_an, only : paw_an_type 85 use m_pawrhoij, only : pawrhoij_type 86 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 87 88 !This section has been created automatically by the script Abilint (TD). 89 !Do not modify the following lines by hand. 90 #undef ABI_FUNC 91 #define ABI_FUNC 'calc_efg' 92 use interfaces_14_hidewrite 93 use interfaces_65_paw 94 use interfaces_67_common, except_this_one => calc_efg 95 !End of the abilint section 96 97 implicit none 98 99 !Arguments ------------------------------------ 100 !scalars 101 integer,intent(in) :: my_natom,natom,nfft,nspden,nsym,ntypat,paral_kgb,prtefg,usepaw 102 integer,optional,intent(in) :: comm_atom 103 real(dp),intent(in) :: ucvol 104 type(MPI_type),intent(in) :: mpi_enreg 105 type(pawang_type),intent(in) :: pawang 106 !arrays 107 integer,intent(in) :: ngfft(18),symrel(3,3,nsym),typat(natom) 108 integer,optional,target,intent(in) :: mpi_atmtab(:) 109 real(dp),intent(in) :: ptcharge(ntypat) 110 real(dp),intent(in) :: quadmom(ntypat),rhor(nfft,nspden),rprimd(3,3) 111 real(dp),intent(in) :: tnons(3,nsym),zion(ntypat) 112 real(dp),intent(inout) :: xred(3,natom) 113 type(paw_an_type),intent(in) :: paw_an(my_natom) 114 type(pawrad_type),intent(in) :: pawrad(ntypat) 115 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 116 type(pawtab_type),intent(in) :: pawtab(ntypat) 117 118 !Local variables------------------------------- 119 !scalars 120 integer :: INFO,LDA,LWORK,N,iatom,my_comm_atom 121 logical :: my_atmtab_allocated,paral_atom 122 real(dp) :: cq,eta,vxx,vyy,vzz 123 character(len=500) :: message 124 !arrays 125 integer,pointer :: my_atmtab(:) 126 real(dp) :: eigval(3),matr(3,3),work(8) 127 real(dp),allocatable :: efg(:,:,:),efg_el(:,:,:),efg_ion(:,:,:),efg_paw(:,:,:) 128 real(dp),allocatable :: efg_point_charge(:,:,:) 129 130 ! ************************************************************************ 131 132 !Compatibility tests 133 if (usepaw /= 1) then 134 message = ' usepaw /= 1 but EFG calculation requires PAW ' 135 MSG_ERROR(message) 136 end if 137 138 !Set up parallelism over atoms 139 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 140 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 141 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 142 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 143 144 ABI_ALLOCATE(efg,(3,3,natom)) 145 ABI_ALLOCATE(efg_el,(3,3,natom)) 146 ABI_ALLOCATE(efg_ion,(3,3,natom)) 147 ABI_ALLOCATE(efg_paw,(3,3,natom)) 148 ABI_ALLOCATE(efg_point_charge,(3,3,natom)) 149 efg_el(:,:,:) = zero 150 efg_ion(:,:,:) = zero 151 efg_paw(:,:,:) = zero 152 efg_point_charge(:,:,:) = zero 153 154 call make_efg_el(efg_el,mpi_enreg,natom,nfft,ngfft,nspden,nsym,paral_kgb,rhor,rprimd,symrel,tnons,xred) 155 156 call make_efg_ion(efg_ion,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,zion) 157 158 if (paral_atom) then 159 call make_efg_onsite(efg_paw,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab,& 160 & rprimd,symrel,tnons,xred,comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 161 else 162 call make_efg_onsite(efg_paw,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab,& 163 & rprimd,symrel,tnons,xred) 164 end if 165 166 !calculate efg due to pure point charges, as input in variable ptcharge(ntypat) 167 !note here all atoms of the same type will have the same valence; in the future this 168 !could be made more flexible by having ptcharge(natom) but that will require a slightly 169 !different version than the existing make_efg_ion routine 170 if(prtefg > 2) then 171 call make_efg_ion(efg_point_charge,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,ptcharge) 172 end if 173 174 efg(:,:,:) = efg_el(:,:,:) + efg_ion(:,:,:) + efg_paw(:,:,:) 175 176 write(message,'(a,a,a)' ) ch10,' Electric Field Gradient Calculation ',ch10 177 call wrtout(ab_out,message,'COLL') 178 179 LDA=3; LWORK=8;N=3 ! these parameters are needed for the LAPACK dsyev routine 180 do iatom = 1, natom 181 matr(:,:) = efg(:,:,iatom) 182 call dsyev('V','U',N,matr,LDA,eigval,work,LWORK,INFO) ! get eigenvalues and eigenvectors 183 if (eigval(3) > abs(eigval(1)) ) then ! In NMR, the convention is that whatever component is 184 ! largest in magnitude is called Vzz, next comes Vxx, then Vyy 185 vzz = eigval(3) 186 vxx = eigval(1) 187 vyy = eigval(2) 188 else 189 vzz = eigval(1) 190 vxx = eigval(3) 191 vyy = eigval(2) 192 end if 193 if (abs(quadmom(typat(iatom))) > tol8 ) then ! only relevant when quadmom > 0 for a given atom 194 ! cq = (eQ)*Vzz/h, where Q is the electric quadrupole moment and Vzz is the largest in magnitude 195 ! principal component of the EFG tensor. Q is input in quadmom in barns, and Vzz is computed in atomic 196 ! units. The factor 2349647.81 Ha^{-1}Bohr^2 fm^{-2} sec^-1 converts from atomic units to frequency (see 197 ! http://www.ismar.org/ISMARpedia/index.php/Nuclear_Quadrupole_Resonance for discussion); we divide by 198 ! 10^6 to convert to MHz from Hz and multiply by 100 to convert from fm^2 to Barns. 199 cq = vzz*quadmom(typat(iatom))*2349647.81/1.0E4 200 if(abs(cq) > tol6 )then ! if Cq is non-zero, eta is meaningful, otherwise it s numerical noise 201 eta = abs(vxx - vyy)/abs(vzz) 202 else 203 eta=zero 204 end if 205 else 206 cq =zero 207 eta =zero 208 end if 209 ! we always write Cq and eta, these are the NMR observables 210 write(message,'(a,i3,a,i3,a,f13.6,a,f13.6)') ' Atom ',iatom,', typat ',typat(iatom),': Cq = ',cq,' MHz eta = ',eta 211 call wrtout(ab_out,message,'COLL') 212 if (prtefg > 1) then ! print detailed results on component EFG's 213 write(message,'(a,a,f13.6,a,a,3f13.6)')ch10,' efg eigval : ',eigval(1),ch10,& 214 & '- eigvec : ',matr(1,1),matr(2,1),matr(3,1) 215 call wrtout(ab_out,message,'COLL') 216 write(message,'(a,f13.6,a,a,3f13.6)')' efg eigval : ',eigval(2),ch10,& 217 & '- eigvec : ',matr(1,2),matr(2,2),matr(3,2) 218 call wrtout(ab_out,message,'COLL') 219 write(message,'(a,f13.6,a,a,3f13.6)')' efg eigval : ',eigval(3),ch10,& 220 & '- eigvec : ',matr(1,3),matr(2,3),matr(3,3) 221 call wrtout(ab_out,message,'COLL') 222 write(message,'(a,a,3f13.6)')ch10,' total efg : ',efg(1,1,iatom),efg(1,2,iatom),efg(1,3,iatom) 223 call wrtout(ab_out,message,'COLL') 224 write(message,'(a,3f13.6)')' total efg : ',efg(2,1,iatom),efg(2,2,iatom),efg(2,3,iatom) 225 call wrtout(ab_out,message,'COLL') 226 write(message,'(a,3f13.6,a)')' total efg : ',efg(3,1,iatom),efg(3,2,iatom),efg(3,3,iatom),ch10 227 call wrtout(ab_out,message,'COLL') 228 write(message,'(a,a,3f13.6)')ch10,' efg_el : ',efg_el(1,1,iatom),efg_el(1,2,iatom),efg_el(1,3,iatom) 229 call wrtout(ab_out,message,'COLL') 230 write(message,'(a,3f13.6)')' efg_el : ',efg_el(2,1,iatom),efg_el(2,2,iatom),efg_el(2,3,iatom) 231 call wrtout(ab_out,message,'COLL') 232 write(message,'(a,3f13.6,a)')' efg_el : ',efg_el(3,1,iatom),efg_el(3,2,iatom),efg_el(3,3,iatom),ch10 233 call wrtout(ab_out,message,'COLL') 234 write(message,'(a,3f13.6)')' efg_ion : ',efg_ion(1,1,iatom),efg_ion(1,2,iatom),efg_ion(1,3,iatom) 235 call wrtout(ab_out,message,'COLL') 236 write(message,'(a,3f13.6)')' efg_ion : ',efg_ion(2,1,iatom),efg_ion(2,2,iatom),efg_ion(2,3,iatom) 237 call wrtout(ab_out,message,'COLL') 238 write(message,'(a,3f13.6,a)')' efg_ion : ',efg_ion(3,1,iatom),efg_ion(3,2,iatom),efg_ion(3,3,iatom),ch10 239 call wrtout(ab_out,message,'COLL') 240 write(message,'(a,3f13.6)')' efg_paw : ',efg_paw(1,1,iatom),efg_paw(1,2,iatom),efg_paw(1,3,iatom) 241 call wrtout(ab_out,message,'COLL') 242 write(message,'(a,3f13.6)')' efg_paw : ',efg_paw(2,1,iatom),efg_paw(2,2,iatom),efg_paw(2,3,iatom) 243 call wrtout(ab_out,message,'COLL') 244 write(message,'(a,3f13.6,a)')' efg_paw : ',efg_paw(3,1,iatom),efg_paw(3,2,iatom),efg_paw(3,3,iatom),ch10 245 call wrtout(ab_out,message,'COLL') 246 end if 247 if (prtefg > 2) then ! write output of pure pointcharge calculation 248 matr(:,:) = efg_point_charge(:,:,iatom) 249 call dsyev('V','U',N,matr,LDA,eigval,work,LWORK,INFO) ! get eigenvalues and eigenvectors 250 if (eigval(3) > abs(eigval(1)) ) then ! In NMR, the convention is that whatever component is 251 ! largest in magnitude is called Vzz, next comes Vxx, then Vyy 252 vzz = eigval(3) 253 vxx = eigval(1) 254 vyy = eigval(2) 255 else 256 vzz = eigval(1) 257 vxx = eigval(3) 258 vyy = eigval(2) 259 end if 260 if (abs(quadmom(typat(iatom))) > tol8 ) then ! only relevant when quadmom > 0 for a given atom 261 ! cq = e2Qq/h, where Vzz = eq and quadmom = Q; the other factors convert from atomic units to MHz 262 cq = vzz*quadmom(typat(iatom))*2349647.81/1.0E4 263 if(abs(cq) > tol6 )then ! if Cq is non-zero, eta is meaningful, otherwise it s numerical noise 264 eta = abs(vxx - vyy)/abs(vzz) 265 else 266 eta=zero 267 end if 268 else 269 cq =zero 270 eta =zero 271 end if 272 ! we always write Cq and eta, these are the NMR observables 273 write(message,'(a,i3,a,i3,a,f13.6,a,f13.6)') ' Atom ',iatom,', typat ',typat(iatom),& 274 & ': Point charge Cq = ',cq,' MHz eta = ',eta 275 call wrtout(ab_out,message,'COLL') 276 write(message,'(a,a,f13.6,a,a,3f13.6)')ch10,' point charge efg eigval : ',eigval(1),ch10,& 277 & '- eigvec : ',matr(1,1),matr(2,1),matr(3,1) 278 call wrtout(ab_out,message,'COLL') 279 write(message,'(a,f13.6,a,a,3f13.6)')' point charge efg eigval : ',eigval(2),ch10,& 280 & '- eigvec : ',matr(1,2),matr(2,2),matr(3,2) 281 call wrtout(ab_out,message,'COLL') 282 write(message,'(a,f13.6,a,a,3f13.6)')' point charge efg eigval : ',eigval(3),ch10,& 283 & '- eigvec : ',matr(1,3),matr(2,3),matr(3,3) 284 call wrtout(ab_out,message,'COLL') 285 write(message,'(a,a,3f13.6)')ch10,' point charge efg : ',efg_point_charge(1,1,iatom),& 286 & efg_point_charge(1,2,iatom),efg_point_charge(1,3,iatom) 287 call wrtout(ab_out,message,'COLL') 288 write(message,'(a,3f13.6)')' point charge efg : ',efg_point_charge(2,1,iatom),& 289 & efg_point_charge(2,2,iatom),efg_point_charge(2,3,iatom) 290 call wrtout(ab_out,message,'COLL') 291 write(message,'(a,3f13.6,a)')' point charge efg : ',efg_point_charge(3,1,iatom),& 292 & efg_point_charge(3,2,iatom),efg_point_charge(3,3,iatom),ch10 293 call wrtout(ab_out,message,'COLL') 294 end if 295 end do 296 write(message,'(3a)')ch10,ch10,ch10 297 call wrtout(ab_out,message,'COLL') 298 299 ABI_DEALLOCATE(efg) 300 ABI_DEALLOCATE(efg_el) 301 ABI_DEALLOCATE(efg_ion) 302 ABI_DEALLOCATE(efg_paw) 303 ABI_DEALLOCATE(efg_point_charge) 304 305 !Destroy atom table used for parallelism 306 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 307 308 !DEBUG 309 !write(std_out,*)' calc_efg : exit ' 310 !stop 311 !ENDDEBUG 312 313 end subroutine calc_efg