Line data Source code
1 : !!****m* ABINIT/m_paw_nmr
2 : !! NAME
3 : !! m_paw_nmr
4 : !!
5 : !! FUNCTION
6 : !! This module contains routines related to Nuclear Magnetic Resonance (NMR)
7 : !! observables (PAW approach).
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2018-2026 ABINIT group (JWZ, MT)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
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 !========================================================================================
50 : !!***
51 :
52 : !--------------------------------------------------------------------------------------------------
53 :
54 : !!****f* m_paw_nmr/make_efg_onsite
55 : !! NAME
56 : !! make_efg_onsite
57 : !!
58 : !! FUNCTION
59 : !! Compute the electric field gradient due to onsite densities
60 : !!
61 : !! INPUTS
62 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
63 : !! comm_atom=--optional-- MPI communicator over atoms
64 : !! my_natom=number of atoms treated by current processor
65 : !! natom=number of atoms in cell.
66 : !! nsym=number of symmetries in space group
67 : !! ntypat=number of atom types
68 : !! paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh
69 : !! pawang <type(pawang_type)>=paw angular mesh and related data
70 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
71 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
72 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
73 : !! rprimd(3,3), conversion from crystal coordinates to cartesian coordinates
74 : !! symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
75 : !! tnons(3,nsym) = nonsymmorphic translations
76 : !! xred(3,natom), location of atoms in crystal coordinates.
77 : !!
78 : !! OUTPUT
79 : !! efg(3,3,natom), the 3x3 efg tensor at each site due to nhat
80 :
81 : !! NOTES
82 : !! This routine computes the electric field gradient, specifically the components
83 : !! $\partial^2 V/\partial x_\alpha \partial x_\beta$ of the potential generated by the valence
84 : !! electrons, at each atomic site in the unit cell. Key references: Kresse and Joubert, ``From
85 : !! ultrasoft pseudopotentials to the projector augmented wave method'', Phys. Rev. B. 59, 1758--1775 (1999) [[cite:Kresse1999]],
86 : !! and Profeta, Mauri, and Pickard, ``Accurate first principles prediction of $^{17}$O NMR parameters in
87 : !! SiO$_2$: Assignment of the zeolite ferrierite spectrum'', J. Am. Chem. Soc. 125, 541--548 (2003) [[cite:Profeta2003]]. See in particular
88 : !! Eq. 11 and 12 of Profeta et al., but note that their sum over occupied states times 2 for occupation number is
89 : !! replaced in the Kresse and Joubert formulation by the sum over $\rho_{ij}$ occupations for each basis element pair.
90 : !!
91 : !! SOURCE
92 :
93 18 : subroutine make_efg_onsite(efg,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab, &
94 9 : & rprimd,symrel,tnons,xred,&
95 9 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
96 :
97 : !Arguments ------------------------------------
98 : !scalars
99 : integer,intent(in) :: my_natom,natom,nsym,ntypat
100 : integer,optional,intent(in) :: comm_atom
101 : type(pawang_type),intent(in) :: pawang
102 : !arrays
103 : integer,intent(in) :: symrel(3,3,nsym)
104 : integer,optional,target,intent(in) :: mpi_atmtab(:)
105 : real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym),xred(3,natom)
106 : real(dp),intent(out) :: efg(3,3,natom)
107 : type(paw_an_type),intent(in) :: paw_an(my_natom)
108 : type(pawrad_type),intent(in) :: pawrad(ntypat)
109 : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
110 : type(pawtab_type),intent(in) :: pawtab(ntypat)
111 :
112 : !Local variables-------------------------------
113 : !scalars
114 : integer :: cplex,iatom,iatom_tot,ictr,ierr,imesh_size,ispden,itypat
115 : integer :: lm,lm_size,local_paw_print_vol
116 : integer :: mesh_size,my_comm_atom,nzlmopt,nspden
117 : integer :: opt_compch,opt_dens,opt_l,opt_print
118 : logical :: my_atmtab_allocated,paral_atom
119 : real(dp) :: c1,c2,c3,compch_sph,intg
120 : !arrays
121 9 : integer,pointer :: my_atmtab(:)
122 9 : logical,allocatable :: lmselectin(:),lmselectout(:)
123 9 : real(dp),allocatable :: ff(:),nhat1(:,:,:),rho1(:,:,:),trho1(:,:,:)
124 :
125 : ! ************************************************************************
126 :
127 : DBG_ENTER("COLL")
128 :
129 9 : if (my_natom>0) then
130 9 : ABI_CHECK(pawrhoij(1)%qphase==1,'make_efg_onsite: not supposed to be called with qphqse=2!')
131 : end if
132 :
133 542 : efg(:,:,:) = zero
134 :
135 : !Set up parallelism over atoms
136 9 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
137 9 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
138 9 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
139 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
140 9 : & my_natom_ref=my_natom)
141 :
142 : !the following factors arise in expanding the angular dependence of the electric field gradient tensor in
143 : !terms of real spherical harmonics. The real spherical harmonics are as in the routine initylmr.F90; see
144 : !in particular also http://www.unioviedo.es/qcg/art/Theochem419-19-ov-BF97-rotation-matrices.pdf
145 9 : c1 = sqrt(16.0*pi/5.0)
146 9 : c2 = sqrt(4.0*pi/5.0)
147 9 : c3 = sqrt(12.0*pi/5.0)
148 :
149 : !loop over atoms in cell
150 50 : do iatom = 1, my_natom
151 41 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
152 41 : itypat=pawrhoij(iatom)%itypat
153 :
154 41 : lm_size = paw_an(iatom)%lm_size
155 41 : if (lm_size < 5) cycle ! if lm_size < 5 then the on-site densities for this atom have no L=2 component
156 : ! and therefore nothing to contribute to the on-site electric field gradient
157 :
158 41 : mesh_size=pawtab(itypat)%mesh_size
159 123 : ABI_MALLOC(ff,(mesh_size))
160 :
161 41 : cplex = pawrhoij(iatom)%qphase
162 41 : nspden = pawrhoij(iatom)%nspden
163 123 : ABI_MALLOC(lmselectin,(lm_size))
164 82 : ABI_MALLOC(lmselectout,(lm_size))
165 490 : lmselectin = .true. ! compute all moments of densities
166 41 : nzlmopt = -1
167 41 : opt_compch = 0
168 41 : compch_sph = zero
169 41 : opt_dens = 0 ! compute all densities
170 41 : opt_l = -1 ! all moments contribute
171 41 : opt_print = 0 ! do not print out moments
172 41 : local_paw_print_vol = 0 ! standard amount of printing in pawdensities
173 :
174 205 : ABI_MALLOC(nhat1,(cplex*mesh_size,lm_size,nspden*(1-((opt_dens+1)/2))))
175 164 : ABI_MALLOC(rho1,(cplex*mesh_size,lm_size,nspden))
176 164 : ABI_MALLOC(trho1,(cplex*mesh_size,lm_size,nspden*(1-((opt_dens+1)/2))))
177 :
178 : ! construct multipole expansion of on-site charge densities for this atom
179 : call pawdensities(compch_sph,cplex,iatom_tot,lmselectin,lmselectout,lm_size,&
180 : & nhat1,nspden,nzlmopt,opt_compch,opt_dens,opt_l,opt_print,&
181 : & pawang,local_paw_print_vol,pawrad(itypat),pawrhoij(iatom),pawtab(itypat),&
182 41 : & rho1,trho1)
183 :
184 : ! spin components:
185 : ! nspden(1) contains total in all cases
186 41 : ispden = 1
187 :
188 246 : do lm = 5, 9 ! loop on L=2 components of multipole expansion
189 :
190 205 : if(.not. lmselectout(lm)) cycle ! skip moments that contributes zero
191 :
192 : ! the following is r^2*(n1-tn1-nhat)/r^3 for this multipole moment
193 : ! use only the real part of the density in case of cplex==2
194 78327 : do imesh_size = 2, mesh_size
195 78206 : ictr = cplex*(imesh_size - 1) + 1
196 78206 : ff(imesh_size)=rho1(ictr,lm,ispden)-trho1(ictr,lm,ispden)-nhat1(ictr,lm,ispden)
197 78327 : ff(imesh_size)=ff(imesh_size)/pawrad(itypat)%rad(imesh_size)
198 : end do
199 121 : call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
200 121 : call simp_gen(intg,ff,pawrad(itypat))
201 162 : select case (lm)
202 : case (5) ! S_{2,-2}
203 34 : efg(1,2,iatom_tot) = efg(1,2,iatom_tot) - c3*intg ! xy case
204 : case (6) ! S_{2,-1}
205 4 : efg(2,3,iatom_tot) = efg(2,3,iatom_tot) - c3*intg ! yz case
206 : case (7) ! S_{2, 0}
207 41 : efg(1,1,iatom_tot) = efg(1,1,iatom_tot) + c2*intg ! xx case
208 41 : efg(2,2,iatom_tot) = efg(2,2,iatom_tot) + c2*intg ! yy case
209 41 : efg(3,3,iatom_tot) = efg(3,3,iatom_tot) - c1*intg ! zz case
210 : case (8) ! S_{2,+1}
211 4 : efg(1,3,iatom_tot) = efg(1,3,iatom_tot) - c3*intg ! xz case
212 : case (9) ! S_{2,+2}
213 38 : efg(1,1,iatom_tot) = efg(1,1,iatom_tot) - c3*intg ! xx case
214 121 : efg(2,2,iatom_tot) = efg(2,2,iatom_tot) + c3*intg ! yy case
215 : end select
216 :
217 : end do ! end loop over LM components with L=2
218 :
219 :
220 : ! Symmetrization of EFG
221 41 : efg(2,1,iatom_tot) = efg(1,2,iatom_tot)
222 41 : efg(3,1,iatom_tot) = efg(1,3,iatom_tot)
223 41 : efg(3,2,iatom_tot) = efg(2,3,iatom_tot)
224 :
225 41 : ABI_FREE(lmselectin)
226 41 : ABI_FREE(lmselectout)
227 41 : ABI_FREE(ff)
228 41 : ABI_FREE(nhat1)
229 41 : ABI_FREE(rho1)
230 50 : ABI_FREE(trho1)
231 :
232 : end do ! Loop on atoms
233 :
234 : !Reduction in case of parallelisation over atoms
235 9 : if (paral_atom) then
236 0 : call xmpi_sum(efg,my_comm_atom,ierr)
237 : end if
238 :
239 : ! symmetrize tensor at each atomic site using point symmetry operations
240 50 : do iatom = 1, natom
241 50 : call matpointsym(iatom,efg(:,:,iatom),natom,nsym,rprimd,symrel,tnons,xred)
242 : end do
243 :
244 : !Destroy atom table used for parallelism
245 9 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
246 :
247 : DBG_EXIT("COLL")
248 :
249 18 : end subroutine make_efg_onsite
250 : !!***
251 :
252 : !----------------------------------------------------------------------
253 :
254 : !!****f* m_paw_nmr/make_fc_paw
255 : !! NAME
256 : !! make_fc_paw
257 : !!
258 : !! FUNCTION
259 : !! Compute the Fermi-contact term due to the PAW cores
260 : !!
261 : !! INPUTS
262 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
263 : !! comm_atom=--optional-- MPI communicator over atoms
264 : !! my_natom=number of atoms treated by current processor
265 : !! natom=number of atoms in cell.
266 : !! nspden=number of spin ensity component
267 : !! ntypat=number of atom types
268 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
269 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
270 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
271 : !!
272 : !! OUTPUT
273 : !! fc(nspden,natom)=the Fermi-contact interaction at each site due to PAW for each spin density
274 : !!
275 : !! NOTES
276 : !! The Fermi contact interaction is the electron density evaluated exactly at the nuclear site.
277 : !! For a nuclear site at R, we are thus computing the expectation value of $\delta^3(R)$, the
278 : !! the three-dimensional delta function at vector position $R$. In terms of the radial variable only
279 : !! the delta function is $\delta(r)/4\pi r^2$. Because this observable is
280 : !! absolutely confined within the PAW radius, only the response due to the AE PAW functions is
281 : !! needed, the pseudo wavefunctions and pseudo PAW functions cancel each other out. We then
282 : !! must compute the integral of $u_i/r times u_j/r \delta(R)d^3r$, for the $l=0$ angular momentum
283 : !! states only. This is simplified with the use of L'H\^{o}spital's theorem to take the limit
284 : !! as $r\rightarrow 0$, yielding $u_i'(r) u_j'(r)$. To compute the derivatives we just fit the
285 : !! first 5 points of the $u$ functions to a line through the origin, using the least squares
286 : !! procedure resulting from $\chi = sum_i (y_i - m*x_i)^2$ . This is more stable than
287 : !! computing the derivative of the whole function and extrapolating it to zero.
288 : !! See Zwanziger, J. Phys. Conden. Matt. 21, 15024-15036 (2009) [[cite:Zwanziger2009]].
289 : !!
290 : !! SOURCE
291 :
292 10 : subroutine make_fc_paw(fc,my_natom,natom,nspden,ntypat,pawrhoij,pawrad,pawtab,&
293 5 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
294 :
295 : !Arguments ------------------------------------
296 : !scalars
297 : integer,intent(in) :: my_natom,natom,nspden,ntypat
298 : integer,optional,intent(in) :: comm_atom
299 : !arrays
300 : integer,optional,target,intent(in) :: mpi_atmtab(:)
301 : real(dp),intent(out) :: fc(nspden,natom)
302 : type(pawrad_type),intent(in) :: pawrad(ntypat)
303 : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
304 : type(pawtab_type),target,intent(in) :: pawtab(ntypat)
305 :
306 : !Local variables-------------------------------
307 : !scalars
308 : integer :: iatom,iatom_tot,ierr,irhoij,islope,ispden,itypat
309 : integer :: ilmn,il,iln,ilm,im,jl,jlm,jlmn,jln,jm,j0lmn,jrhoij
310 : integer :: klmn,kln,mesh_size,my_comm_atom,nslope
311 : logical :: my_atmtab_allocated,paral_atom
312 : real(dp) :: mi,mj,xi,xxsum,xysumi,xysumj,yi,yj
313 : !arrays
314 5 : integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
315 5 : integer,pointer :: my_atmtab(:)
316 :
317 : ! ************************************************************************
318 :
319 : DBG_ENTER("COLL")
320 :
321 5 : if (my_natom>0) then
322 5 : ABI_CHECK(pawrhoij(1)%qphase==1,'make_fc_paw: not supposed to be called with qphqse=2!')
323 : end if
324 :
325 : !Set up parallelism over atoms
326 5 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
327 5 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
328 5 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
329 5 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
330 :
331 : !Initialization
332 48 : fc(:,:)=zero
333 :
334 : !number of points to use in computing initial slopes of radial functions
335 : nslope = 5
336 :
337 : !loop over atoms in cell
338 26 : do iatom = 1, my_natom
339 21 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
340 21 : itypat = pawrhoij(iatom)%itypat
341 21 : mesh_size=pawtab(itypat)%mesh_size
342 21 : indlmn => pawtab(itypat)%indlmn
343 :
344 : ! loop over spin components
345 48 : do ispden=1,nspden
346 :
347 : ! loop over basis elements for this atom
348 : ! ----
349 259 : do jlmn=1,pawtab(itypat)%lmn_size
350 216 : jl=indlmn(1,jlmn)
351 216 : jm=indlmn(2,jlmn)
352 216 : jlm=indlmn(4,jlmn)
353 216 : jln=indlmn(5,jlmn)
354 216 : j0lmn=jlmn*(jlmn-1)/2
355 1520 : do ilmn=1,jlmn
356 1282 : il=indlmn(1,ilmn)
357 1282 : im=indlmn(2,ilmn)
358 1282 : iln=indlmn(5,ilmn)
359 1282 : ilm=indlmn(4,ilmn)
360 1282 : klmn=j0lmn+ilmn
361 1282 : kln = pawtab(itypat)%indklmn(2,klmn) ! need this for mesh selection below
362 :
363 1498 : if (jl==0 .and. il==0) then ! select only s-states
364 :
365 : ! Loop over non-zero elements of rhoij
366 66 : jrhoij=1
367 1584 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
368 1518 : if (klmn==pawrhoij(iatom)%rhoijselect(irhoij)) then ! rho_ij /= 0 for this klmn
369 : xxsum = 0 ! these three variables will be used to compute the slopes
370 : xysumi = 0
371 : xysumj = 0
372 396 : do islope=1, nslope
373 330 : xi=0
374 330 : if(pawrad(itypat)%mesh_type == 1) xi = (islope - 1)*pawrad(itypat)%rstep
375 330 : if(pawrad(itypat)%mesh_type == 2) xi = pawrad(itypat)%rstep * &
376 330 : & (exp(pawrad(itypat)%lstep * (islope - 1)) - 1)
377 330 : if(pawrad(itypat)%mesh_type == 3) then
378 0 : if (islope == 1) then
379 : xi = 0
380 : else
381 0 : xi = pawrad(itypat)%rstep * exp(pawrad(itypat)%lstep*(islope-1))
382 : end if
383 : end if
384 330 : if(pawrad(itypat)%mesh_type == 4) xi = &
385 0 : & -pawrad(itypat)%rstep*log(1.0-(islope-1)/pawrad(itypat)%mesh_size)
386 330 : yi = pawtab(itypat)%phi(islope,iln) ! function value for u_i
387 330 : yj = pawtab(itypat)%phi(islope,jln) ! function value for u_j
388 330 : xxsum = xxsum + xi*xi
389 330 : xysumi = xysumi + xi*yi
390 396 : xysumj = xysumj + xi*yj
391 : end do
392 : ! the slopes of the radial functions are obtained by minimizing
393 : ! chi = sum(y_i - m*x_i)^2 (in other words, a linear least squares
394 : ! fit constrained to go through the origin)
395 : ! the result is m = sum(y_i*x_i)/sum(x_i*x_i)
396 66 : mi = xysumi/xxsum
397 66 : mj = xysumj/xxsum
398 : ! accumulate the rho_ij contribution to the fermi contact for this spin density:
399 : fc(ispden,iatom_tot)=fc(ispden,iatom_tot)+&
400 66 : & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(jrhoij,ispden)*mi*mj/four_pi
401 : end if ! end selection on klmn for nonzero rho_ij
402 1584 : jrhoij=jrhoij+pawrhoij(iatom)%cplex_rhoij
403 : end do ! end loop over nonzero rho_ij
404 : end if ! end l=l'=0 selection
405 : end do ! end loop over ilmn
406 : end do ! end loop over jlmn
407 : end do ! end loop over spin densities
408 : end do ! Loop on atoms
409 :
410 : !Reduction in case of parallelisation over atoms
411 5 : if (paral_atom) then
412 0 : call xmpi_sum(fc,my_comm_atom,ierr)
413 : end if
414 :
415 : !Destroy atom table used for parallelism
416 5 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
417 :
418 : DBG_EXIT("COLL")
419 :
420 5 : end subroutine make_fc_paw
421 : !!***
422 :
423 : !----------------------------------------------------------------------
424 :
425 : END MODULE m_paw_nmr
426 : !!***
|
