TABLE OF CONTENTS


ABINIT/m_paw_nmr [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_nmr

FUNCTION

  This module contains routines related to Nuclear Magnetic Resonance (NMR)
   observables (PAW approach).

COPYRIGHT

 Copyright (C) 2018-2018 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

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 MODULE m_paw_nmr
25 
26 
27  use defs_basis
28  use m_errors
29  use m_abicore
30  use m_xmpi
31 
32  use m_symtk,      only : matpointsym
33  use m_pawang,     only : pawang_type
34  use m_pawtab,     only : pawtab_type
35  use m_pawrad,     only : pawrad_type,pawrad_deducer0,simp_gen
36  use m_pawtab,     only : pawtab_type
37  use m_paw_an,     only : paw_an_type
38  use m_pawrhoij,   only : pawrhoij_type
39  use m_paw_denpot, only : pawdensities
40  use m_paral_atom, only : get_my_atmtab,free_my_atmtab
41 
42  implicit none
43 
44  private
45 
46 !public procedures.
47  public :: make_efg_onsite ! Compute the electric field gradient due to PAW on-site densities
48  public :: make_fc_paw     ! Compute the PAW on-site contribution to the Fermi-contact
49 
50 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.

PARENTS

      calc_efg

CHILDREN

      free_my_atmtab,get_my_atmtab,matpointsym,pawdensities,pawrad_deducer0
      simp_gen,xmpi_sum

SOURCE

101 subroutine make_efg_onsite(efg,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab, &
102 &                          rprimd,symrel,tnons,xred,&
103 &                          mpi_atmtab,comm_atom) ! optional arguments (parallelism)
104 
105 
106 !This section has been created automatically by the script Abilint (TD).
107 !Do not modify the following lines by hand.
108 #undef ABI_FUNC
109 #define ABI_FUNC 'make_efg_onsite'
110 !End of the abilint section
111 
112  implicit none
113 
114 !Arguments ------------------------------------
115 !scalars
116  integer,intent(in) :: my_natom,natom,nsym,ntypat
117  integer,optional,intent(in) :: comm_atom
118  type(pawang_type),intent(in) :: pawang
119 !arrays
120  integer,intent(in) :: symrel(3,3,nsym)
121  integer,optional,target,intent(in) :: mpi_atmtab(:)
122  real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym),xred(3,natom)
123  real(dp),intent(out) :: efg(3,3,natom)
124  type(paw_an_type),intent(in) :: paw_an(my_natom)
125  type(pawrad_type),intent(in) :: pawrad(ntypat)
126  type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
127  type(pawtab_type),intent(in) :: pawtab(ntypat)
128 
129 !Local variables-------------------------------
130 !scalars
131  integer :: cplex,iatom,iatom_tot,ictr,ierr,imesh_size,ispden,itypat
132  integer :: lm,lm_size,lnspden,mesh_size,my_comm_atom,nzlmopt,nspden
133  integer :: opt_compch,opt_dens,opt_l,opt_print
134  logical :: my_atmtab_allocated,paral_atom
135  real(dp) :: c1,c2,c3,compch_sph,intg
136 !arrays
137  integer,pointer :: my_atmtab(:)
138  logical,allocatable :: lmselectin(:),lmselectout(:)
139  real(dp),allocatable :: ff(:),nhat1(:,:,:),rho1(:,:,:),trho1(:,:,:)
140 
141 ! ************************************************************************
142 
143  DBG_ENTER("COLL")
144 
145  efg(:,:,:) = zero
146 
147 !Set up parallelism over atoms
148  paral_atom=(present(comm_atom).and.(my_natom/=natom))
149  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
150  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
151  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
152 & my_natom_ref=my_natom)
153 
154 !the following factors arise in expanding the angular dependence of the electric field gradient tensor in
155 !terms of real spherical harmonics. The real spherical harmonics are as in the routine initylmr.F90; see
156 !in particular also http://www.unioviedo.es/qcg/art/Theochem419-19-ov-BF97-rotation-matrices.pdf
157  c1 = sqrt(16.0*pi/5.0)
158  c2 = sqrt(4.0*pi/5.0)
159  c3 = sqrt(12.0*pi/5.0)
160 
161 !loop over atoms in cell
162  do iatom = 1, my_natom
163    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
164    itypat=pawrhoij(iatom)%itypat
165 
166    lm_size = paw_an(iatom)%lm_size
167    if (lm_size < 5) cycle ! if lm_size < 5 then the on-site densities for this atom have no L=2 component
168 !  and therefore nothing to contribute to the on-site electric field gradient
169 
170    mesh_size=pawtab(itypat)%mesh_size
171    ABI_ALLOCATE(ff,(mesh_size))
172 
173    cplex = pawrhoij(iatom)%cplex
174    nspden = pawrhoij(iatom)%nspden
175    ABI_ALLOCATE(lmselectin,(lm_size))
176    ABI_ALLOCATE(lmselectout,(lm_size))
177    lmselectin = .true. ! compute all moments of densities
178    nzlmopt = -1
179    opt_compch = 0
180    compch_sph = zero
181    opt_dens = 0 ! compute all densities
182    opt_l = -1 ! all moments contribute
183    opt_print = 0 ! do not print out moments
184 
185    ABI_ALLOCATE(nhat1,(cplex*mesh_size,lm_size,nspden))
186    ABI_ALLOCATE(rho1,(cplex*mesh_size,lm_size,nspden))
187    ABI_ALLOCATE(trho1,(cplex*mesh_size,lm_size,nspden))
188 
189 !  loop over spin components
190 !  nspden = 1: just a single component
191 !  nspden = 2: two components, loop over them
192 !  nspden = 4: total is in component 1, only one of interest
193    if ( nspden == 2 ) then
194      lnspden = 2
195    else
196      lnspden = 1
197    end if
198    do ispden=1,lnspden
199 
200 !    construct multipole expansion of on-site charge densities for this atom
201      call pawdensities(compch_sph,cplex,iatom_tot,lmselectin,lmselectout,lm_size,&
202 &     nhat1,nspden,nzlmopt,opt_compch,opt_dens,opt_l,opt_print,&
203 &     pawang,0,pawrad(itypat),pawrhoij(iatom),pawtab(itypat),&
204 &     rho1,trho1)
205 
206      do lm = 5, 9 ! loop on L=2 components of multipole expansion
207 
208        if(.not. lmselectout(lm)) cycle ! skip moments that contributes zero
209 
210 !      the following is r^2*(n1-tn1-nhat)/r^3 for this multipole moment
211 !      use only the real part of the density in case of cplex==2
212        do imesh_size = 2, mesh_size
213          ictr = cplex*(imesh_size - 1) + 1
214          ff(imesh_size)=(rho1(ictr,lm,ispden)-trho1(ictr,lm,ispden)-&
215 &         nhat1(ictr,lm,ispden))/&
216 &         pawrad(itypat)%rad(imesh_size)
217        end do
218        call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
219        call simp_gen(intg,ff,pawrad(itypat))
220        select case (lm)
221        case (5) ! S_{2,-2}
222          efg(1,2,iatom_tot) = efg(1,2,iatom_tot) - c3*intg ! xy case
223        case (6) ! S_{2,-1}
224          efg(2,3,iatom_tot) = efg(2,3,iatom_tot) - c3*intg ! yz case
225        case (7) ! S_{2, 0}
226          efg(1,1,iatom_tot) = efg(1,1,iatom_tot) + c2*intg ! xx case
227          efg(2,2,iatom_tot) = efg(2,2,iatom_tot) + c2*intg ! yy case
228          efg(3,3,iatom_tot) = efg(3,3,iatom_tot) - c1*intg ! zz case
229        case (8) ! S_{2,+1}
230          efg(1,3,iatom_tot) = efg(1,3,iatom_tot) - c3*intg ! xz case
231        case (9) ! S_{2,+2}
232          efg(1,1,iatom_tot) = efg(1,1,iatom_tot) - c3*intg ! xx case
233          efg(2,2,iatom_tot) = efg(2,2,iatom_tot) + c3*intg ! yy case
234        end select
235 
236      end do  ! end loop over LM components with L=2
237 
238    end do    ! Loop on spin components
239 
240 !  Symmetrization of EFG
241    efg(2,1,iatom_tot) = efg(1,2,iatom_tot)
242    efg(3,1,iatom_tot) = efg(1,3,iatom_tot)
243    efg(3,2,iatom_tot) = efg(2,3,iatom_tot)
244 
245    ABI_DEALLOCATE(lmselectin)
246    ABI_DEALLOCATE(lmselectout)
247    ABI_DEALLOCATE(ff)
248    ABI_DEALLOCATE(nhat1)
249    ABI_DEALLOCATE(rho1)
250    ABI_DEALLOCATE(trho1)
251 
252  end do     ! Loop on atoms
253 
254 !Reduction in case of parallelisation over atoms
255  if (paral_atom) then
256    call xmpi_sum(efg,my_comm_atom,ierr)
257  end if
258 
259 ! symmetrize tensor at each atomic site using point symmetry operations
260  do iatom = 1, natom
261    call matpointsym(iatom,efg(:,:,iatom),natom,nsym,rprimd,symrel,tnons,xred)
262  end do
263 
264 !Destroy atom table used for parallelism
265  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
266 
267  DBG_EXIT("COLL")
268 
269  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]].

PARENTS

      calc_fc

CHILDREN

      free_my_atmtab,get_my_atmtab,xmpi_sum

SOURCE

318 subroutine make_fc_paw(fc,my_natom,natom,nspden,ntypat,pawrhoij,pawrad,pawtab,&
319 &                      mpi_atmtab,comm_atom) ! optional arguments (parallelism)
320 
321 
322 !This section has been created automatically by the script Abilint (TD).
323 !Do not modify the following lines by hand.
324 #undef ABI_FUNC
325 #define ABI_FUNC 'make_fc_paw'
326 !End of the abilint section
327 
328  implicit none
329 
330 !Arguments ------------------------------------
331 !scalars
332  integer,intent(in) :: my_natom,natom,nspden,ntypat
333  integer,optional,intent(in) :: comm_atom
334 !arrays
335  integer,optional,target,intent(in) :: mpi_atmtab(:)
336  real(dp),intent(out) :: fc(nspden,natom)
337  type(pawrad_type),intent(in) :: pawrad(ntypat)
338  type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
339  type(pawtab_type),target,intent(in) :: pawtab(ntypat)
340 
341 !Local variables-------------------------------
342 !scalars
343  integer :: iatom,iatom_tot,ierr,irhoij,islope,ispden,itypat
344  integer :: ilmn,il,iln,ilm,im,jl,jlm,jlmn,jln,jm,j0lmn
345  integer :: klmn,kln,mesh_size,my_comm_atom,nslope
346  logical :: my_atmtab_allocated,paral_atom
347  real(dp) :: mi,mj,xi,xxsum,xysumi,xysumj,yi,yj
348 !arrays
349  integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
350  integer,pointer :: my_atmtab(:)
351 
352 ! ************************************************************************
353 
354  DBG_ENTER("COLL")
355 
356 !Set up parallelism over atoms
357  paral_atom=(present(comm_atom).and.(my_natom/=natom))
358  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
359  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
360  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
361 
362 !Initialization
363  fc(:,:)=zero
364 
365 !number of points to use in computing initial slopes of radial functions
366  nslope = 5
367 
368 !loop over atoms in cell
369  do iatom = 1, my_natom
370    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
371    itypat = pawrhoij(iatom)%itypat
372    mesh_size=pawtab(itypat)%mesh_size
373    indlmn => pawtab(itypat)%indlmn
374 
375 !  loop over spin components
376    do ispden=1,nspden
377 
378 !    loop over basis elements for this atom
379 !    ----
380      do jlmn=1,pawtab(itypat)%lmn_size
381        jl=indlmn(1,jlmn)
382        jm=indlmn(2,jlmn)
383        jlm=indlmn(4,jlmn)
384        jln=indlmn(5,jlmn)
385        j0lmn=jlmn*(jlmn-1)/2
386        do ilmn=1,jlmn
387          il=indlmn(1,ilmn)
388          im=indlmn(2,ilmn)
389          iln=indlmn(5,ilmn)
390          ilm=indlmn(4,ilmn)
391          klmn=j0lmn+ilmn
392          kln = pawtab(itypat)%indklmn(2,klmn) ! need this for mesh selection below
393 
394          if (jl==0 .and. il==0) then ! select only s-states
395 
396 !          Loop over non-zero elements of rhoij
397            do irhoij=1,pawrhoij(iatom)%nrhoijsel
398              if (klmn==pawrhoij(iatom)%rhoijselect(irhoij)) then ! rho_ij /= 0 for this klmn
399                xxsum = 0 ! these three variables will be used to compute the slopes
400                xysumi = 0
401                xysumj = 0
402                do islope=1, nslope
403                  xi=0
404                  if(pawrad(itypat)%mesh_type == 1) xi = (islope - 1)*pawrad(itypat)%rstep
405                  if(pawrad(itypat)%mesh_type == 2) xi = pawrad(itypat)%rstep * &
406 &                 (exp(pawrad(itypat)%lstep * (islope - 1)) - 1)
407                  if(pawrad(itypat)%mesh_type == 3) then
408                    if (islope == 1) then
409                      xi = 0
410                    else
411                      xi = pawrad(itypat)%rstep * exp(pawrad(itypat)%lstep*(islope-1))
412                    end if
413                  end if
414                  if(pawrad(itypat)%mesh_type == 4) xi = &
415 &                 -pawrad(itypat)%rstep*log(1.0-(islope-1)/pawrad(itypat)%mesh_size)
416                  yi = pawtab(itypat)%phi(islope,iln) ! function value for u_i
417                  yj = pawtab(itypat)%phi(islope,jln) ! function value for u_j
418                  xxsum =  xxsum + xi*xi
419                  xysumi = xysumi + xi*yi
420                  xysumj = xysumj + xi*yj
421                end do
422 !              the slopes of the radial functions are obtained by minimizing
423 !              chi = sum(y_i - m*x_i)^2 (in other words, a linear least squares
424 !              fit constrained to go through the origin)
425 !              the result is m = sum(y_i*x_i)/sum(x_i*x_i)
426                mi = xysumi/xxsum
427                mj = xysumj/xxsum
428 !              accumulate the rho_ij contribution to the fermi contact for this spin density:
429                if (pawrhoij(iatom)%cplex == 1) then
430                  fc(ispden,iatom_tot)=fc(ispden,iatom_tot)+&
431 &                 pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*mi*mj/four_pi
432                else
433                  fc(ispden,iatom_tot)=fc(ispden,iatom_tot)+&
434 &                 pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*mi*mj/four_pi
435                end if
436              end if ! end selection on klmn for nonzero rho_ij
437            end do ! end loop over nonzero rho_ij
438          end if ! end l=l'=0 selection
439        end do ! end loop over ilmn
440      end do ! end loop over jlmn
441    end do ! end loop over spin densities
442  end do     ! Loop on atoms
443 
444 !Reduction in case of parallelisation over atoms
445  if (paral_atom) then
446    call xmpi_sum(fc,my_comm_atom,ierr)
447  end if
448 
449 !Destroy atom table used for parallelism
450  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
451 
452  DBG_EXIT("COLL")
453 
454  end subroutine make_fc_paw