TABLE OF CONTENTS
ABINIT/m_paw_nmr [ Modules ]
NAME
m_paw_nmr
FUNCTION
This module contains routines related to Nuclear Magnetic Resonance (NMR) observables (PAW approach).
COPYRIGHT
Copyright (C) 2018-2024 ABINIT group (JWZ, 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 .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_paw_nmr 24 25 26 use defs_basis 27 use m_errors 28 use m_abicore 29 use m_xmpi 30 31 use m_symtk, only : matpointsym 32 use m_pawang, only : pawang_type 33 use m_pawtab, only : pawtab_type 34 use m_pawrad, only : pawrad_type,pawrad_deducer0,simp_gen 35 use m_pawtab, only : pawtab_type 36 use m_paw_an, only : paw_an_type 37 use m_pawrhoij, only : pawrhoij_type 38 use m_paw_denpot, only : pawdensities 39 use m_paral_atom, only : get_my_atmtab,free_my_atmtab 40 41 implicit none 42 43 private 44 45 !public procedures. 46 public :: make_efg_onsite ! Compute the electric field gradient due to PAW on-site densities 47 public :: make_fc_paw ! Compute the PAW on-site contribution to the Fermi-contact 48 49 CONTAINS !========================================================================================
m_paw_nmr/make_efg_onsite [ Functions ]
[ Top ] [ m_paw_nmr ] [ Functions ]
NAME
make_efg_onsite
FUNCTION
Compute the electric field gradient due to onsite densities
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. nsym=number of symmetries in space group ntypat=number of atom types 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 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 site due to nhat
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]]. See in particular
Eq. 11 and 12 of Profeta et al., but note that their sum over occupied states times 2 for occupation number is
replaced in the Kresse and Joubert formulation by the sum over $\rho_{ij}$ occupations for each basis element pair.
SOURCE
93 subroutine make_efg_onsite(efg,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab, & 94 & rprimd,symrel,tnons,xred,& 95 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 96 97 implicit none 98 99 !Arguments ------------------------------------ 100 !scalars 101 integer,intent(in) :: my_natom,natom,nsym,ntypat 102 integer,optional,intent(in) :: comm_atom 103 type(pawang_type),intent(in) :: pawang 104 !arrays 105 integer,intent(in) :: symrel(3,3,nsym) 106 integer,optional,target,intent(in) :: mpi_atmtab(:) 107 real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym),xred(3,natom) 108 real(dp),intent(out) :: efg(3,3,natom) 109 type(paw_an_type),intent(in) :: paw_an(my_natom) 110 type(pawrad_type),intent(in) :: pawrad(ntypat) 111 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 112 type(pawtab_type),intent(in) :: pawtab(ntypat) 113 114 !Local variables------------------------------- 115 !scalars 116 integer :: cplex,iatom,iatom_tot,ictr,ierr,imesh_size,ispden,itypat 117 integer :: lm,lm_size,local_paw_print_vol 118 integer :: mesh_size,my_comm_atom,nzlmopt,nspden 119 integer :: opt_compch,opt_dens,opt_l,opt_print 120 logical :: my_atmtab_allocated,paral_atom 121 real(dp) :: c1,c2,c3,compch_sph,intg 122 !arrays 123 integer,pointer :: my_atmtab(:) 124 logical,allocatable :: lmselectin(:),lmselectout(:) 125 real(dp),allocatable :: ff(:),nhat1(:,:,:),rho1(:,:,:),trho1(:,:,:) 126 127 ! ************************************************************************ 128 129 DBG_ENTER("COLL") 130 131 if (my_natom>0) then 132 ABI_CHECK(pawrhoij(1)%qphase==1,'make_efg_onsite: not supposed to be called with qphqse=2!') 133 end if 134 135 efg(:,:,:) = zero 136 137 !Set up parallelism over atoms 138 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 139 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 140 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 141 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,& 142 & my_natom_ref=my_natom) 143 144 !the following factors arise in expanding the angular dependence of the electric field gradient tensor in 145 !terms of real spherical harmonics. The real spherical harmonics are as in the routine initylmr.F90; see 146 !in particular also http://www.unioviedo.es/qcg/art/Theochem419-19-ov-BF97-rotation-matrices.pdf 147 c1 = sqrt(16.0*pi/5.0) 148 c2 = sqrt(4.0*pi/5.0) 149 c3 = sqrt(12.0*pi/5.0) 150 151 !loop over atoms in cell 152 do iatom = 1, my_natom 153 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 154 itypat=pawrhoij(iatom)%itypat 155 156 lm_size = paw_an(iatom)%lm_size 157 if (lm_size < 5) cycle ! if lm_size < 5 then the on-site densities for this atom have no L=2 component 158 ! and therefore nothing to contribute to the on-site electric field gradient 159 160 mesh_size=pawtab(itypat)%mesh_size 161 ABI_MALLOC(ff,(mesh_size)) 162 163 cplex = pawrhoij(iatom)%qphase 164 nspden = pawrhoij(iatom)%nspden 165 ABI_MALLOC(lmselectin,(lm_size)) 166 ABI_MALLOC(lmselectout,(lm_size)) 167 lmselectin = .true. ! compute all moments of densities 168 nzlmopt = -1 169 opt_compch = 0 170 compch_sph = zero 171 opt_dens = 0 ! compute all densities 172 opt_l = -1 ! all moments contribute 173 opt_print = 0 ! do not print out moments 174 local_paw_print_vol = 0 ! standard amount of printing in pawdensities 175 176 ABI_MALLOC(nhat1,(cplex*mesh_size,lm_size,nspden*(1-((opt_dens+1)/2)))) 177 ABI_MALLOC(rho1,(cplex*mesh_size,lm_size,nspden)) 178 ABI_MALLOC(trho1,(cplex*mesh_size,lm_size,nspden*(1-((opt_dens+1)/2)))) 179 180 ! construct multipole expansion of on-site charge densities for this atom 181 call pawdensities(compch_sph,cplex,iatom_tot,lmselectin,lmselectout,lm_size,& 182 & nhat1,nspden,nzlmopt,opt_compch,opt_dens,opt_l,opt_print,& 183 & pawang,local_paw_print_vol,pawrad(itypat),pawrhoij(iatom),pawtab(itypat),& 184 & rho1,trho1) 185 186 ! spin components: 187 ! nspden(1) contains total in all cases 188 ispden = 1 189 190 do lm = 5, 9 ! loop on L=2 components of multipole expansion 191 192 if(.not. lmselectout(lm)) cycle ! skip moments that contributes zero 193 194 ! the following is r^2*(n1-tn1-nhat)/r^3 for this multipole moment 195 ! use only the real part of the density in case of cplex==2 196 do imesh_size = 2, mesh_size 197 ictr = cplex*(imesh_size - 1) + 1 198 ff(imesh_size)=rho1(ictr,lm,ispden)-trho1(ictr,lm,ispden)-nhat1(ictr,lm,ispden) 199 ff(imesh_size)=ff(imesh_size)/pawrad(itypat)%rad(imesh_size) 200 end do 201 call pawrad_deducer0(ff,mesh_size,pawrad(itypat)) 202 call simp_gen(intg,ff,pawrad(itypat)) 203 select case (lm) 204 case (5) ! S_{2,-2} 205 efg(1,2,iatom_tot) = efg(1,2,iatom_tot) - c3*intg ! xy case 206 case (6) ! S_{2,-1} 207 efg(2,3,iatom_tot) = efg(2,3,iatom_tot) - c3*intg ! yz case 208 case (7) ! S_{2, 0} 209 efg(1,1,iatom_tot) = efg(1,1,iatom_tot) + c2*intg ! xx case 210 efg(2,2,iatom_tot) = efg(2,2,iatom_tot) + c2*intg ! yy case 211 efg(3,3,iatom_tot) = efg(3,3,iatom_tot) - c1*intg ! zz case 212 case (8) ! S_{2,+1} 213 efg(1,3,iatom_tot) = efg(1,3,iatom_tot) - c3*intg ! xz case 214 case (9) ! S_{2,+2} 215 efg(1,1,iatom_tot) = efg(1,1,iatom_tot) - c3*intg ! xx case 216 efg(2,2,iatom_tot) = efg(2,2,iatom_tot) + c3*intg ! yy case 217 end select 218 219 end do ! end loop over LM components with L=2 220 221 222 ! Symmetrization of EFG 223 efg(2,1,iatom_tot) = efg(1,2,iatom_tot) 224 efg(3,1,iatom_tot) = efg(1,3,iatom_tot) 225 efg(3,2,iatom_tot) = efg(2,3,iatom_tot) 226 227 ABI_FREE(lmselectin) 228 ABI_FREE(lmselectout) 229 ABI_FREE(ff) 230 ABI_FREE(nhat1) 231 ABI_FREE(rho1) 232 ABI_FREE(trho1) 233 234 end do ! Loop on atoms 235 236 !Reduction in case of parallelisation over atoms 237 if (paral_atom) then 238 call xmpi_sum(efg,my_comm_atom,ierr) 239 end if 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 !Destroy atom table used for parallelism 247 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 248 249 DBG_EXIT("COLL") 250 251 end subroutine make_efg_onsite
m_paw_nmr/make_fc_paw [ Functions ]
[ Top ] [ m_paw_nmr ] [ Functions ]
NAME
make_fc_paw
FUNCTION
Compute the Fermi-contact term due to the PAW cores
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 ensity component ntypat=number of atom types 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
OUTPUT
fc(nspden,natom)=the Fermi-contact interaction at each site due to PAW for each spin density
NOTES
The Fermi contact interaction is the electron density evaluated exactly at the nuclear site.
For a nuclear site at R, we are thus computing the expectation value of $\delta^3(R)$, the
the three-dimensional delta function at vector position $R$. In terms of the radial variable only
the delta function is $\delta(r)/4\pi r^2$. Because this observable is
absolutely confined within the PAW radius, only the response due to the AE PAW functions is
needed, the pseudo wavefunctions and pseudo PAW functions cancel each other out. We then
must compute the integral of $u_i/r times u_j/r \delta(R)d^3r$, for the $l=0$ angular momentum
states only. This is simplified with the use of L'H\^{o}spital's theorem to take the limit
as $r\rightarrow 0$, yielding $u_i'(r) u_j'(r)$. To compute the derivatives we just fit the
first 5 points of the $u$ functions to a line through the origin, using the least squares
procedure resulting from $\chi = sum_i (y_i - m*x_i)^2$ . This is more stable than
computing the derivative of the whole function and extrapolating it to zero.
See Zwanziger, J. Phys. Conden. Matt. 21, 15024-15036 (2009) [[cite:Zwanziger2009]].
SOURCE
294 subroutine make_fc_paw(fc,my_natom,natom,nspden,ntypat,pawrhoij,pawrad,pawtab,& 295 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 296 297 implicit none 298 299 !Arguments ------------------------------------ 300 !scalars 301 integer,intent(in) :: my_natom,natom,nspden,ntypat 302 integer,optional,intent(in) :: comm_atom 303 !arrays 304 integer,optional,target,intent(in) :: mpi_atmtab(:) 305 real(dp),intent(out) :: fc(nspden,natom) 306 type(pawrad_type),intent(in) :: pawrad(ntypat) 307 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 308 type(pawtab_type),target,intent(in) :: pawtab(ntypat) 309 310 !Local variables------------------------------- 311 !scalars 312 integer :: iatom,iatom_tot,ierr,irhoij,islope,ispden,itypat 313 integer :: ilmn,il,iln,ilm,im,jl,jlm,jlmn,jln,jm,j0lmn,jrhoij 314 integer :: klmn,kln,mesh_size,my_comm_atom,nslope 315 logical :: my_atmtab_allocated,paral_atom 316 real(dp) :: mi,mj,xi,xxsum,xysumi,xysumj,yi,yj 317 !arrays 318 integer,ABI_CONTIGUOUS pointer :: indlmn(:,:) 319 integer,pointer :: my_atmtab(:) 320 321 ! ************************************************************************ 322 323 DBG_ENTER("COLL") 324 325 if (my_natom>0) then 326 ABI_CHECK(pawrhoij(1)%qphase==1,'make_fc_paw: not supposed to be called with qphqse=2!') 327 end if 328 329 !Set up parallelism over atoms 330 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 331 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 332 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 333 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 334 335 !Initialization 336 fc(:,:)=zero 337 338 !number of points to use in computing initial slopes of radial functions 339 nslope = 5 340 341 !loop over atoms in cell 342 do iatom = 1, my_natom 343 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 344 itypat = pawrhoij(iatom)%itypat 345 mesh_size=pawtab(itypat)%mesh_size 346 indlmn => pawtab(itypat)%indlmn 347 348 ! loop over spin components 349 do ispden=1,nspden 350 351 ! loop over basis elements for this atom 352 ! ---- 353 do jlmn=1,pawtab(itypat)%lmn_size 354 jl=indlmn(1,jlmn) 355 jm=indlmn(2,jlmn) 356 jlm=indlmn(4,jlmn) 357 jln=indlmn(5,jlmn) 358 j0lmn=jlmn*(jlmn-1)/2 359 do ilmn=1,jlmn 360 il=indlmn(1,ilmn) 361 im=indlmn(2,ilmn) 362 iln=indlmn(5,ilmn) 363 ilm=indlmn(4,ilmn) 364 klmn=j0lmn+ilmn 365 kln = pawtab(itypat)%indklmn(2,klmn) ! need this for mesh selection below 366 367 if (jl==0 .and. il==0) then ! select only s-states 368 369 ! Loop over non-zero elements of rhoij 370 jrhoij=1 371 do irhoij=1,pawrhoij(iatom)%nrhoijsel 372 if (klmn==pawrhoij(iatom)%rhoijselect(irhoij)) then ! rho_ij /= 0 for this klmn 373 xxsum = 0 ! these three variables will be used to compute the slopes 374 xysumi = 0 375 xysumj = 0 376 do islope=1, nslope 377 xi=0 378 if(pawrad(itypat)%mesh_type == 1) xi = (islope - 1)*pawrad(itypat)%rstep 379 if(pawrad(itypat)%mesh_type == 2) xi = pawrad(itypat)%rstep * & 380 & (exp(pawrad(itypat)%lstep * (islope - 1)) - 1) 381 if(pawrad(itypat)%mesh_type == 3) then 382 if (islope == 1) then 383 xi = 0 384 else 385 xi = pawrad(itypat)%rstep * exp(pawrad(itypat)%lstep*(islope-1)) 386 end if 387 end if 388 if(pawrad(itypat)%mesh_type == 4) xi = & 389 & -pawrad(itypat)%rstep*log(1.0-(islope-1)/pawrad(itypat)%mesh_size) 390 yi = pawtab(itypat)%phi(islope,iln) ! function value for u_i 391 yj = pawtab(itypat)%phi(islope,jln) ! function value for u_j 392 xxsum = xxsum + xi*xi 393 xysumi = xysumi + xi*yi 394 xysumj = xysumj + xi*yj 395 end do 396 ! the slopes of the radial functions are obtained by minimizing 397 ! chi = sum(y_i - m*x_i)^2 (in other words, a linear least squares 398 ! fit constrained to go through the origin) 399 ! the result is m = sum(y_i*x_i)/sum(x_i*x_i) 400 mi = xysumi/xxsum 401 mj = xysumj/xxsum 402 ! accumulate the rho_ij contribution to the fermi contact for this spin density: 403 fc(ispden,iatom_tot)=fc(ispden,iatom_tot)+& 404 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(jrhoij,ispden)*mi*mj/four_pi 405 end if ! end selection on klmn for nonzero rho_ij 406 jrhoij=jrhoij+pawrhoij(iatom)%cplex_rhoij 407 end do ! end loop over nonzero rho_ij 408 end if ! end l=l'=0 selection 409 end do ! end loop over ilmn 410 end do ! end loop over jlmn 411 end do ! end loop over spin densities 412 end do ! Loop on atoms 413 414 !Reduction in case of parallelisation over atoms 415 if (paral_atom) then 416 call xmpi_sum(fc,my_comm_atom,ierr) 417 end if 418 419 !Destroy atom table used for parallelism 420 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 421 422 DBG_EXIT("COLL") 423 424 end subroutine make_fc_paw