Branch data 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, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
719 : 9 : integer, 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 : : !!***
|
