Line data Source code
1 : !!****m* ABINIT/m_nucprop
2 : !! NAME
3 : !! m_nucprop
4 : !!
5 : !! FUNCTION
6 : !! routines used to compute properties at the nuclear sites, including
7 : !! electric field gradient and Fermi contact
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 1998-2026 ABINIT group (MT, JWZ)
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_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_matrix, only : matr3inv
43 : use m_symtk, only : matpointsym
44 : use m_fft, only : fourdp
45 :
46 : implicit none
47 :
48 : private
49 : !!***
50 :
51 : public :: calc_efg
52 : public :: calc_fc
53 : public :: make_efg_ion
54 : public :: make_efg_el
55 : !!***
56 :
57 : contains
58 :
59 : !!****f* ABINIT/calc_efg
60 : !! NAME
61 : !! calc_efg
62 : !!
63 : !! FUNCTION
64 : !! calculation and output of electric field gradient tensor at each atomic site
65 : !!
66 : !! INPUTS
67 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
68 : !! comm_atom=--optional-- MPI communicator over atoms
69 : !! mpi_enreg=information about MPI parallelization
70 : !! my_natom=number of atoms treated by current processor
71 : !! natom=number of atoms in cell.
72 : !! nfft=number of points on fft grid
73 : !! ngfft(18)=details of fft
74 : !! nhat(nfft,nspden)=compensation charge density
75 : !! nspden=number of spin densities
76 : !! nsym=number of symmetries in space group
77 : !! ntypat=number of atom types
78 : !! nucefg=1 to print summary output, 2 for detailed output
79 : !! ptcharge(ntypat)=user input charges on atoms to make simple point charge calc
80 : !! paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh
81 : !! pawang <type(pawang_type)>=paw angular mesh and related data
82 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
83 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
84 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
85 : !! quadmom(ntypat)=quadrupole moments in barns of different atomic nuclei
86 : !! rhor(nfft,nspden)=electron density on grid (strictly $\tilde{n}+\hat{n}$)
87 : !! rprimd(3,3)=matrix relating cartesian coordinates to crystal coordinates
88 : !! symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
89 : !! tnons(3,nsym)=nonsymmorphic translations
90 : !! typat(natom)=type (integer) for each atom
91 : !! ucvol=unit cell volume in Bohr^3
92 : !! usepaw=1 if we are using PAW formalism, 0 else
93 : !! xred(3,natom)=vectors locating each atom in the unit cell, in crystal coords
94 : !! zion(ntypat)=net core charge on each type of atom
95 : !!
96 : !! OUTPUT
97 : !! (only writing, printing)
98 : !!
99 : !! SOURCE
100 :
101 18 : subroutine calc_efg(efg,mpi_enreg,my_natom,natom,nfft,ngfft,nhat,nspden,nsym,nucefg,ntypat,&
102 9 : paw_an,pawang,pawrad,pawrhoij,pawtab,&
103 9 : ptcharge,quadmom,rhor,rprimd,symrel,tnons,typat,ucvol,usepaw,xred,zion,&
104 9 : mpi_atmtab,comm_atom) ! optional arguments (parallelism)
105 :
106 : !Arguments ------------------------------------
107 : !scalars
108 : integer,intent(in) :: my_natom,natom,nfft,nspden,nsym,nucefg,ntypat,usepaw
109 : integer,optional,intent(in) :: comm_atom
110 : real(dp),intent(in) :: ucvol
111 : type(MPI_type),intent(in) :: mpi_enreg
112 : type(pawang_type),intent(in) :: pawang
113 : !arrays
114 : integer,intent(in) :: ngfft(18),symrel(3,3,nsym),typat(natom)
115 : integer,optional,target,intent(in) :: mpi_atmtab(:)
116 : real(dp),intent(in) :: nhat(nfft,nspden),ptcharge(ntypat)
117 : real(dp),intent(in) :: quadmom(ntypat),rhor(nfft,nspden),rprimd(3,3)
118 : real(dp),intent(in) :: tnons(3,nsym),zion(ntypat)
119 : real(dp),intent(inout) :: xred(3,natom)
120 : real(dp),intent(out) :: efg(3,3,natom)
121 : type(paw_an_type),intent(in) :: paw_an(my_natom)
122 : type(pawrad_type),intent(in) :: pawrad(ntypat)
123 : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
124 : type(pawtab_type),intent(in) :: pawtab(ntypat)
125 :
126 : !Local variables-------------------------------
127 : !scalars
128 : integer :: ii,INFO,LDA,LWORK,N,iatom,my_comm_atom
129 : logical :: my_atmtab_allocated,paral_atom
130 : real(dp) :: cq,efgsi21,eta,vxx,vyy,vzz
131 : character(len=500) :: message
132 : !arrays
133 9 : integer,pointer :: my_atmtab(:)
134 : real(dp) :: eigval(3),matr(3,3),work(8)
135 9 : real(dp),allocatable :: efg_el(:,:,:),efg_ion(:,:,:),efg_paw(:,:,:)
136 9 : real(dp),allocatable :: efg_point_charge(:,:,:)
137 :
138 : ! ************************************************************************
139 :
140 : !Compatibility tests
141 9 : if (usepaw /= 1) then
142 0 : message = ' usepaw /= 1 but EFG calculation requires PAW '
143 0 : ABI_ERROR(message)
144 : end if
145 :
146 9 : efgsi21=efg_si*1.0E-21 ! efg_si is electric field gradient in SI units, defined in
147 : ! defs_basis. Multiply by 10E-21 for nice printing below.
148 : ! 1 EFG in au is 9.725E21 volts/m^2
149 :
150 : !Set up parallelism over atoms
151 9 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
152 9 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
153 9 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
154 9 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
155 :
156 27 : ABI_MALLOC(efg_el,(3,3,natom))
157 18 : ABI_MALLOC(efg_ion,(3,3,natom))
158 18 : ABI_MALLOC(efg_paw,(3,3,natom))
159 18 : ABI_MALLOC(efg_point_charge,(3,3,natom))
160 542 : efg_el(:,:,:) = zero
161 542 : efg_ion(:,:,:) = zero
162 542 : efg_paw(:,:,:) = zero
163 542 : efg_point_charge(:,:,:) = zero
164 :
165 9 : call make_efg_el(efg_el,mpi_enreg,natom,nfft,ngfft,nhat,nspden,nsym,rhor,rprimd,symrel,tnons,xred)
166 :
167 9 : call make_efg_ion(efg_ion,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,zion)
168 :
169 9 : if (paral_atom) then
170 : call make_efg_onsite(efg_paw,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab,&
171 0 : & rprimd,symrel,tnons,xred,comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
172 : else
173 : call make_efg_onsite(efg_paw,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab,&
174 9 : & rprimd,symrel,tnons,xred)
175 : end if
176 :
177 : !calculate efg due to pure point charges, as input in variable ptcharge(ntypat)
178 : !note here all atoms of the same type will have the same valence; in the future this
179 : !could be made more flexible by having ptcharge(natom) but that will require a slightly
180 : !different version than the existing make_efg_ion routine
181 9 : if(nucefg > 2) then
182 1 : call make_efg_ion(efg_point_charge,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,ptcharge)
183 : end if
184 :
185 542 : efg(:,:,:) = efg_el(:,:,:) + efg_ion(:,:,:) + efg_paw(:,:,:)
186 :
187 9 : write(message,'(a,a,a)' ) ch10,' Electric Field Gradient Calculation ',ch10
188 9 : call wrtout(ab_out,message,'COLL')
189 :
190 9 : LDA=3; LWORK=8;N=3 ! these parameters are needed for the LAPACK dsyev routine
191 50 : do iatom = 1, natom
192 533 : matr(:,:) = efg(:,:,iatom)
193 41 : call dsyev('V','U',N,matr,LDA,eigval,work,LWORK,INFO) ! get eigenvalues and eigenvectors
194 41 : if (eigval(3) > abs(eigval(1)) ) then ! In NMR, the convention is that whatever component is
195 : ! largest in magnitude is called Vzz, next comes Vxx, then Vyy
196 16 : vzz = eigval(3)
197 16 : vxx = eigval(1)
198 16 : vyy = eigval(2)
199 : else
200 25 : vzz = eigval(1)
201 25 : vxx = eigval(3)
202 25 : vyy = eigval(2)
203 : end if
204 : ! Cq = vzz*(eQ)/h, where Q is the nuclear quad moment in barns (10E-28 m2)
205 : ! Multiply Q by 1E-8 * Ang_Bohr**2 to get nuclear moment in Bohr^2
206 : ! resulting vzz*(eQ) is energy in Ha (recall e = 1 in au)
207 : ! then convert to MHz with Ha_THz*1.E6
208 : ! resulting factors are Ang_Bohr**2*Ha_THz*1.0D-2
209 41 : cq = vzz*quadmom(typat(iatom))*Ang_Bohr**2*Ha_THz*1.0D-2
210 41 : if (abs(cq) > tol8) then
211 29 : eta = abs(vxx-vyy)/abs(vzz)
212 : else
213 12 : cq = zero
214 12 : eta = zero ! if Cq is small then eta is meaningless
215 : end if
216 :
217 41 : write(message,'(a,a,i4,a,i4)')ch10,' atom : ',iatom,' typat : ',typat(iatom)
218 41 : call wrtout(ab_out,message,'COLL')
219 41 : if (nucefg > 1) then
220 41 : write(message,'(2a,f9.4,a,f9.4,a,f9.4)') ch10,' Nuclear quad. mom. (barns) : ',quadmom(typat(iatom)),&
221 82 : & ' Cq (MHz) : ',cq,' eta : ',eta
222 41 : call wrtout(ab_out,message,'COLL')
223 : end if
224 :
225 : ! for printing and test portability, it's better to simply set very small eigvals to zero
226 164 : do ii=1,3
227 164 : if (abs(eigval(ii))<tol8) eigval(ii)=zero
228 : end do
229 41 : write(message,'(2a,f13.6,a,f16.8,a,a,3f13.6)')ch10,' efg eigval (au) : ',eigval(1),' ; (1.0E+21 V/m^2) : ',eigval(1)*efgsi21,ch10,&
230 82 : & '- eigvec : ',matr(1,1),matr(2,1),matr(3,1)
231 41 : call wrtout(ab_out,message,'COLL')
232 41 : write(message,'(2a,f13.6,a,f16.8,a,a,3f13.6)')ch10,' efg eigval (au) : ',eigval(2),' ; (1.0E+21 V/m^2) : ',eigval(2)*efgsi21,ch10,&
233 82 : & '- eigvec : ',matr(1,2),matr(2,2),matr(3,2)
234 41 : call wrtout(ab_out,message,'COLL')
235 41 : write(message,'(2a,f13.6,a,f16.8,a,a,3f13.6)')ch10,' efg eigval (au) : ',eigval(3),' ; (1.0E+21 V/m^2) : ',eigval(3)*efgsi21,ch10,&
236 82 : & '- eigvec : ',matr(1,3),matr(2,3),matr(3,3)
237 41 : call wrtout(ab_out,message,'COLL')
238 41 : write(message,'(a,a,3f13.6)')ch10,' total efg : ',efg(1,1,iatom),efg(1,2,iatom),efg(1,3,iatom)
239 41 : call wrtout(ab_out,message,'COLL')
240 41 : write(message,'(a,3f13.6)')' total efg : ',efg(2,1,iatom),efg(2,2,iatom),efg(2,3,iatom)
241 41 : call wrtout(ab_out,message,'COLL')
242 41 : write(message,'(a,3f13.6,a)')' total efg : ',efg(3,1,iatom),efg(3,2,iatom),efg(3,3,iatom),ch10
243 41 : call wrtout(ab_out,message,'COLL')
244 41 : write(message,'(a,a,3f13.6)')ch10,' efg_el : ',efg_el(1,1,iatom),efg_el(1,2,iatom),efg_el(1,3,iatom)
245 41 : call wrtout(ab_out,message,'COLL')
246 41 : write(message,'(a,3f13.6)')' efg_el : ',efg_el(2,1,iatom),efg_el(2,2,iatom),efg_el(2,3,iatom)
247 41 : call wrtout(ab_out,message,'COLL')
248 41 : write(message,'(a,3f13.6,a)')' efg_el : ',efg_el(3,1,iatom),efg_el(3,2,iatom),efg_el(3,3,iatom),ch10
249 41 : call wrtout(ab_out,message,'COLL')
250 41 : write(message,'(a,3f13.6)')' efg_ion : ',efg_ion(1,1,iatom),efg_ion(1,2,iatom),efg_ion(1,3,iatom)
251 41 : call wrtout(ab_out,message,'COLL')
252 41 : write(message,'(a,3f13.6)')' efg_ion : ',efg_ion(2,1,iatom),efg_ion(2,2,iatom),efg_ion(2,3,iatom)
253 41 : call wrtout(ab_out,message,'COLL')
254 41 : write(message,'(a,3f13.6,a)')' efg_ion : ',efg_ion(3,1,iatom),efg_ion(3,2,iatom),efg_ion(3,3,iatom),ch10
255 41 : call wrtout(ab_out,message,'COLL')
256 41 : write(message,'(a,3f13.6)')' efg_paw : ',efg_paw(1,1,iatom),efg_paw(1,2,iatom),efg_paw(1,3,iatom)
257 41 : call wrtout(ab_out,message,'COLL')
258 41 : write(message,'(a,3f13.6)')' efg_paw : ',efg_paw(2,1,iatom),efg_paw(2,2,iatom),efg_paw(2,3,iatom)
259 41 : call wrtout(ab_out,message,'COLL')
260 41 : write(message,'(a,3f13.6,a)')' efg_paw : ',efg_paw(3,1,iatom),efg_paw(3,2,iatom),efg_paw(3,3,iatom),ch10
261 41 : call wrtout(ab_out,message,'COLL')
262 50 : if (nucefg > 2) then ! write output of pure pointcharge calculation
263 78 : matr(:,:) = efg_point_charge(:,:,iatom)
264 6 : call dsyev('V','U',N,matr,LDA,eigval,work,LWORK,INFO) ! get eigenvalues and eigenvectors
265 6 : if (eigval(3) > abs(eigval(1)) ) then ! In NMR, the convention is that whatever component is
266 : ! largest in magnitude is called Vzz, next comes Vxx, then Vyy
267 0 : vzz = eigval(3)
268 0 : vxx = eigval(1)
269 0 : vyy = eigval(2)
270 : else
271 6 : vzz = eigval(1)
272 6 : vxx = eigval(3)
273 6 : vyy = eigval(2)
274 : end if
275 6 : cq = vzz*quadmom(typat(iatom))*Ang_Bohr**2*Ha_THz*1.0D-2
276 6 : if (abs(cq) > tol8) then
277 4 : eta = abs(vxx-vyy)/abs(vzz)
278 : else
279 2 : eta = zero ! if Cq is small then eta is meaningless
280 : end if
281 6 : write(message,'(a,f9.4,a,f9.4)') ' Point charge Cq = ',cq,' MHz eta = ',eta
282 6 : call wrtout(ab_out,message,'COLL')
283 6 : write(message,'(2a,f13.6,a,f16.8,a,a,3f13.6)')ch10,' point charge eigval (au) : ',&
284 6 : & eigval(1),' ; (1.0E+21 V/m^2) : ',eigval(1)*efgsi21,ch10,&
285 12 : & '- eigvec : ',matr(1,1),matr(2,1),matr(3,1)
286 6 : call wrtout(ab_out,message,'COLL')
287 6 : write(message,'(2a,f13.6,a,f16.8,a,a,3f13.6)')ch10,' point charge eigval (au) : ',&
288 6 : & eigval(2),' ; (1.0E+21 V/m^2) : ',eigval(2)*efgsi21,ch10,&
289 12 : & '- eigvec : ',matr(1,2),matr(2,2),matr(3,2)
290 6 : call wrtout(ab_out,message,'COLL')
291 6 : write(message,'(2a,f13.6,a,f16.8,a,a,3f13.6)')ch10,' point charge eigval (au) : ',&
292 6 : & eigval(3),' ; (1.0E+21 V/m^2) : ',eigval(3)*efgsi21,ch10,&
293 12 : & '- eigvec : ',matr(1,3),matr(2,3),matr(3,3)
294 6 : call wrtout(ab_out,message,'COLL')
295 6 : write(message,'(a,a,3f13.6)')ch10,' point charge efg : ',efg_point_charge(1,1,iatom),&
296 12 : & efg_point_charge(1,2,iatom),efg_point_charge(1,3,iatom)
297 6 : call wrtout(ab_out,message,'COLL')
298 6 : write(message,'(a,3f13.6)')' point charge efg : ',efg_point_charge(2,1,iatom),&
299 12 : & efg_point_charge(2,2,iatom),efg_point_charge(2,3,iatom)
300 6 : call wrtout(ab_out,message,'COLL')
301 6 : write(message,'(a,3f13.6,a)')' point charge efg : ',efg_point_charge(3,1,iatom),&
302 12 : & efg_point_charge(3,2,iatom),efg_point_charge(3,3,iatom),ch10
303 6 : call wrtout(ab_out,message,'COLL')
304 : end if
305 : end do
306 9 : write(message,'(3a)')ch10,ch10,ch10
307 9 : call wrtout(ab_out,message,'COLL')
308 :
309 9 : ABI_FREE(efg_el)
310 9 : ABI_FREE(efg_ion)
311 9 : ABI_FREE(efg_paw)
312 9 : ABI_FREE(efg_point_charge)
313 :
314 : !Destroy atom table used for parallelism
315 9 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
316 :
317 : !DEBUG
318 : !write(std_out,*)' calc_efg : exit '
319 : !stop
320 : !ENDDEBUG
321 :
322 9 : end subroutine calc_efg
323 : !!***
324 :
325 : !!***
326 : !!****f* ABINIT/calc_fc
327 : !! NAME
328 : !! calc_fc
329 : !!
330 : !! FUNCTION
331 : !! calculation and output of Fermi-contact term at each atomic site
332 : !!
333 : !! INPUTS
334 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
335 : !! comm_atom=--optional-- MPI communicator over atoms
336 : !! my_natom=number of atoms treated by current processor
337 : !! natom=number of atoms in cell.
338 : !! nspden=number of spin density components
339 : !! ntypat=number of atom types
340 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
341 : !! pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
342 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
343 : !! typat(natom)=type (integer) for each atom
344 : !! usepaw=1 if PAW is activated
345 : !!
346 : !! OUTPUT
347 : !! (only writing, printing)
348 : !!
349 : !! SIDE EFFECTS
350 : !!
351 : !!
352 : !! NOTES
353 : !!
354 : !! SOURCE
355 :
356 10 : subroutine calc_fc(my_natom,natom,nspden,ntypat,pawrad,pawrhoij,pawtab,typat,usepaw,&
357 5 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
358 :
359 : !Arguments ------------------------------------
360 : !scalars
361 : integer,intent(in) :: my_natom,natom,nspden,ntypat,usepaw
362 : integer,optional,intent(in) :: comm_atom
363 : !arrays
364 : integer,intent(in) :: typat(natom)
365 : integer,optional,target,intent(in) :: mpi_atmtab(:)
366 : type(pawrad_type),intent(in) :: pawrad(ntypat)
367 : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
368 : type(pawtab_type),intent(in) :: pawtab(ntypat)
369 :
370 : !Local variables-------------------------------
371 : !scalars
372 : integer :: iatom,my_comm_atom
373 : logical :: my_atmtab_allocated,paral_atom
374 : character(len=500) :: message
375 : !arrays
376 5 : integer,pointer :: my_atmtab(:)
377 5 : real(dp),allocatable :: fc(:,:)
378 :
379 : !***********************************************************************
380 :
381 : !Compatibility tests
382 5 : if (usepaw /= 1) then
383 0 : message = ' usepaw /= 1 but Fermi-contact calculation requires PAW '
384 0 : ABI_ERROR(message)
385 : end if
386 :
387 : !Set up parallelism over atoms
388 5 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
389 5 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
390 5 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
391 5 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
392 :
393 : !Initialization
394 20 : ABI_MALLOC(fc,(nspden,natom))
395 :
396 : !Computation
397 5 : if (paral_atom) then
398 : call make_fc_paw(fc,my_natom,natom,nspden,ntypat,pawrhoij,pawrad,pawtab,&
399 0 : & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
400 : else
401 5 : call make_fc_paw(fc,my_natom,natom,nspden,ntypat,pawrhoij,pawrad,pawtab)
402 : end if
403 :
404 : !Printing
405 5 : write(message,'(a,a,a)' ) ch10,' Fermi-contact Term Calculation ',ch10
406 5 : call wrtout(ab_out,message,'COLL')
407 :
408 26 : do iatom = 1, natom
409 26 : if (nspden == 2) then
410 1 : write(message,'(a,i3,a,i3,a,f12.4)') ' Atom ',iatom,', typat ',typat(iatom),': FC total = ',&
411 2 : & fc(1,iatom)+fc(2,iatom)
412 1 : call wrtout(ab_out,message,'COLL')
413 1 : write(message,'(a,i3,a,i3,a,f12.4)') ' Atom ',iatom,', typat ',typat(iatom),': FC up - down = ',&
414 2 : & fc(1,iatom)-fc(2,iatom)
415 1 : call wrtout(ab_out,message,'COLL')
416 : else
417 20 : write(message,'(a,i3,a,i3,a,f12.4)') ' Atom ',iatom,', typat ',typat(iatom),': FC = ',&
418 40 : & fc(1,iatom)
419 20 : call wrtout(ab_out,message,'COLL')
420 : end if
421 : end do
422 :
423 5 : write(message,'(3a)')ch10,ch10,ch10
424 5 : call wrtout(ab_out,message,'COLL')
425 :
426 : !Memory deallocation
427 5 : ABI_FREE(fc)
428 :
429 : !Destroy atom table used for parallelism
430 5 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
431 :
432 5 : end subroutine calc_fc
433 : !!***
434 :
435 : !!****f* ABINIT/make_efg_ion
436 : !! NAME
437 : !! make_efg_ion
438 : !!
439 : !! FUNCTION
440 : !! compute the electric field gradient due to ionic cores
441 : !!
442 : !! INPUTS
443 : !! natom, number of atoms in the unit cell
444 : !! nsym=number of symmetries in space group
445 : !! ntypat, the number of types of atoms in the unit cell
446 : !! rprimd(3,3), the matrix giving the transformation from crystal to cartesian coordinates
447 : !! symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
448 : !! tnons(3,nsym) = nonsymmorphic translations
449 : !! typat(natom), the type of each atom in the unit cell
450 : !! ucvol, the volume of the unit cell in atomic units
451 : !! xred(3,natom) the location of each atom in the cell in crystallographic coordinates
452 : !! zion(ntypat) the net charge on each type of atom
453 : !!
454 : !! OUTPUT
455 : !! efg(3,3,natom), the 3x3 efg tensors at each atomic site
456 : !!
457 : !! SIDE EFFECTS
458 : !!
459 : !! NOTES
460 : !! This routine computes the electric field gradient, specifically the components
461 : !! $\partial^2 V/\partial x_\alpha \partial x_\beta$ of the potential generated by the ionic cores,
462 : !! at each atomic site in the unit cell.
463 : !! Key references:
464 : !! Profeta, Mauri, and Pickard, ``Accurate first principles prediction of $^{17}$O NMR parameters in
465 : !! SiO$_2$: Assignment of the zeolite ferrierite spectrum'', J. Am. Chem. Soc. 125, 541--548 (2003) [[cite:Profeta2003]];
466 : !! A. Honma, ``Dipolar lattice-sums with applications to the exciton bands of anthracene crystal and
467 : !! the crystal field due to point charges'', J. Phys. Soc. Jpn. 42, 1129--1135 (1977) [[cite:Honma1977]];
468 : !! and Kresse and Joubert, ``From ultrasoft pseudopotentials to the projector augmented wave method'',
469 : !! Phys. Rev. B. 59, 1758--1775 (1999) [[cite:Kresse1999]]. In Kresse and Joubert's notation, the ionic cores are $n_{Zc}$;
470 : !! these charges are given by the net core charges on the pseudoatoms. Due to otherwise slow convergence,
471 : !! the sum over atoms is carried out by an Ewald method as detailed in the Honma reference, specifically
472 : !! his Eq. 4.8.
473 : !!
474 : !! SOURCE
475 :
476 10 : subroutine make_efg_ion(efg,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,zion)
477 :
478 : !Arguments ------------------------------------
479 : !scalars
480 : integer,intent(in) :: natom,nsym,ntypat
481 : real(dp) :: ucvol
482 : !arrays
483 : integer,intent(in) :: symrel(3,3,nsym),typat(natom)
484 : real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym)
485 : real(dp),intent(in) :: zion(ntypat)
486 : real(dp),intent(inout) :: xred(3,natom)
487 : real(dp),intent(out) :: efg(3,3,natom)
488 : !Local variables-------------------------------
489 : !scalars
490 : integer :: iatom,ishell,ii,jatom,jj,nshell,sx,sy,sz
491 : real(dp) :: cph,dampfac,derfc_karg,derivs,gsq,karg
492 : real(dp) :: lenrho,phase,qk,rlkcut,trace,xi0
493 : real(dp) :: glkcut
494 : !arrays
495 : real(dp) :: cvec(3),gvec(3),gpl(3),gprimd(3,3)
496 : real(dp) :: rhok(3),rhored(3),rpl(3)
497 10 : real(dp),allocatable :: efg_g(:,:,:),efg_r(:,:,:)
498 10 : real(dp),allocatable :: xcart(:,:)
499 :
500 : ! ************************************************************************
501 :
502 : !DEBUG
503 : !write(std_out,*)' make_efg_ion : enter'
504 : !ENDDEBUG
505 :
506 30 : ABI_MALLOC(efg_g,(3,3,natom))
507 20 : ABI_MALLOC(efg_r,(3,3,natom))
508 30 : ABI_MALLOC(xcart,(3,natom))
509 621 : efg(:,:,:) = zero ! final efg tensor
510 621 : efg_g(:,:,:) = zero ! part of tensor accumulated in G space
511 621 : efg_r(:,:,:) = zero ! part of tensor accumulated in R space
512 :
513 10 : call xred2xcart(natom,rprimd,xcart,xred) ! get atomic locations in cartesian coords
514 :
515 40 : do ii = 1, 3 ! generate the lengths of the unit cell edges in atomic units
516 40 : rpl(ii) = sqrt(rprimd(1,ii)**2+rprimd(2,ii)**2+rprimd(3,ii)**2)
517 : end do
518 90 : xi0 = sqrt(pi/(maxval(rpl)*minval(rpl))) ! this estimate for xi0 is from Honma's paper
519 :
520 10 : call matr3inv(rprimd,gprimd) ! gprimd holds the inverse transpose of rprimd
521 : !remember ordering: rprimd( (x_comp,y_comp,z_comp), (edge 1, edge 2, edge 3) )
522 : !while gprimd( (edge 1, edge 2, edge 3),(x_comp, y_comp, z_comp) )
523 40 : do ii = 1, 3 ! generate the lengths of the reciprocal cell edges
524 40 : gpl(ii) = sqrt(gprimd(ii,1)**2+gprimd(ii,2)**2+gprimd(ii,3)**2)
525 : end do
526 :
527 : !go out enough shells such that g**2/4*xi0**2 is of order 30
528 50 : nshell = int(anint(sqrt(30.0)*xi0/(pi*minval(gpl))))
529 50 : glkcut = (0.95*nshell*two*pi*minval(gpl))**2
530 :
531 56 : do ishell = 0, nshell ! loop over shells
532 270 : do sx = -ishell, ishell
533 1586 : do sy = -ishell, ishell
534 10874 : do sz = -ishell, ishell
535 10660 : if ( .not. (sx==0 .and. sy==0 .and. sz==0) ) then ! avoid origin
536 : ! constrain to be on shell surface, not interior
537 9288 : if ( abs(sx)==ishell .or. abs(sy)==ishell .or. abs(sz)==ishell ) then
538 5736 : cvec(1)=sx;cvec(2)=sy;cvec(3)=sz
539 : ! make the g vector in cartesian coords
540 5736 : gvec(:) = zero
541 22944 : do ii = 1, 3
542 74568 : do jj = 1, 3
543 68832 : gvec(ii) = gvec(ii) + gprimd(ii,jj)*cvec(jj)*two*pi
544 : end do
545 : end do
546 22944 : gsq = dot_product(gvec,gvec)
547 5736 : if(gsq < glkcut) then
548 1200 : dampfac = exp(-gsq/(4.0*xi0*xi0)) ! see Honma eq. 4.8
549 7364 : do iatom = 1, natom
550 42576 : do jatom = 1, natom
551 35212 : qk = zion(typat(jatom)) ! charge on neighbor atom
552 140848 : rhok = xcart(:,jatom)-xcart(:,iatom)
553 140848 : phase = dot_product(gvec,rhok)
554 35212 : cph = cos(phase)
555 147012 : do ii = 1, 3
556 457756 : do jj = 1, 3
557 316908 : derivs = -3.0*gvec(ii)*gvec(jj)/gsq
558 316908 : if (ii == jj) derivs = 1.0 + derivs
559 : efg_g(ii,jj,iatom) = efg_g(ii,jj,iatom) + &
560 422544 : & qk*cph*derivs*dampfac
561 : end do ! end loop over jj
562 : end do ! end loop over ii
563 : end do ! end loop over jatom
564 : end do ! end loop over iatom
565 : end if ! constrain to gsq < glkcut
566 : end if ! end selection on shell edge
567 : end if ! end avoidance of origin
568 : end do ! end loop over sz
569 : end do ! end loop over sy
570 : end do ! end loop over sx
571 : end do ! end loop over ishell
572 :
573 : !sum in real space begins here
574 :
575 : !go out enough shells such that (r*xi0)**2 is of order 30
576 50 : nshell = int(anint(sqrt(30.)/(minval(rpl)*xi0)))
577 50 : rlkcut = nshell*minval(rpl)*0.95
578 : !
579 : !go out enough shells so that rlkcut is of order 30 bohr
580 : !nshell=int(anint(30.0/minval(rpl)))
581 : !rlkcut = 0.95*nshell*minval(rpl)
582 :
583 56 : do ishell = 0, nshell ! total set of cells to loop over
584 270 : do sx = -ishell, ishell ! loop over all cells in each dimension
585 1586 : do sy = -ishell, ishell
586 10874 : do sz = -ishell, ishell
587 : ! constrain to shell surface, not interior
588 10660 : if ( abs(sx)==ishell .or. abs(sy)==ishell .or. abs(sz)==ishell ) then
589 35763 : do jatom = 1, natom ! loop over atoms in shell cell
590 212250 : do iatom = 1, natom ! loop over atoms in central unit cell
591 202916 : if (.NOT. (jatom == iatom .AND. sx == 0 .AND. sy == 0 .AND. sz == 0)) then ! avoid self term
592 172852 : qk = zion(typat(jatom)) ! charge on each neighbor atom
593 : ! ! rhored is the vector in crystal coords from neighbor to target
594 172852 : rhored(1) = xred(1,jatom) + sx - xred(1,iatom)
595 172852 : rhored(2) = xred(2,jatom) + sy - xred(2,iatom)
596 172852 : rhored(3) = xred(3,jatom) + sz - xred(3,iatom)
597 : ! ! rhok is rhored in cartesian coords
598 172852 : rhok(1) = rprimd(1,1)*rhored(1)+rprimd(1,2)*rhored(2)+rprimd(1,3)*rhored(3)
599 172852 : rhok(2) = rprimd(2,1)*rhored(1)+rprimd(2,2)*rhored(2)+rprimd(2,3)*rhored(3)
600 172852 : rhok(3) = rprimd(3,1)*rhored(1)+rprimd(3,2)*rhored(2)+rprimd(3,3)*rhored(3)
601 691408 : trace = dot_product(rhok,rhok)
602 172852 : lenrho = sqrt(trace)
603 172852 : if (lenrho < rlkcut) then ! this restriction is critical as it ensures
604 : ! ! that we sum over a sphere of atoms in real space
605 : ! ! no matter what shape the unit cell has
606 24896 : karg = xi0*lenrho
607 24896 : derfc_karg = abi_derfc(karg)
608 : ! see Honma eq. 2.10 for derivation of the following damping factor
609 24896 : dampfac = (1.0+3.0/(2.0*karg*karg))*exp(-karg*karg)+3.0*sqrt(pi)*derfc_karg/(4.0*karg**3)
610 99584 : do ii = 1, 3 ! loop over tensor elements
611 323648 : do jj = 1, 3 ! loop over tensor elements
612 224064 : derivs = -3.0*rhok(ii)*rhok(jj)/trace
613 224064 : if(ii == jj) derivs = derivs + 1.0 ! see Honma eq 4.8 re: sign
614 : ! accumulate real space tensor element,
615 : ! weighted by charge of neighbor and Ewald damping factor
616 298752 : efg_r(ii,jj,iatom) = efg_r(ii,jj,iatom) + qk*derivs*dampfac
617 : end do ! end loop over jj in efg(ii,jj,iatom)
618 : end do ! end loop over ii in efg(ii,jj,iatom)
619 : end if ! end if statement restricting to a sphere of radius rlkcut
620 : end if ! end if statement avoiding the self atom term
621 : end do ! end loop over i atoms in cell
622 : end do ! end loop over j atoms in cell
623 : end if ! end selection on outer shell of cells only
624 : end do ! end loop over sz cells
625 : end do ! end loop over sy cells
626 : end do ! end loop over sx cells
627 : end do ! end loop over shells
628 :
629 : !now combine the g-space and r-space parts, properly weighted (see Honma)
630 57 : do iatom = 1, natom
631 198 : do ii = 1, 3
632 611 : do jj = 1, 3
633 : efg(ii,jj,iatom) = four_pi*efg_g(ii,jj,iatom)/(three*ucvol)-&
634 564 : & four*xi0**3*efg_r(ii,jj,iatom)/(three*sqrt(pi))
635 : ! note extra factor of two: compare Honma eq. 4.6
636 : end do
637 : end do
638 : end do
639 :
640 : ! symmetrize tensor at each atomic site using point symmetry operations
641 57 : do iatom = 1, natom
642 57 : call matpointsym(iatom,efg(:,:,iatom),natom,nsym,rprimd,symrel,tnons,xred)
643 : end do
644 :
645 10 : ABI_FREE(efg_g)
646 10 : ABI_FREE(efg_r)
647 10 : ABI_FREE(xcart)
648 :
649 : !DEBUG
650 : !write(std_out,*)' make_efg_ion : exit '
651 : !stop
652 : !ENDDEBUG
653 :
654 10 : end subroutine make_efg_ion
655 : !!***
656 :
657 : !!****f* ABINIT/make_efg_el
658 : !! NAME
659 : !! make_efg_el
660 : !!
661 : !! FUNCTION
662 : !! compute the electric field gradient due to electron density
663 : !!
664 : !! INPUTS
665 : !! mpi_enreg=information about MPI parallelization
666 : !! natom, number of atoms in unit cell
667 : !! nfft,ngfft(18), number of FFT points and details of FFT
668 : !! nhat(nfft,nspden) compensation charge density
669 : !! nspden, number of spin components
670 : !! nsym=number of symmetries in space group
671 : !! rhor(nfft,nspden), valence electron density, here $\tilde{n} + \hat{n}$
672 : !! rprimd(3,3), conversion from crystal coordinates to cartesian coordinates
673 : !! symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
674 : !! tnons(3,nsym) = nonsymmorphic translations
675 : !! xred(3,natom), location of atoms in crystal coordinates.
676 : !!
677 : !! OUTPUT
678 : !! efg(3,3,natom), the 3x3 efg tensor at each atomic site due to rhor
679 : !!
680 : !! NOTES
681 : !! This routine computes the electric field gradient, specifically the components
682 : !! $\partial^2 V/\partial x_\alpha \partial x_\beta$ of the potential generated by the valence
683 : !! electrons, at each atomic site in the unit cell. Key references: Kresse and Joubert, ``From
684 : !! ultrasoft pseudopotentials to the projector augmented wave method'', Phys. Rev. B. 59, 1758--1775 (1999) [[cite:Kresse1999]],
685 : !! and Profeta, Mauri, and Pickard, ``Accurate first principles prediction of $^{17}$O NMR parameters in
686 : !! SiO$_2$: Assignment of the zeolite ferrierite spectrum'', J. Am. Chem. Soc. 125, 541--548 (2003) [[cite:Profeta2003]]. This
687 : !! routine computes the second derivatives of the potential generated by $\tilde{n}$ (see Kresse and Joubert
688 : !! for notation, Fourier-transforming the density, doing the sum in G space, and then transforming back at
689 : !! each atomic site. The final formula is
690 : !! \begin{displaymath}
691 : !! \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)
692 : !! \left(\frac{\tilde{n}(G)}{\pi G^2}\right)e^{2\pi i G\cdot R}
693 : !! \end{displaymath}
694 : !!
695 : !!
696 : !! SOURCE
697 :
698 9 : subroutine make_efg_el(efg,mpi_enreg,natom,nfft,ngfft,nhat,nspden,nsym,rhor,rprimd,symrel,tnons,xred)
699 :
700 : !Arguments ------------------------------------
701 : !scalars
702 : integer,intent(in) :: natom,nfft,nspden,nsym
703 : type(MPI_type),intent(in) :: mpi_enreg
704 : !arrays
705 : integer,intent(in) :: ngfft(18),symrel(3,3,nsym)
706 : real(dp),intent(in) :: nhat(nfft,nspden),rhor(nfft,nspden),rprimd(3,3),tnons(3,nsym),xred(3,natom)
707 : real(dp),intent(out) :: efg(3,3,natom)
708 :
709 : !Local variables-------------------------------
710 : !scalars
711 : integer :: cplex,fftdir,fofg_index,iatom,i1,i2,i2_local,i23,i3,id1,id2,id3
712 : integer :: ierr,ig,ig2,ig3,ii,ii1,ing,jj
713 : integer :: me_fft,n1,n2,n3,nproc_fft,tim_fourdp
714 : real(dp) :: cph,derivs,phase,sph,trace
715 : ! type(MPI_type) :: mpi_enreg_seq
716 : !arrays
717 : integer :: id(3)
718 9 : integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
719 9 : integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
720 : real(dp) :: gprimd(3,3),gqred(3),gvec(3),ratom(3)
721 9 : real(dp),allocatable :: fofg(:,:),fofr(:),gq(:,:),xcart(:,:)
722 :
723 : ! ************************************************************************
724 :
725 : !DEBUG
726 : !write(std_out,*)' make_efg_el : enter'
727 : !ENDDEBUG
728 :
729 27 : ABI_MALLOC(fofg,(2,nfft))
730 27 : ABI_MALLOC(fofr,(nfft))
731 27 : ABI_MALLOC(xcart,(3,natom))
732 :
733 542 : efg(:,:,:) = zero
734 9 : call xred2xcart(natom,rprimd,xcart,xred) ! get atomic locations in cartesian coords
735 9 : call matr3inv(rprimd,gprimd)
736 :
737 9 : tim_fourdp = 0 ! timing code, not using
738 9 : fftdir = -1 ! FT from R to G
739 9 : cplex = 1 ! fofr is real
740 : !here we are only interested in the valence pseudo charge density, which is rhor(:,1)-nhat(:,1)
741 : !regardless of the value of nspden. This may change in the future depending on
742 : !developments with noncollinear magnetization and so forth. Such a change will
743 : !require an additional loop over nspden.
744 : !Multiply by -1 to convert the electron particle density to the charge density
745 644537 : fofr(:) = -(rhor(:,1)-nhat(:,1))
746 :
747 : ! Get the distrib associated with this fft_grid See hartre.F90 for another example where
748 : ! this is done
749 9 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
750 9 : nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
751 9 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
752 :
753 9 : call fourdp(cplex,fofg,fofr,fftdir,mpi_enreg,nfft,1,ngfft,tim_fourdp) ! construct charge density in G space
754 :
755 : ! the following loops over G vectors has been copied from hartre.F90 in order to be compatible with
756 : ! possible FFT parallelism
757 :
758 : ! In order to speed the routine, precompute the components of g
759 : ! Also check if the booked space was large enough...
760 27 : ABI_MALLOC(gq,(3,max(n1,n2,n3)))
761 36 : do ii=1,3
762 27 : id(ii)=ngfft(ii)/2+2
763 1084 : do ing=1,ngfft(ii)
764 1048 : ig=ing-(ing/id(ii))*ngfft(ii)-1
765 1075 : gq(ii,ing)=ig
766 : end do
767 : end do
768 9 : id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
769 :
770 : ! Triple loop on each dimension
771 309 : do i3=1,n3
772 300 : ig3=i3-(i3/id3)*n3-1
773 300 : gqred(3) = gq(3,i3)
774 :
775 13853 : do i2=1,n2
776 13544 : ig2=i2-(i2/id2)*n2-1
777 13844 : if (fftn2_distrib(i2) == me_fft) then
778 :
779 13544 : gqred(2) = gq(2,i2)
780 13544 : i2_local = ffti2_local(i2)
781 13544 : i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1))
782 : ! Do the test that eliminates the Gamma point outside of the inner loop
783 13544 : ii1=1
784 13544 : if(i23==0 .and. ig2==0 .and. ig3==0) ii1=2
785 :
786 : ! Final inner loop on the first dimension (note the lower limit)
787 658063 : do i1=ii1,n1
788 : ! gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
789 644519 : gqred(1) = gq(1,i1)
790 8378747 : gvec(1:3) = MATMUL(gprimd,gqred)
791 644519 : fofg_index=i1+i23
792 2578076 : trace = dot_product(gvec,gvec)
793 2591620 : do ii = 1, 3 ! sum over components of efg tensor
794 8378747 : do jj = 1, 3 ! sum over components of efg tensor
795 5800671 : derivs = gvec(ii)*gvec(jj) ! This term is $G_\alpha G_\beta$
796 5800671 : if (ii == jj) derivs = derivs - trace/three
797 34391427 : do iatom = 1, natom ! sum over atoms in unit cell
798 106628796 : ratom(:) = xcart(:,iatom) ! extract location of atom iatom
799 106628796 : phase = two_pi*dot_product(gvec,ratom) ! argument of $e^{2\pi i G\cdot R}$
800 26657199 : cph = cos(phase)
801 26657199 : sph = sin(phase)
802 : efg(ii,jj,iatom) = efg(ii,jj,iatom) - &
803 32457870 : & four_pi*derivs*(fofg(1,fofg_index)*cph-fofg(2,fofg_index)*sph)/trace ! real part of efg tensor
804 : end do ! end loop over atoms in cell
805 : end do ! end loop over jj in V_ij
806 : end do ! end loop over ii in V_ij
807 : end do ! End loop on i1
808 : end if
809 : end do ! End loop on i2
810 : end do ! End loop on i3
811 :
812 9 : call xmpi_sum(efg,mpi_enreg%comm_fft,ierr)
813 :
814 : ! symmetrize tensor at each atomic site using point symmetry operations
815 50 : do iatom = 1, natom
816 50 : call matpointsym(iatom,efg(:,:,iatom),natom,nsym,rprimd,symrel,tnons,xred)
817 : end do
818 :
819 9 : ABI_FREE(fofg)
820 9 : ABI_FREE(fofr)
821 9 : ABI_FREE(xcart)
822 9 : ABI_FREE(gq)
823 :
824 : !DEBUG
825 : !write(std_out,*)' make_efg_el : exit '
826 : !stop
827 : !ENDDEBUG
828 :
829 9 : end subroutine make_efg_el
830 : !!***
831 :
832 : end module m_nucprop
833 : !!***
|
