Branch data Line data Source code
1 : : !!****m* ABINIT/m_mkcore
2 : : !! NAME
3 : : !! m_mkcore
4 : : !!
5 : : !! FUNCTION
6 : : !! Routines related to non-linear core correction.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, TRangel, MT)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_mkcore
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_xmpi
27 : : use m_errors
28 : : use m_linalg_interfaces
29 : :
30 : : use defs_abitypes, only : mpi_type
31 : : use m_geometry, only : strconv
32 : : use m_time, only : timab
33 : : use m_mpinfo, only : ptabs_fourdp
34 : : use m_sort, only : sort_dp
35 : : use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free
36 : : use m_pawtab, only : pawtab_type
37 : : use m_paw_numeric, only : paw_splint
38 : :
39 : : implicit none
40 : :
41 : : private
42 : : !!***
43 : :
44 : : public :: mkcore
45 : : public :: mkcore_alt
46 : : public :: dfpt_mkcore ! Derivative of the core electron density with respect to one specific atom displacement
47 : : !!***
48 : :
49 : : contains
50 : : !!***
51 : :
52 : : !!****f* ABINIT/mkcore
53 : : !! NAME
54 : : !! mkcore
55 : : !!
56 : : !! FUNCTION
57 : : !! Optionally compute:
58 : : !! (1) pseudo core electron density throughout unit cell
59 : : !! (2) pseudo-core contribution to forces
60 : : !! (3) pseudo-core contribution to stress tensor
61 : : !! (4) pseudo-core contrib. to frozen-wf part the dynamical matrix (part 2)
62 : : !!
63 : : !! INPUTS
64 : : !! natom=number of atoms in cell.
65 : : !! nfft=(effective) number of FFT grid points (for this processor)
66 : : !! nspden=number of spin-density components
67 : : !! ntypat=number of types of atoms in cell.
68 : : !! n1,n2,n3=fft grid dimensions.
69 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
70 : : !! option: 1 for computing xccc3d (core charge density),
71 : : !! 2 for computing core charge contribution to $d(E_{xc})/d(tau)$,
72 : : !! 3 for computing core charge contribution to stress tensor corstr,
73 : : !! 4 for contribution to frozen-wavefunction part of dynamical matrix
74 : : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
75 : : !! typat(natom)=integer type for each atom in cell
76 : : !! ucvol=unit cell volume (bohr**3).
77 : : !! vxc(nfft,nspden)=exchange-correlation potential (hartree) in real
78 : : !! space--only used when option=2,3, or 4, else ignored
79 : : !! xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
80 : : !! xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives,
81 : : !! for each type of atom, from psp
82 : : !! xred(3,natom)=reduced coordinates for atoms in unit cell
83 : : !!
84 : : !! OUTPUT
85 : : !! corstr(6)=core charge contribution to stress tensor, only if option=3
86 : : !! dyfrx2(3,3,natom)=non-linear xc core correction part of the
87 : : !! frozen-wavefunction part of the dynamical matrix, only for option=4
88 : : !! grxc(3,natom)=d(Exc)/d(xred), hartree (only computed when option=2, else
89 : : !! ignored)
90 : : !!
91 : : !! SIDE EFFECTS
92 : : !! xccc3d(n1*n2*n3)=3D core electron density for XC core correction, bohr^-3
93 : : !! (computed and returned when option=1, needed as input when option=3)
94 : : !!
95 : : !! NOTES
96 : : !! Note that this routine is tightly connected to the dfpt_mkcore.f routine
97 : : !!
98 : : !! SOURCE
99 : :
100 : 9539 : subroutine mkcore(corstr,dyfrx2,grxc,mpi_enreg,natom,nfft,nspden,ntypat,n1,n1xccc,&
101 : 9539 : & n2,n3,option,rprimd,typat,ucvol,vxc,xcccrc,xccc1d,xccc3d,xred)
102 : :
103 : : !Arguments ------------------------------------
104 : : !scalars
105 : : integer,intent(in) :: n1,n1xccc,n2,n3,natom,nfft,nspden,ntypat,option
106 : : real(dp),intent(in) :: ucvol
107 : : type(mpi_type),intent(in) :: mpi_enreg
108 : : !arrays
109 : : integer,intent(in) :: typat(natom)
110 : : real(dp),intent(in) :: rprimd(3,3),vxc(nfft,nspden),xccc1d(n1xccc,6,ntypat)
111 : : real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom)
112 : : real(dp),intent(inout) :: xccc3d(nfft)
113 : : real(dp),intent(out) :: corstr(6),dyfrx2(3,3,natom)
114 : : real(dp),intent(inout) :: grxc(3,natom)
115 : :
116 : : !Local variables-------------------------------
117 : : !scalars
118 : : integer :: i1,i2,i3,iatom,ier,ifft,ishift,ishift1,ishift2
119 : : integer :: ishift3,itypat,ixp,jj,me_fft,mrange,mu,nfftot,nu
120 : : real(dp) :: dd,delta,delta2div6,deltam1,diff,difmag,difmag2
121 : : real(dp) :: difmag2_fact,difmag2_part,fact,func,grxc1,grxc2,grxc3,range,range2
122 : : real(dp) :: rangem1,rdiff1,rdiff2,rdiff3,strdia,t1,t2,t3,term,term1,term2
123 : : character(len=500) :: message
124 : : !arrays
125 : : integer :: igrid(3),irange(3),ngfft(3)
126 : 9539 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
127 : 9539 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
128 : 9539 : integer,allocatable :: ii(:,:)
129 : : real(dp) :: yy,aa,bb,cc
130 : : real(dp) :: corfra(3,3),lencp(3),rmet(3,3),scale(3),tau(3),tsec(2),tt(3)
131 : 9539 : real(dp),allocatable :: rrdiff(:,:),work(:,:,:)
132 : : !************************************************************************
133 : :
134 : 9539 : call timab(12,1,tsec)
135 : :
136 : : !Make sure option is acceptable
137 [ - + ]: 9539 : if (option<0.or.option>4) then
138 : : write(message, '(a,i12,a,a,a)' )&
139 : 0 : 'option=',option,' is not allowed.',ch10,&
140 : 0 : 'Must be 1, 2, 3 or 4.'
141 : 0 : ABI_BUG(message)
142 : : end if
143 : :
144 : : !Zero out only the appropriate array according to option:
145 : : !others are dummies with no storage
146 : :
147 [ + + ]: 9539 : if (option==1) then
148 : : ! Zero out array to permit accumulation over atom types below:
149 [ + + ]: 49544811 : xccc3d(:)=zero
150 [ + + ]: 6164 : else if (option==2) then
151 : : ! Zero out gradient of Exc array
152 [ + + + + ]: 63309 : grxc(:,:)=zero
153 [ + + ]: 2063 : else if (option==3) then
154 : : ! Zero out locally defined stress array
155 : 9539 : corfra(:,:)=zero
156 : 9539 : strdia=zero
157 [ + - ]: 246 : else if (option==4) then
158 : : ! Zero out fr-wf part of the dynamical matrix
159 [ + + + + : 7175 : dyfrx2(:,:,:)=zero
+ + ]
160 : : else
161 : 0 : ABI_BUG(" Can't be here! (bad option)")
162 : : end if
163 : :
164 : : !Compute lengths of cross products for pairs of primitive
165 : : !translation vectors (used in setting index search range below)
166 : : lencp(1)=cross_mkcore(rprimd(1,2),rprimd(2,2),rprimd(3,2),&
167 : 9539 : & rprimd(1,3),rprimd(2,3),rprimd(3,3))
168 : : lencp(2)=cross_mkcore(rprimd(1,3),rprimd(2,3),rprimd(3,3),&
169 : 9539 : & rprimd(1,1),rprimd(2,1),rprimd(3,1))
170 : : lencp(3)=cross_mkcore(rprimd(1,1),rprimd(2,1),rprimd(3,1),&
171 : 9539 : & rprimd(1,2),rprimd(2,2),rprimd(3,2))
172 : :
173 : : !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3
174 : : !(recall ucvol=R1.(R2xR3))
175 [ + + ]: 38156 : scale(:)=ucvol/lencp(:)
176 : :
177 : : !Compute metric tensor in real space rmet
178 [ + + ]: 38156 : do nu=1,3
179 [ + + ]: 124007 : rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+rprimd(3,:)*rprimd(3,nu)
180 : : end do
181 : :
182 : 9539 : ngfft(1)=n1
183 : 9539 : ngfft(2)=n2
184 : 9539 : ngfft(3)=n3
185 : 9539 : nfftot=n1*n2*n3
186 : 9539 : me_fft = mpi_enreg%me_fft
187 : :
188 : : ! Get the distrib associated with this fft_grid
189 : 9539 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
190 : :
191 : 9539 : delta=one/(n1xccc-1)
192 : 9539 : deltam1=n1xccc-1
193 : 9539 : delta2div6=delta**2/6.0d0
194 : :
195 [ + + ]: 9539 : if (option>=2) then
196 [ - + + - : 30820 : ABI_MALLOC(work,(n1,n2,n3))
- + + - -
+ + - + -
+ - - + -
+ ]
197 : : ! For spin-polarization, replace vxc by (1/2)*(vxc(up)+vxc(down))
198 : : ! For non-collinear magnetism, replace vxc by (1/2)*(vxc^{11}+vxc^{22})
199 [ + + ]: 6164 : if (nspden>=2) then
200 : : ifft=1
201 [ + + ]: 25140 : do i3=1,n3
202 [ + - ]: 25140 : if(me_fft==fftn3_distrib(i3)) then
203 [ + + ]: 582472 : do i2=1,n2
204 [ + + ]: 14456818 : do i1=1,n1
205 : 13874346 : work(i1,i2,i3)=half*(vxc(ifft,1)+vxc(ifft,2))
206 : 14432725 : ifft=ifft+1
207 : : end do
208 : : end do
209 : : end if
210 : : end do
211 : : else
212 : : ifft=1
213 [ + + ]: 130196 : do i3=1,n3
214 [ + + ]: 130196 : if(me_fft==fftn3_distrib(i3)) then
215 [ + + ]: 3442173 : do i2=1,n2
216 [ + + ]: 126854617 : do i1=1,n1
217 : 123412444 : work(i1,i2,i3)=vxc(ifft,1)
218 : 126730338 : ifft=ifft+1
219 : : end do
220 : : end do
221 : : end if
222 : : end do
223 : : ! call DCOPY(nfft,vxc,1,work,1)
224 : : end if
225 : : end if
226 : :
227 : : !Loop over atoms in unit cell
228 [ + + ]: 43716 : do iatom=1,natom
229 : :
230 [ + + ]: 34177 : if(option==2)then
231 : 14802 : grxc1=zero
232 : 14802 : grxc2=zero
233 : 14802 : grxc3=zero
234 : : end if
235 : :
236 : : ! Set search range (density cuts off perfectly beyond range)
237 : 34177 : itypat=typat(iatom)
238 : 34177 : range=xcccrc(itypat)
239 : :
240 : : ! Skip loop if this atom has no core charge
241 [ + + ]: 34177 : if (abs(range)<1.d-16) cycle
242 : :
243 : 31831 : range2=range**2
244 : 31831 : rangem1=one/range
245 : :
246 : : ! Consider each component in turn : compute range
247 [ + + ]: 127324 : do mu=1,3
248 : :
249 : : ! Convert reduced coord of given atom to [0,1)
250 : 95493 : tau(mu)=mod(xred(mu,iatom)+one-aint(xred(mu,iatom)),one)
251 : :
252 : : ! Use tau to find nearest grid point along R(mu)
253 : : ! (igrid=0 is the origin; shift by 1 to agree with usual index)
254 : 95493 : igrid(mu)=nint(tau(mu)*dble(ngfft(mu)))
255 : :
256 : : ! Use range to compute an index range along R(mu)
257 : : ! (add 1 to make sure it covers full range)
258 : 127324 : irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu)))
259 : :
260 : : end do
261 : :
262 : : ! Allocate arrays that depends on the range
263 [ + + ]: 127324 : mrange=maxval(irange(1:3))
264 [ + - - + : 127324 : ABI_MALLOC(ii,(2*mrange+1,3))
+ - + - +
- + - - +
- + ]
265 [ - + + - : 95493 : ABI_MALLOC(rrdiff,(2*mrange+1,3))
+ - + - -
+ - + ]
266 : :
267 : : ! Set up counters that explore the relevant range
268 : : ! of points around the atom
269 [ + + ]: 127324 : do mu=1,3
270 : 95493 : ishift=0
271 [ + + ]: 2927119 : do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
272 : 2799795 : ishift=ishift+1
273 : 2799795 : ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu))
274 : 2895288 : rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu)
275 : : end do
276 : : end do
277 : :
278 : : ! Conduct triple loop over restricted range of grid points for iatom
279 [ + + ]: 955404 : do ishift3=1,1+2*irange(3)
280 : : ! map back to [1,ngfft(3)] for usual fortran index in unit cell
281 : 923573 : i3=ii(ishift3,3)
282 [ + + ]: 923573 : if(fftn3_distrib(i3)/=mpi_enreg%me_fft) cycle
283 : : ! find vector from atom location to grid point (reduced)
284 : 915083 : rdiff3=rrdiff(ishift3,3)
285 : :
286 [ + + ]: 32957825 : do ishift2=1,1+2*irange(2)
287 : 32010911 : i2=ii(ishift2,2)
288 : 32010911 : rdiff2=rrdiff(ishift2,2)
289 : : ! Prepare the computation of difmag2
290 : : difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2&
291 : 32010911 : & +2.0d0*rmet(3,2)*rdiff3*rdiff2
292 : 32010911 : difmag2_fact=2.0d0*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2)
293 : :
294 [ + + ]: 1415101127 : do ishift1=1,1+2*irange(1)
295 : 1382166643 : rdiff1=rrdiff(ishift1,1)
296 : :
297 : : ! Compute (rgrid-tau-Rprim)**2
298 : 1382166643 : difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1)
299 : :
300 : : ! Only accept contribution inside defined range
301 [ + + ]: 1414177554 : if (difmag2<range2-tol12) then
302 : :
303 : : ! Prepare computation of core charge function and derivative,
304 : : ! using splines
305 : 502758795 : i1=ii(ishift1,1)
306 : 502758795 : difmag=sqrt(difmag2)
307 : 502758795 : yy=difmag*rangem1
308 : :
309 : : ! Compute index of yy over 1 to n1xccc scale
310 : 502758795 : jj=1+int(yy*(n1xccc-1))
311 : 502758795 : diff=yy-(jj-1)*delta
312 : :
313 : : ! Will evaluate spline fit (p. 86 Numerical Recipes, Press et al;
314 : : ! NOTE error in book for sign of "aa" term in derivative;
315 : : ! also see splfit routine).
316 : 502758795 : bb = diff*deltam1
317 : 502758795 : aa = one-bb
318 : 502758795 : cc = aa*(aa**2-one)*delta2div6
319 : 502758795 : dd = bb*(bb**2-one)*delta2div6
320 : :
321 : :
322 : : ! Test first for option 2, the most frequently used
323 [ + + ]: 502758795 : if (option==2) then
324 : :
325 : : ! Accumulate contributions to Exc gradients
326 : :
327 [ + + ]: 169150370 : if (difmag>1.0d-10) then
328 : :
329 : : ! Evaluate spline fit of 1st der of core charge density
330 : : ! from xccc1d(:,2,:) and (:,4,:)
331 : : func=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +&
332 : 169147441 : & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat)
333 : 169147441 : term=work(i1,i2,i3)*func/difmag
334 : 169147441 : grxc1=grxc1+rdiff1*term
335 : 169147441 : grxc2=grxc2+rdiff2*term
336 : 169147441 : grxc3=grxc3+rdiff3*term
337 : : end if
338 : :
339 [ + + ]: 333608425 : else if (option==1) then
340 : :
341 : : ! Evaluate spline fit of core charge density
342 : : ! from xccc1d(:,1,:) and (:,3,:)
343 : : func=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +&
344 : 218189533 : & cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat)
345 : :
346 : : ! Accumulate contributions to core electron density
347 : : ! throughout unit cell
348 : 218189533 : ifft=i1+n1*(i2-1+n2*(ffti3_local(i3)-1))
349 : 218189533 : xccc3d(ifft)=xccc3d(ifft)+func
350 : :
351 [ + + ]: 115418892 : else if (option==3) then
352 : :
353 : : ! Accumulate contributions to stress tensor
354 : : ! in reduced coordinates
355 : :
356 [ + + ]: 93681706 : if (difmag>1.0d-10) then
357 : :
358 : : ! Evaluate spline fit of 1st der of core charge density
359 : : ! from xccc1d(:,2,:) and (:,4,:)
360 : : func=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +&
361 : 93680102 : & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat)
362 : 93680102 : term=work(i1,i2,i3)*func*rangem1/difmag/dble(n1*n2*n3)
363 : : ! Write out the 6 symmetric components
364 : 93680102 : corfra(1,1)=corfra(1,1)+term*rdiff1**2
365 : 93680102 : corfra(2,2)=corfra(2,2)+term*rdiff2**2
366 : 93680102 : corfra(3,3)=corfra(3,3)+term*rdiff3**2
367 : 93680102 : corfra(3,2)=corfra(3,2)+term*rdiff3*rdiff2
368 : 93680102 : corfra(3,1)=corfra(3,1)+term*rdiff3*rdiff1
369 : 93680102 : corfra(2,1)=corfra(2,1)+term*rdiff2*rdiff1
370 : : ! (the above still needs to be transformed to cartesian coords)
371 : :
372 : : end if
373 : :
374 : : ! Also compute a diagonal term
375 : : ! Evaluate spline fit of core charge density
376 : : ! from xccc1d(:,1,:) and (:,3,:)
377 : : func=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +&
378 : 93681706 : & cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat)
379 : 93681706 : strdia=strdia+work(i1,i2,i3)*func
380 : :
381 [ + - ]: 21737186 : else if (option==4) then
382 : :
383 : : ! Compute frozen-wf contribution to Dynamical matrix
384 : :
385 : 21737186 : tt(1)=rmet(1,1)*rdiff1+rmet(1,2)*rdiff2+rmet(1,3)*rdiff3
386 : 21737186 : tt(2)=rmet(2,1)*rdiff1+rmet(2,2)*rdiff2+rmet(2,3)*rdiff3
387 : 21737186 : tt(3)=rmet(3,1)*rdiff1+rmet(3,2)*rdiff2+rmet(3,3)*rdiff3
388 : :
389 [ + + ]: 21737186 : if (difmag>1.d-10) then
390 : :
391 : : ! Accumulate contributions to dynamical matrix
392 : 21736924 : term=(ucvol/dble(nfftot))*work(i1,i2,i3)*rangem1/difmag
393 : : ! Evaluate spline fit of 1st der of core charge density
394 : : ! from xccc1d(:,2,:) and (:,4,:)
395 : : func=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +&
396 : 21736924 : & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat)
397 : 21736924 : term1=term*func
398 : : ! Evaluate spline fit of 2nd der of core charge density
399 : : ! from xccc1d(:,3,:) and (:,5,:)
400 : : func=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +&
401 : 21736924 : & cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat)
402 : 21736924 : term2=term*func*rangem1/difmag
403 [ + + ]: 86947696 : do mu=1,3
404 [ + + ]: 282580012 : do nu=1,3
405 : : dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)&
406 : : & +(term2-term1/difmag**2)*tt(mu)*tt(nu)&
407 : 260843088 : & +term1*rmet(mu,nu)
408 : : end do
409 : : end do
410 : :
411 : : else
412 : :
413 : : ! There is a contribution from difmag=zero !
414 : : ! Evaluate spline fit of 2nd der of core charge density
415 : : ! from xccc1d(:,3,:) and (:,5,:)
416 : : func=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +&
417 : 262 : & cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat)
418 : 262 : term=(ucvol/dble(nfftot))*work(i1,i2,i3)*func*rangem1**2
419 [ + + ]: 1048 : do mu=1,3
420 [ + + ]: 3406 : do nu=1,3
421 : 3144 : dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)+term*rmet(mu,nu)
422 : : end do
423 : : end do
424 : :
425 : : ! End of condition not to be precisely on the point (difmag=zero)
426 : : end if
427 : :
428 : : ! If option is not 1, 2, 3, or 4.
429 : : else
430 : 0 : ABI_BUG("Can't be here in mkcore")
431 : : ! End of choice of option
432 : : end if
433 : :
434 : : end if ! End of condition on the range
435 : : end do ! End loop on ishift1
436 : : end do ! End loop on ishift2
437 : : end do ! End loop on ishift3
438 : :
439 [ - + ]: 31831 : ABI_FREE(ii)
440 [ - + ]: 31831 : ABI_FREE(rrdiff)
441 : :
442 [ + + ]: 41370 : if(option==2)then
443 : 13391 : fact=-(ucvol/dble(nfftot))/range
444 : 13391 : grxc(1,iatom)=grxc1*fact
445 : 13391 : grxc(2,iatom)=grxc2*fact
446 : 13391 : grxc(3,iatom)=grxc3*fact
447 : : end if
448 : :
449 : : end do ! End big loop on atoms
450 : :
451 [ + + ]: 9539 : if (option==2) then
452 : :
453 : : ! Apply rmet as needed to get reduced coordinate gradients
454 [ + + ]: 18903 : do iatom=1,natom
455 : 14802 : t1=grxc(1,iatom)
456 : 14802 : t2=grxc(2,iatom)
457 : 14802 : t3=grxc(3,iatom)
458 [ + + ]: 63309 : grxc(:,iatom)=rmet(:,1)*t1+rmet(:,2)*t2+rmet(:,3)*t3
459 : :
460 : : end do
461 : : end if
462 : :
463 [ + + ]: 9539 : if (option==3) then
464 : :
465 : : ! Transform stress tensor from full storage mode to symmetric storage mode
466 : 1817 : corstr(1)=corfra(1,1)
467 : 1817 : corstr(2)=corfra(2,2)
468 : 1817 : corstr(3)=corfra(3,3)
469 : 1817 : corstr(4)=corfra(3,2)
470 : 1817 : corstr(5)=corfra(3,1)
471 : 1817 : corstr(6)=corfra(2,1)
472 : :
473 : : ! Transform stress tensor from reduced coordinates to cartesian coordinates
474 : 1817 : call strconv(corstr,rprimd,corstr)
475 : :
476 : : ! Compute diagonal contribution to stress tensor (need input xccc3d)
477 : : ! strdia = (1/N) Sum(r) [mu_xc_avg(r) * rho_core(r)]
478 : 1817 : ifft=0 ; strdia=zero
479 [ + + ]: 41529 : do i3=1,n3
480 [ + + ]: 41529 : if(me_fft==fftn3_distrib(i3)) then
481 [ + + ]: 949478 : do i2=1,n2
482 [ + + ]: 27660881 : do i1=1,n1
483 : 26711403 : ifft=ifft+1
484 : 27621569 : strdia=strdia+work(i1,i2,i3)*xccc3d(ifft)
485 : : end do
486 : : end do
487 : : end if
488 : : end do
489 : 1817 : strdia=strdia/dble(nfftot)
490 : : ! strdia=DDOT(nfft,work,1,xccc3d,1)/dble(nfftot)
491 : :
492 : : ! Add diagonal term to stress tensor
493 : 1817 : corstr(1)=corstr(1)+strdia
494 : 1817 : corstr(2)=corstr(2)+strdia
495 : 1817 : corstr(3)=corstr(3)+strdia
496 : : end if
497 : :
498 [ + + ]: 9539 : if(option>=2) then
499 [ - + ]: 6164 : ABI_FREE(work)
500 : : end if
501 : :
502 [ + + ]: 9539 : if(mpi_enreg%nproc_fft > 1) then
503 : 72 : call timab(539,1,tsec)
504 [ + + ]: 72 : if(option==2) then
505 : 24 : call xmpi_sum(grxc,mpi_enreg%comm_fft,ier)
506 : : end if
507 [ + + ]: 72 : if(option==3) then
508 : 24 : call xmpi_sum(corstr,mpi_enreg%comm_fft,ier)
509 : : end if
510 [ - + ]: 72 : if(option==4) then
511 : 0 : call xmpi_sum(dyfrx2,mpi_enreg%comm_fft,ier)
512 : : end if
513 : 72 : call timab(539,2,tsec)
514 : : end if
515 : :
516 [ - + ]: 19078 : call timab(12,2,tsec)
517 : :
518 : : contains
519 : : !!***
520 : :
521 : : !!****f* ABINIT/cross_mkcore
522 : : !! NAME
523 : : !! cross_mkcore
524 : : !!
525 : : !! FUNCTION
526 : : !! Define magnitude of cross product of two vectors
527 : : !!
528 : : !! SOURCE
529 : :
530 : 28617 : function cross_mkcore(xx,yy,zz,aa,bb,cc)
531 : :
532 : : real(dp) :: cross_mkcore
533 : : real(dp),intent(in) :: xx,yy,zz,aa,bb,cc
534 : : ! *************************************************************************
535 : 28617 : cross_mkcore=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2)
536 : 28617 : end function cross_mkcore
537 : :
538 : : end subroutine mkcore
539 : : !!***
540 : :
541 : : !--------------------------------------------------------------------------------------
542 : :
543 : : !!****f* ABINIT/mkcore_alt
544 : : !! NAME
545 : : !! mkcore_alt
546 : : !!
547 : : !! FUNCTION
548 : : !! Optionally compute:
549 : : !! (1) pseudo core electron density throughout unit cell
550 : : !! (2) pseudo-core contribution to forces
551 : : !! (3) pseudo-core contribution to stress tensor
552 : : !! (4) pseudo-core contrib. to frozen-wf part the dynamical matrix (part 2)
553 : : !! This routine is an alternative to mkcore, to be used for PAW and/or WVL.
554 : : !!
555 : : !! INPUTS
556 : : !! atindx1(natom)=index table for atoms, inverse of atindx
557 : : !! icoulomb= periodic treatment of Hartree potential: 0=periodic, 1=free BC, 2=surface BC
558 : : !! mpi_enreg=information about MPI parallelization
559 : : !! natom=number of atoms in cell.
560 : : !! nfft=(effective) number of FFT grid points (for this processor)
561 : : !! nspden=number of spin-density components
562 : : !! ntypat=number of types of atoms in cell
563 : : !! n1,n2,n3=fft grid dimensions.
564 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
565 : : !! option: 1 for computing core charge density
566 : : !! 2 for computing core charge contribution to forces
567 : : !! 3 for computing core charge contribution to stress tensor
568 : : !! 4 for computing contribution to frozen-wf part of dynamical matrix
569 : : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
570 : : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
571 : : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
572 : : !! ucvol=unit cell volume (bohr**3)
573 : : !! usepaw=flag for PAW method
574 : : !! vxc(nfft,nspden)=exchange-correlation potential (hartree) in real space
575 : : !! xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
576 : : !! xccc1d(n1xccc,6,ntypat)=1D core charge function and 5 derivatives for each atom type
577 : : !! xred(3,natom)=reduced coordinates for atoms in unit cell
578 : : !! [usekden]= --optional-- if TRUE, output the kinetic energy density instead of the density
579 : : !!
580 : : !! OUTPUT
581 : : !! === if option==1 ===
582 : : !! xccc3d(n1*n2*n3)=3D core electron density for XC core correction (bohr^-3)
583 : : !! === if option==2 ===
584 : : !! grxc(3,natom)=core charge contribution to forces
585 : : !! === if option==3 ===
586 : : !! corstr(6)=core charge contribution to stress tensor
587 : : !! === if option==4 ===
588 : : !! dyfrx2(3,3,natom)=non-linear xc core correction part of the
589 : : !! frozen-wavefunction part of the dynamical matrix
590 : : !!
591 : : !! SIDE EFFECTS
592 : : !! xccc3d(n1*n2*n3)=3D core electron (event. kinetic energy) density for XC core correction (bohr^-3)
593 : : !! (computed and returned when option=1, needed as input when option=3)
594 : : !!
595 : : !! NOTES
596 : : !! Based on mkcore.F90
597 : : !!
598 : : !! SOURCE
599 : :
600 : 92 : subroutine mkcore_alt(atindx1,corstr,dyfrx2,grxc,icoulomb,mpi_enreg,natom,nfft,nspden,&
601 : 92 : & nattyp,ntypat,n1,n1xccc,n2,n3,option,rprimd,ucvol,vxc,xcccrc,xccc1d,&
602 [ - + - + ]: 46 : & xccc3d,xred,pawrad,pawtab,usepaw,&
603 : : & usekden) ! optional argument
604 : :
605 : : !Arguments ------------------------------------
606 : : !scalars
607 : : integer,intent(in) :: icoulomb,n1,n1xccc,n2,n3,natom,nfft,nspden,ntypat,option,usepaw
608 : : logical,intent(in),optional :: usekden
609 : : real(dp),intent(in) :: ucvol
610 : : type(mpi_type),intent(in) :: mpi_enreg
611 : : type(pawrad_type),intent(in) :: pawrad(:)
612 : : type(pawtab_type),target,intent(in) :: pawtab(:)
613 : : !arrays
614 : : integer,intent(in) :: atindx1(natom),nattyp(ntypat)
615 : : real(dp),intent(in) :: rprimd(3,3),xccc1d(n1xccc,6,ntypat)
616 : : real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom)
617 : : real(dp),intent(in),target :: vxc(nfft,nspden)
618 : : real(dp),intent(out) :: corstr(6),grxc(3,natom),dyfrx2(3,3,natom)
619 : : real(dp),intent(inout) :: xccc3d(nfft)
620 : :
621 : : !Local variables-------------------------------
622 : : !scalars
623 : : integer :: i1,i2,i3,iat,iatm,iatom,ier,ipts
624 : : integer :: ishift,ishift1,ishift2,ishift3
625 : : integer :: itypat,ixp,jj,jpts,me_fft,mrange,msz,mu
626 : : integer :: nfftot,npts,npts12,nu
627 : : logical :: letsgo,usekden_
628 : : real(dp) :: aa,bb,cc,dd,delta,delta2div6,deltam1
629 : : real(dp) :: diff,difmag,fact,range,range2
630 : : real(dp) :: rangem1,rdiff1,rdiff2,rdiff3
631 : : real(dp) :: rnorm2,rnorm2_fact,rnorm2_part
632 : : real(dp) :: strdia,t1,t2,t3,term,term1,term2,yy
633 : : character(len=1) :: geocode
634 : : character(len=500) :: message
635 : 46 : type(pawrad_type) :: core_mesh
636 : : !arrays
637 : : integer :: igrid(3),irange(3),ishiftmax(3),ngfft(3)
638 : 46 : integer,allocatable :: ii(:,:),iindex(:),indx1(:),indx2(:)
639 : 46 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
640 : 46 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
641 : : logical :: per(3)
642 : : real(dp) :: corfra(3,3),corgr(3),lencp(3),rmet(3,3)
643 : : real(dp) :: scale(3),tau(3),tsec(2),tt(3)
644 : 46 : real(dp),allocatable :: dtcore(:),d2tcore(:),rnorm(:)
645 : 46 : real(dp),allocatable :: rrdiff(:,:),tcore(:)
646 : 46 : real(dp),allocatable,target :: tcoretau(:,:) ! only needed in PAW case.
647 : 46 : real(dp), contiguous, pointer :: corespl(:,:),vxc_eff(:)
648 : : !************************************************************************
649 : :
650 : 46 : call timab(12,1,tsec)
651 : :
652 : : !Make sure options are acceptable
653 [ - + ]: 46 : if (option<0.or.option>4) then
654 : : write(message, '(a,i12,a,a,a)' )&
655 : 0 : 'option=',option,' is not allowed.',ch10,&
656 : 0 : 'Must be 1, 2, 3 or 4.'
657 : 0 : ABI_BUG(message)
658 : : end if
659 : :
660 [ + - ]: 46 : usekden_=.false.;if (present(usekden)) usekden_=usekden
661 [ + - ]: 46 : if (usekden_) then
662 : 46 : message='usekden=1 mkcore_alt not yet in production. You have been warned! May not work with PAW or NC'
663 : 46 : ABI_WARNING(message)
664 : : end if
665 : :
666 : :
667 : : !Zero out only the appropriate array according to option:
668 [ + + ]: 46 : if (option==1) then
669 [ + + ]: 2003084 : xccc3d(:)=zero
670 [ + + ]: 6 : else if (option==2) then
671 [ + + + + ]: 23 : grxc(:,:)=zero
672 [ - + ]: 3 : else if (option==3) then
673 : 46 : corfra(:,:)=zero
674 : 46 : strdia=zero
675 [ # # ]: 0 : else if (option==4) then
676 [ # # # # : 0 : dyfrx2(:,:,:)=zero
# # ]
677 : : end if
678 : :
679 : : !Conditions for periodicity in the three directions
680 : 46 : geocode='P'
681 [ - + ]: 46 : if (icoulomb==1) geocode='F'
682 [ - + ]: 46 : if (icoulomb==2) geocode='S'
683 : 46 : per(1)=(geocode /= 'F')
684 : 46 : per(2)=(geocode == 'P')
685 : 46 : per(3)=(geocode /= 'F')
686 : :
687 : : !Compute lengths of cross products for pairs of primitive
688 : : !translation vectors (used in setting index search range below)
689 : : lencp(1)=cross_mkcore_alt(rprimd(1,2),rprimd(2,2),rprimd(3,2),&
690 : 46 : & rprimd(1,3),rprimd(2,3),rprimd(3,3))
691 : : lencp(2)=cross_mkcore_alt(rprimd(1,3),rprimd(2,3),rprimd(3,3),&
692 : 46 : & rprimd(1,1),rprimd(2,1),rprimd(3,1))
693 : : lencp(3)=cross_mkcore_alt(rprimd(1,1),rprimd(2,1),rprimd(3,1),&
694 : 46 : & rprimd(1,2),rprimd(2,2),rprimd(3,2))
695 : :
696 : : !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3
697 : : !(recall ucvol=R1.(R2xR3))
698 [ + + ]: 184 : scale(:)=ucvol/lencp(:)
699 : :
700 : : !Compute metric tensor in real space rmet
701 [ + + ]: 184 : do nu=1,3
702 [ + + ]: 598 : rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+rprimd(3,:)*rprimd(3,nu)
703 : : end do
704 : :
705 : : !Get the distrib associated with this fft_grid
706 : 46 : ngfft(1)=n1;ngfft(2)=n2;ngfft(3)=n3
707 : 46 : nfftot=n1*n2*n3 ; me_fft=mpi_enreg%me_fft
708 : 46 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
709 : :
710 : 46 : delta=one/(n1xccc-1)
711 : 46 : deltam1=n1xccc-1
712 : 46 : delta2div6=delta**2/6.0_dp
713 : :
714 [ + + ]: 46 : if (option>=2) then
715 : : ! For spin-polarization, replace vxc by (1/2)*(vxc(up)+vxc(down))
716 : : ! For non-collinear magnetism, replace vxc by (1/2)*(vxc^{11}+vxc^{22})
717 [ - + ]: 6 : if (nspden>=2) then
718 [ # # # # : 0 : ABI_MALLOC(vxc_eff,(nfft))
# # # # #
# # # ]
719 [ # # ]: 0 : do jj=1,nfft
720 : 0 : vxc_eff(jj)=half*(vxc(jj,1)+vxc(jj,2))
721 : : end do
722 : : else
723 : 6 : vxc_eff => vxc(1:nfft,1)
724 : : end if
725 : : end if
726 : :
727 : : !Loop over atom types
728 : 46 : iatm=0
729 [ + + ]: 100 : do itypat=1,ntypat
730 : :
731 : : ! Set search range (density cuts off perfectly beyond range)
732 [ - + ]: 54 : range=xcccrc(itypat);if (usepaw==1) range=pawtab(itypat)%rcore
733 : :
734 : : ! Skip loop if this type has no core charge
735 [ + + ]: 54 : if (abs(range)<1.d-16) cycle
736 : :
737 : 15 : range2=range**2 ; rangem1=one/range
738 : :
739 : :
740 : : ! PAW: select core density type and create mesh
741 [ - + ]: 15 : if (usepaw==1) then
742 [ # # ]: 0 : if (usekden_) then
743 : 0 : msz=pawtab(itypat)%coretau_mesh_size
744 [ # # # # : 0 : ABI_MALLOC(tcoretau,(msz,1))
# # # # #
# # # #
# ]
745 [ # # ]: 0 : tcoretau(:,1)=pawtab(itypat)%coretau(:)
746 : 0 : corespl => tcoretau
747 : : else
748 : 0 : msz=pawtab(itypat)%core_mesh_size
749 : 0 : corespl => pawtab(itypat)%tcoredens
750 : : end if
751 : : call pawrad_init(core_mesh,mesh_size=msz,&
752 : : & mesh_type=pawrad(itypat)%mesh_type,&
753 : 0 : & rstep=pawrad(itypat)%rstep,lstep=pawrad(itypat)%lstep)
754 : : end if
755 : :
756 : : ! Loop over atoms of the type
757 [ + + ]: 30 : do iat=1,nattyp(itypat)
758 : 15 : iatm=iatm+1;iatom=atindx1(iatm)
759 : :
760 [ + + ]: 15 : if(option==2) corgr(:)=zero
761 : :
762 : : ! Consider each component in turn : compute range
763 [ + + ]: 60 : do mu=1,3
764 : : ! Convert reduced coord of given atom to [0,1)
765 : 45 : tau(mu)=mod(xred(mu,iatom)+one-aint(xred(mu,iatom)),one)
766 : : ! Use tau to find nearest grid point along R(mu)
767 : : ! (igrid=0 is the origin; shift by 1 to agree with usual index)
768 : 45 : igrid(mu)=nint(tau(mu)*real(ngfft(mu),dp))
769 : : ! Use range to compute an index range along R(mu)
770 : : ! (add 1 to make sure it covers full range)
771 : 60 : irange(mu)=1+nint((range/scale(mu))*real(ngfft(mu),dp))
772 : : end do
773 : :
774 : : ! Allocate arrays that depends on the range
775 [ + + ]: 60 : mrange=maxval(irange(1:3))
776 [ + - - + : 60 : ABI_MALLOC(ii,(2*mrange+1,3))
+ - + - +
- + - - +
- + ]
777 [ - + + - : 45 : ABI_MALLOC(rrdiff,(2*mrange+1,3))
+ - + - -
+ - + ]
778 : :
779 : : ! Set up counters that explore the relevant range of points around the atom
780 [ + - ]: 15 : if (geocode=='P') then
781 : : ! Fully periodic version
782 [ + + ]: 60 : do mu=1,3
783 : 45 : ishift=0
784 [ + + ]: 2790 : do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
785 : 2745 : ishift=ishift+1
786 : 2745 : ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu))
787 : 2790 : rrdiff(ishift,mu)=real(ixp,dp)/real(ngfft(mu),dp)-tau(mu)
788 : : end do
789 : 60 : ishiftmax(mu)=ishift
790 : : end do
791 : : else
792 : : ! Free or surface conditions
793 [ # # ]: 0 : do mu=1,3
794 : 0 : ishift=0
795 [ # # ]: 0 : do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
796 : 0 : call indpos_mkcore_alt(per(mu),ixp,ngfft(mu),jj,letsgo)
797 [ # # ]: 0 : if (letsgo) then
798 : 0 : ishift=ishift+1;ii(ishift,mu)=1+jj
799 : 0 : rrdiff(ishift,mu)=real(ixp,dp)/real(ngfft(mu),dp)-tau(mu)
800 : : end if
801 : : end do
802 : 0 : ishiftmax(mu)=ishift
803 : : end do
804 : : end if
805 : 15 : npts12=ishiftmax(1)*ishiftmax(2)
806 [ - + + - : 45 : ABI_MALLOC(indx1,(npts12))
+ - + - -
+ - + ]
807 [ - + + - : 30 : ABI_MALLOC(indx2,(npts12))
- + - + ]
808 [ - + + - : 30 : ABI_MALLOC(iindex,(npts12))
- + - + ]
809 [ - + + - : 45 : ABI_MALLOC(rnorm,(npts12))
+ - + - -
+ - + ]
810 [ + + ]: 15 : if (option==1.or.option==3) then
811 [ - + + - : 20 : ABI_MALLOC(tcore,(npts12))
- + - + -
+ ]
812 : : end if
813 [ + + ]: 15 : if (option>=2) then
814 [ - + + - : 20 : ABI_MALLOC(dtcore,(npts12))
- + - + -
+ ]
815 : : end if
816 [ - + ]: 15 : if (option==4) then
817 [ # # # # : 0 : ABI_MALLOC(d2tcore,(npts12))
# # # # #
# ]
818 : : end if
819 : :
820 : : ! Conduct loop over restricted range of grid points for iatom
821 [ + + ]: 930 : do ishift3=1,ishiftmax(3)
822 : 915 : i3=ii(ishift3,3)
823 : 915 : rdiff3=rrdiff(ishift3,3)
824 : :
825 [ + - ]: 915 : if(fftn3_distrib(i3)/=mpi_enreg%me_fft) cycle
826 : :
827 : : ! Select the vectors located around the current atom
828 : : ! TR: all of the following could be done inside or
829 : : ! outside the loops (i2,i1,i3).
830 : : ! Outside: the memory consumption increases.
831 : : ! Inside: the time of calculation increases.
832 : : ! Here, I choose to do it here, somewhere in the middle.
833 : 915 : npts=0
834 [ + + ]: 58698 : do ishift2=1,ishiftmax(2)
835 : 57783 : i2=ii(ishift2,2) ; rdiff2=rrdiff(ishift2,2)
836 : : rnorm2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2 &
837 : 57783 : & +2.0d0*rmet(3,2)*rdiff3*rdiff2
838 : 57783 : rnorm2_fact=2.0d0*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2)
839 [ + + ]: 3832485 : do ishift1=1,ishiftmax(1)
840 : 3773787 : i1=ii(ishift1,1) ; rdiff1=rrdiff(ishift1,1)
841 : 3773787 : rnorm2=rnorm2_part+rdiff1*(rnorm2_fact+rmet(1,1)*rdiff1)
842 : : ! Only accept contributions inside defined range
843 [ + + ]: 3831570 : if (rnorm2<range2-tol12) then
844 : 1315515 : npts=npts+1 ; iindex(npts)=npts
845 : 1315515 : indx1(npts)=ishift1;indx2(npts)=ishift2
846 : 1315515 : rnorm(npts)=sqrt(rnorm2)
847 : : end if
848 : : end do
849 : : end do
850 [ + + ]: 915 : if (npts==0) cycle
851 [ - + ]: 873 : if (npts>npts12) then
852 : 0 : message='npts>npts12!'
853 : 0 : ABI_BUG(message)
854 : : end if
855 : :
856 : : ! Evaluate core density (and derivatives) on the set of selected points
857 [ - + ]: 873 : if (usepaw==1) then
858 : : ! PAW: use splint routine
859 : 0 : call sort_dp(npts,rnorm(1:npts),iindex(1:npts),tol16)
860 [ # # ]: 0 : if (option==1.or.option==3) then
861 : : ! Evaluate fit of core density
862 : : call paw_splint(core_mesh%mesh_size,core_mesh%rad, &
863 : : & corespl(:,1),corespl(:,3),&
864 : 0 : & npts,rnorm(1:npts),tcore(1:npts))
865 : : end if
866 [ # # ]: 0 : if (option>=2) then
867 : : ! Evaluate fit of 1-der of core density
868 : : call paw_splint(core_mesh%mesh_size,core_mesh%rad, &
869 : : & corespl(:,2),corespl(:,4),&
870 : 0 : & npts,rnorm(1:npts),dtcore(1:npts))
871 : : end if
872 [ # # ]: 0 : if (option==4) then
873 : : ! Evaluate fit of 2nd-der of core density
874 : : call paw_splint(core_mesh%mesh_size,core_mesh%rad, &
875 : : & corespl(:,3),corespl(:,5),&
876 : 0 : & npts,rnorm(1:npts),d2tcore(1:npts))
877 : : end if
878 : : else
879 : : ! Norm-conserving PP:
880 : : ! Evaluate spline fit with method from Numerical Recipes
881 : : ! (p. 86 Numerical Recipes, Press et al;
882 : : ! NOTE error in book for sign of "aa" term in derivative)
883 [ + + ]: 1316388 : do ipts=1,npts
884 : 1315515 : yy=rnorm(ipts)*rangem1
885 : 1315515 : jj=1+int(yy*(n1xccc-1))
886 : 1315515 : diff=yy-(jj-1)*delta
887 : 1315515 : bb = diff*deltam1 ; aa = one-bb
888 : 1315515 : cc = aa*(aa**2-one)*delta2div6
889 : 1315515 : dd = bb*(bb**2-one)*delta2div6
890 [ + + ]: 1315515 : if (option==1.or.option==3) then
891 : : tcore(ipts)=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +&
892 : 877010 : & cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat)
893 : : end if
894 [ + + ]: 1315515 : if (option>=2) then
895 : : dtcore(ipts)=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +&
896 : 877010 : & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat)
897 : : end if
898 [ - + ]: 1316388 : if (option==4) then
899 : : d2tcore(ipts)=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +&
900 : 0 : & cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat)
901 : : end if
902 : : end do
903 : : end if
904 : :
905 : : ! Now, perform the loop over selected grid points
906 [ + + ]: 1316403 : do ipts=1,npts
907 : 1315515 : ishift1=indx1(iindex(ipts))
908 : 1315515 : ishift2=indx2(iindex(ipts))
909 : 1315515 : difmag=rnorm(ipts)
910 : :
911 : 1315515 : rdiff1=rrdiff(ishift1,1);rdiff2=rrdiff(ishift2,2)
912 : 1315515 : jpts=ii(ishift1,1)+n1*(ii(ishift2,2)-1+n2*(ffti3_local(i3)-1))
913 : :
914 : : ! === Evaluate charge density
915 [ + + ]: 1316430 : if (option==1) then
916 : 438505 : xccc3d(jpts)=xccc3d(jpts)+tcore(ipts)
917 : :
918 : : ! === Accumulate contributions to forces
919 [ + + ]: 877010 : else if (option==2) then
920 [ + + ]: 438505 : if (difmag>tol10) then
921 : 438500 : term=vxc_eff(jpts)*dtcore(ipts)/difmag
922 : 438500 : corgr(1)=corgr(1)+rdiff1*term
923 : 438500 : corgr(2)=corgr(2)+rdiff2*term
924 : 438500 : corgr(3)=corgr(3)+rdiff3*term
925 : : end if
926 : :
927 : : ! === Accumulate contributions to stress tensor (in red. coordinates)
928 [ + - ]: 438505 : else if (option==3) then
929 [ + + ]: 438505 : if (difmag>tol10) then
930 : 438500 : term=vxc_eff(jpts)*dtcore(ipts)*rangem1/difmag/real(nfftot,dp)
931 : : ! Write out the 6 symmetric components
932 : 438500 : corfra(1,1)=corfra(1,1)+term*rdiff1*rdiff1
933 : 438500 : corfra(2,2)=corfra(2,2)+term*rdiff2*rdiff2
934 : 438500 : corfra(3,3)=corfra(3,3)+term*rdiff3*rdiff3
935 : 438500 : corfra(3,2)=corfra(3,2)+term*rdiff3*rdiff2
936 : 438500 : corfra(3,1)=corfra(3,1)+term*rdiff3*rdiff1
937 : 438500 : corfra(2,1)=corfra(2,1)+term*rdiff2*rdiff1
938 : : ! (the above still needs to be transformed to cartesian coords)
939 : : end if
940 : : ! Also compute a diagonal term
941 : 438505 : strdia=strdia+vxc_eff(jpts)*tcore(ipts)
942 : :
943 : : ! === Compute frozen-wf contribution to Dynamical matrix
944 [ # # ]: 0 : else if (option==4) then
945 : 0 : tt(1)=rmet(1,1)*rdiff1+rmet(1,2)*rdiff2+rmet(1,3)*rdiff3
946 : 0 : tt(2)=rmet(2,1)*rdiff1+rmet(2,2)*rdiff2+rmet(2,3)*rdiff3
947 : 0 : tt(3)=rmet(3,1)*rdiff1+rmet(3,2)*rdiff2+rmet(3,3)*rdiff3
948 [ # # ]: 0 : if (difmag>tol10) then
949 : 0 : term=(ucvol/real(nfftot,dp))*vxc_eff(jpts)*rangem1/difmag
950 : 0 : term1=term*tcore(ipts)
951 : 0 : term2=term*d2tcore(ipts)*rangem1/difmag
952 [ # # ]: 0 : do mu=1,3
953 [ # # ]: 0 : do nu=1,3
954 : : dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)&
955 : : & +(term2-term1/difmag**2)*tt(mu)*tt(nu)&
956 : 0 : & +term1*rmet(mu,nu)
957 : : end do
958 : : end do
959 : : else
960 : : ! There is a contribution from difmag=zero !
961 : 0 : term=(ucvol/real(nfftot,dp))*vxc_eff(jpts)*d2tcore(ipts)*rangem1**2
962 [ # # ]: 0 : do mu=1,3
963 [ # # ]: 0 : do nu=1,3
964 : 0 : dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)+term*rmet(mu,nu)
965 : : end do
966 : : end do
967 : : end if
968 : : end if ! Choice of option
969 : :
970 : : end do ! Loop on ipts (ishift1, ishift2)
971 : :
972 : : end do ! Loop on ishift3
973 : :
974 [ - + ]: 15 : ABI_FREE(ii)
975 [ - + ]: 15 : ABI_FREE(rrdiff)
976 [ - + ]: 15 : ABI_FREE(indx1)
977 [ - + ]: 15 : ABI_FREE(indx2)
978 [ - + ]: 15 : ABI_FREE(iindex)
979 [ - + ]: 15 : ABI_FREE(rnorm)
980 [ + + ]: 15 : if (allocated(tcore)) then
981 : 10 : ABI_FREE(tcore)
982 : : end if
983 [ - + ]: 15 : if (allocated(tcoretau)) then
984 : 0 : ABI_FREE(tcoretau)
985 : : end if
986 [ + + ]: 15 : if (allocated(dtcore)) then
987 : 10 : ABI_FREE(dtcore)
988 : : end if
989 [ - + ]: 15 : if (allocated(d2tcore)) then
990 : 0 : ABI_FREE(d2tcore)
991 : : end if
992 : :
993 [ + + ]: 30 : if (option==2) then
994 : 5 : fact=-(ucvol/real(nfftot,dp))/range
995 [ + + ]: 20 : grxc(:,iatom)=corgr(:)*fact
996 : : end if
997 : :
998 : : ! End loop on atoms
999 : : end do
1000 : :
1001 [ + - ]: 61 : if (usepaw==1) then
1002 : 0 : call pawrad_free(core_mesh)
1003 : : end if
1004 : :
1005 : : !End loop over atom types
1006 : : end do
1007 : :
1008 [ + + - + ]: 46 : if(option>=2.and.nspden>=2) then
1009 [ # # ]: 0 : ABI_FREE(vxc_eff)
1010 : : end if
1011 : :
1012 : : !Forces: translate into reduced coordinates
1013 [ + + ]: 46 : if (option==2) then
1014 [ + + ]: 8 : do iatom=1,natom
1015 : 5 : t1=grxc(1,iatom);t2=grxc(2,iatom);t3=grxc(3,iatom)
1016 [ + + ]: 23 : grxc(:,iatom)=rmet(:,1)*t1+rmet(:,2)*t2+rmet(:,3)*t3
1017 : : end do
1018 : : end if
1019 : :
1020 : : !Stress tensor: symmetrize, translate into cartesian coord., add diagonal part
1021 [ + + ]: 46 : if (option==3) then
1022 : 3 : corstr(1)=corfra(1,1) ; corstr(2)=corfra(2,2)
1023 : 3 : corstr(3)=corfra(3,3) ; corstr(4)=corfra(3,2)
1024 : 3 : corstr(5)=corfra(3,1) ; corstr(6)=corfra(2,1)
1025 : 3 : call strconv(corstr,rprimd,corstr)
1026 : 3 : corstr(1)=corstr(1)+strdia/real(nfftot,dp)
1027 : 3 : corstr(2)=corstr(2)+strdia/real(nfftot,dp)
1028 : 3 : corstr(3)=corstr(3)+strdia/real(nfftot,dp)
1029 : : end if
1030 : :
1031 : : !If needed sum over MPI processes
1032 [ - + ]: 46 : if(mpi_enreg%nproc_fft>1) then
1033 : 0 : call timab(539,1,tsec)
1034 [ # # ]: 0 : if (option==2) then
1035 : 0 : call xmpi_sum(grxc,mpi_enreg%comm_fft,ier)
1036 : : end if
1037 [ # # ]: 0 : if (option==3) then
1038 : 0 : call xmpi_sum(corstr,mpi_enreg%comm_fft,ier)
1039 : : end if
1040 [ # # ]: 0 : if (option==4) then
1041 : 0 : call xmpi_sum(dyfrx2,mpi_enreg%comm_fft,ier)
1042 : : end if
1043 : 0 : call timab(539,2,tsec)
1044 : : end if
1045 : :
1046 [ - + - + : 92 : call timab(12,2,tsec)
- + - + ]
1047 : :
1048 : : contains
1049 : : !!***
1050 : :
1051 : : !--------------------------------------------------------------
1052 : :
1053 : : !!****f* ABINIT/cross_mkcore_alt
1054 : : !! NAME
1055 : : !! cross_mkcore_alt
1056 : : !!
1057 : : !! FUNCTION
1058 : : !! Define magnitude of cross product of two vectors
1059 : : !!
1060 : : !! SOURCE
1061 : :
1062 : 138 : function cross_mkcore_alt(xx,yy,zz,aa,bb,cc)
1063 : :
1064 : : real(dp) :: cross_mkcore_alt
1065 : : real(dp),intent(in) :: xx,yy,zz,aa,bb,cc
1066 : : ! *************************************************************************
1067 : 138 : cross_mkcore_alt=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2)
1068 : 138 : end function cross_mkcore_alt
1069 : : !!***
1070 : :
1071 : : !--------------------------------------------------------------
1072 : :
1073 : : !!****f* ABINIT/indpos_mkcore_alt
1074 : : !! NAME
1075 : : !! indpos_mkcore_alt
1076 : : !!
1077 : : !! FUNCTION
1078 : : !! Find the grid index of a given position in the cell according to the BC
1079 : : !! Determine also whether the index is inside or outside the box for free BC
1080 : : !!
1081 : : !! SOURCE
1082 : :
1083 : 0 : subroutine indpos_mkcore_alt(periodic,ii,nn,jj,inside)
1084 : : ! Find the grid index of a given position in the cell according to the BC
1085 : : ! Determine also whether the index is inside or outside the box for free BC
1086 : : integer, intent(in) :: ii,nn
1087 : : integer, intent(out) :: jj
1088 : : logical, intent(in) :: periodic
1089 : : logical, intent(out) :: inside
1090 : : ! *************************************************************************
1091 [ # # ]: 0 : if (periodic) then
1092 : 0 : inside=.true. ; jj=modulo(ii-1,nn)+1
1093 : : else
1094 [ # # # # ]: 0 : jj=ii ; inside=(ii>=1.and.ii<=nn)
1095 : : end if
1096 : 0 : end subroutine indpos_mkcore_alt
1097 : :
1098 : : end subroutine mkcore_alt
1099 : : !!***
1100 : :
1101 : : !!****f* ABINIT/dfpt_mkcore
1102 : : !! NAME
1103 : : !! dfpt_mkcore
1104 : : !!
1105 : : !! FUNCTION
1106 : : !! Compute the derivative of the core electron density
1107 : : !! with respect to one specific atom displacement
1108 : : !! In case of derivative with respect to k or
1109 : : !! electric (magnetic) field perturbation, the 1st-order core electron density
1110 : : !! vanishes.
1111 : : !!
1112 : : !! INPUTS
1113 : : !! cplex: if 1, real space 1-order functions on FFT grid are REAL,
1114 : : !! if 2, COMPLEX
1115 : : !! idir=direction of atomic displacement (=1,2 or 3 : displacement of
1116 : : !! atom ipert along the 1st, 2nd or 3rd axis) or cartesian coordinate
1117 : : !! pair for strain perturbation
1118 : : !! ipert=number of the atom being displaced or natom+3,4 for strain
1119 : : !! perturbation
1120 : : !! natom=number of atoms in cell.
1121 : : !! ntypat=number of types of atoms in cell.
1122 : : !! n1,n2,n3=fft grid dimensions.
1123 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
1124 : : !! qphon(3)=wavevector of the phonon
1125 : : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
1126 : : !! typat(natom)=integer type for each atom in cell
1127 : : !! ucvol=unit cell volume (bohr**3).
1128 : : !! xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
1129 : : !! xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives,
1130 : : !! for each type of atom, from psp
1131 : : !! xred(3,natom)=reduced coordinates for atoms in unit cell
1132 : : !!
1133 : : !! OUTPUT
1134 : : !! xccc3d1(cplex*n1*n2*n3)=3D core electron density for XC core correction, bohr^-3
1135 : : !!
1136 : : !! NOTES
1137 : : !! Note that this routine is tightly connected to the mkcore.f routine
1138 : : !!
1139 : : !! SOURCE
1140 : :
1141 : 14692 : subroutine dfpt_mkcore(cplex,idir,ipert,natom,ntypat,n1,n1xccc,&
1142 : 14692 : & n2,n3,qphon,rprimd,typat,ucvol,xcccrc,xccc1d,xccc3d1,xred)
1143 : :
1144 : : !Arguments ------------------------------------
1145 : : !scalars
1146 : : integer,intent(in) :: cplex,idir,ipert,n1,n1xccc,n2,n3,natom,ntypat
1147 : : real(dp),intent(in) :: ucvol
1148 : : !arrays
1149 : : integer,intent(in) :: typat(natom)
1150 : : real(dp),intent(in) :: qphon(3),rprimd(3,3),xccc1d(n1xccc,6,ntypat)
1151 : : real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom)
1152 : : real(dp),intent(out) :: xccc3d1(cplex*n1*n2*n3)
1153 : :
1154 : : !Local variables-------------------------------
1155 : : !scalars
1156 : : integer,parameter :: mshift=401
1157 : : integer :: i1,i2,i3,iatom,ifft,ishift,ishift1,ishift2,ishift3,istr
1158 : : integer :: ixp,jj,ka,kb,mrange,mu,nu
1159 : : real(dp) :: aa,bb,cc,dd,delta,delta2div6,deltam1,diff,difmag
1160 : : real(dp) :: difmag2,difmag2_fact,difmag2_part,func,phase,phi,phr,prod
1161 : : real(dp) :: range,range2,rangem1,rdiff1,rdiff2,rdiff3,term
1162 : : real(dp) :: yy
1163 : : character(len=500) :: message
1164 : : !arrays
1165 : : integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1166 : : integer :: igrid(3),irange(3),ngfft(3)
1167 : 14692 : integer,allocatable :: ii(:,:)
1168 : : real(dp) :: drmetds(3,3),lencp(3),rmet(3,3),scale(3),tau(3)
1169 : 14692 : real(dp),allocatable :: rrdiff(:,:)
1170 : : ! *************************************************************************
1171 : :
1172 : : ! if( ipert<1 .or. ipert> natom+7) then
1173 : : ! write(message,'(a,i0,a,a,a,i0,a)')&
1174 : : !& ' The argument ipert must be between 1 and natom+7=',natom+7,',',ch10,&
1175 : : !& ' while it is ipert=',ipert,'.'
1176 : : ! ABI_BUG(message)
1177 : : ! end if
1178 : :
1179 [ + + + + : 14692 : if( (ipert==natom+3 .or. ipert==natom+4) .and. cplex/=1) then
- + ]
1180 : : write(message,'(3a,i4,a)')&
1181 : 0 : & 'The argument cplex must be 1 for strain perturbationh',ch10,&
1182 : 0 : & 'while it is cplex=',cplex,'.'
1183 : 0 : ABI_BUG(message)
1184 : : end if
1185 : :
1186 : : !Zero out array
1187 [ + + ]: 274243017 : xccc3d1(:)=0.0_dp
1188 : :
1189 : : !For a non-linear XC core correction, the perturbation must be phonon-type or strain type
1190 [ + + + + : 14692 : if(ipert<=natom .or. ipert==natom+3 .or. ipert==natom+4) then
+ + ]
1191 : :
1192 [ - + ]: 14347 : if( idir<1 .or. idir> 3) then
1193 : : write(message,'(a,a,a,i4,a)')&
1194 : 0 : & 'The argument idir must be between 1 and 3,',ch10,&
1195 : 0 : & 'while it is idir=',idir,'.'
1196 : 0 : ABI_BUG(message)
1197 : : end if
1198 : :
1199 : : ! Compute lengths of cross products for pairs of primitive
1200 : : ! translation vectors (used in setting index search range below)
1201 : : lencp(1)=cross_mk(rprimd(1,2),rprimd(2,2),rprimd(3,2),&
1202 : 14347 : & rprimd(1,3),rprimd(2,3),rprimd(3,3))
1203 : : lencp(2)=cross_mk(rprimd(1,3),rprimd(2,3),rprimd(3,3),&
1204 : 14347 : & rprimd(1,1),rprimd(2,1),rprimd(3,1))
1205 : : lencp(3)=cross_mk(rprimd(1,1),rprimd(2,1),rprimd(3,1),&
1206 : 14347 : & rprimd(1,2),rprimd(2,2),rprimd(3,2))
1207 : :
1208 : : ! Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3
1209 : : ! (recall ucvol=R1.(R2xR3))
1210 [ + + ]: 57388 : scale(:)=ucvol/lencp(:)
1211 : :
1212 : : ! Compute metric tensor in real space rmet
1213 [ + + ]: 57388 : do nu=1,3
1214 : : rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+&
1215 [ + + ]: 186511 : & rprimd(3,:)*rprimd(3,nu)
1216 : : end do
1217 : :
1218 : : ! Section to be executed only for strain perturbation
1219 : : ! Compute derivative of metric tensor wrt strain component istr
1220 [ + + + + ]: 14347 : if(ipert==natom+3 .or. ipert==natom+4) then
1221 : 502 : istr=idir + 3*(ipert-natom-3)
1222 : :
1223 : 502 : ka=idx(2*istr-1);kb=idx(2*istr)
1224 [ + + ]: 2008 : do jj = 1,3
1225 [ + + ]: 6526 : drmetds(:,jj)=(rprimd(ka,:)*rprimd(kb,jj)+rprimd(kb,:)*rprimd(ka,jj))
1226 : : end do
1227 : : ! For historical reasons:
1228 [ + + + + ]: 20371 : drmetds(:,:)=0.5_dp*drmetds(:,:)
1229 : :
1230 : : ! end of strain perturbation section
1231 : : end if
1232 : :
1233 : 14347 : ngfft(1)=n1
1234 : 14347 : ngfft(2)=n2
1235 : 14347 : ngfft(3)=n3
1236 : :
1237 : 14347 : delta=1.0_dp/(n1xccc-1)
1238 : 14347 : deltam1=n1xccc-1
1239 : 14347 : delta2div6=delta**2/6.0_dp
1240 : :
1241 : : ! Loop over atoms in unit cell
1242 : : ! Note that we cycle immediately for all except the displaced atom
1243 : : ! for such a perturbation. The loop is executed over all the
1244 : : ! atoms for a strain peturbation.
1245 [ + + ]: 56095 : do iatom=1,natom
1246 [ + + + + ]: 41748 : if(ipert<=natom .and. iatom/=ipert) cycle
1247 : : ! Set search range (density cuts off perfectly beyond range)
1248 : : ! Cycle if no range.
1249 : 14878 : range=0.0_dp
1250 : 14878 : range=xcccrc(typat(iatom))
1251 [ + + ]: 14878 : if(range<1.d-16) cycle
1252 : :
1253 : 13644 : range2=range**2
1254 : 13644 : rangem1=1.0_dp/range
1255 : :
1256 : : ! compute mrange for ii(:,3), rrdiff(:,3), inserted by MM (2005/12/06)
1257 : : ! Consider each component in turn : compute range
1258 [ + + ]: 54576 : do mu=1,3
1259 : :
1260 : : ! Convert reduced coord of given atom to [0,1)
1261 : 40932 : tau(mu)=mod(xred(mu,iatom)+1._dp-aint(xred(mu,iatom)),1._dp)
1262 : :
1263 : : ! Use tau to find nearest grid point along R(mu)
1264 : : ! (igrid=0 is the origin; shift by 1 to agree with usual index)
1265 : 40932 : igrid(mu)=nint(tau(mu)*dble(ngfft(mu)))
1266 : :
1267 : : ! Use range to compute an index range along R(mu)
1268 : : ! (add 1 to make sure it covers full range)
1269 : 54576 : irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu)))
1270 : :
1271 : : end do
1272 : :
1273 : : ! Allocate arrays that depends on the range
1274 [ + + ]: 54576 : mrange=maxval(irange(1:3))
1275 [ + - - + : 54576 : ABI_MALLOC(ii,(2*mrange+1,3))
+ - + - +
- + - - +
- + ]
1276 [ - + + - : 40932 : ABI_MALLOC(rrdiff,(2*mrange+1,3))
+ - + - -
+ - + ]
1277 : :
1278 : : ! Consider each component in turn
1279 [ + + ]: 54576 : do mu=1,3
1280 : :
1281 : : ! temporarily suppressed by MM (2005/12/02)
1282 : : ! Convert reduced coord of given atom to [0,1)
1283 : : ! tau(mu)=mod(xred(mu,iatom)+1._dp-aint(xred(mu,iatom)),1._dp)
1284 : :
1285 : : ! Use tau to find nearest grid point along R(mu)
1286 : : ! (igrid=0 is the origin; shift by 1 to agree with usual index)
1287 : : ! igrid(mu)=nint(tau(mu)*dble(ngfft(mu)))
1288 : :
1289 : : ! Use range to compute an index range along R(mu)
1290 : : ! (add 1 to make sure it covers full range)
1291 : : ! irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu)))
1292 : :
1293 : : ! Check that the largest range is smallest than the maximum
1294 : : ! allowed one
1295 : : ! if(2*irange(mu)+1 > mshift)then
1296 : : ! write(message, '(a,a,a,a,i6,a)' ) ch10,&
1297 : : ! & ' dfpt_mkcore : BUG -',ch10,&
1298 : : ! & ' The range around atom',iatom,' is too large.'
1299 : : ! ABI_BUG(message)
1300 : : ! end if
1301 : :
1302 : : ! Set up a counter that explore the relevant range
1303 : : ! of points around the atom
1304 : 40932 : ishift=0
1305 [ + + ]: 1688176 : do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
1306 : 1633600 : ishift=ishift+1
1307 : 1633600 : ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu))
1308 : 1674532 : rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu)
1309 : : end do
1310 : :
1311 : : ! End loop on mu
1312 : : end do
1313 : :
1314 : : ! Conduct triple loop over restricted range of grid points for iatom
1315 : :
1316 [ + + ]: 545288 : do ishift3=1,1+2*irange(3)
1317 : : ! map back to [1,ngfft(3)] for usual fortran index in unit cell
1318 : 531644 : i3=ii(ishift3,3)
1319 : : ! find vector from atom location to grid point (reduced)
1320 : 531644 : rdiff3=rrdiff(ishift3,3)
1321 : :
1322 [ + + ]: 28801076 : do ishift2=1,1+2*irange(2)
1323 : 28255788 : i2=ii(ishift2,2)
1324 : 28255788 : rdiff2=rrdiff(ishift2,2)
1325 : : ! Prepare the computation of difmag2
1326 : : difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2&
1327 : 28255788 : & +2.0_dp*rmet(3,2)*rdiff3*rdiff2
1328 : 28255788 : difmag2_fact=2.0_dp*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2)
1329 : :
1330 [ + + ]: 1914049192 : do ishift1=1,1+2*irange(1)
1331 : 1885261760 : rdiff1=rrdiff(ishift1,1)
1332 : :
1333 : : ! Compute (rgrid-tau-Rprim)**2
1334 : 1885261760 : difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1)
1335 : :
1336 : : ! Only accept contribution inside defined range
1337 [ + + ]: 1913517548 : if (difmag2<range2) then
1338 : :
1339 : : ! Prepare computation of core charge function and derivative,
1340 : : ! using splines
1341 : 677162857 : difmag=sqrt(difmag2)
1342 [ + + ]: 677162857 : if (difmag>=1.0d-10) then
1343 : 677154932 : i1=ii(ishift1,1)
1344 : 677154932 : yy=difmag*rangem1
1345 : :
1346 : : ! Compute index of yy over 1 to n1xccc scale
1347 : 677154932 : jj=1+int(yy*(n1xccc-1))
1348 : 677154932 : diff=yy-(jj-1)*delta
1349 : :
1350 : : ! Will evaluate spline fit (p. 86 Numerical Recipes, Press et al;
1351 : : ! NOTE error in book for sign of "aa" term in derivative;
1352 : : ! also see splfit routine).
1353 : 677154932 : bb = diff*deltam1
1354 : 677154932 : aa = 1.0_dp-bb
1355 : 677154932 : cc = aa*(aa**2-1.0_dp)*delta2div6
1356 : 677154932 : dd = bb*(bb**2-1.0_dp)*delta2div6
1357 : :
1358 : : ! Evaluate spline fit of 1st der of core charge density
1359 : : ! from xccc1d(:,2,:) and (:,4,:)
1360 : : func=aa*xccc1d(jj,2,typat(iatom))+bb*xccc1d(jj+1,2,typat(iatom)) +&
1361 : 677154932 : & cc*xccc1d(jj,4,typat(iatom))+dd* xccc1d(jj+1,4,typat(iatom))
1362 : :
1363 [ + + ]: 677154932 : if(ipert<=natom) then
1364 : : phase=2*pi*(qphon(1)*(rdiff1+xred(1,iatom)) &
1365 : : & +qphon(2)*(rdiff2+xred(2,iatom)) &
1366 : 664960548 : & +qphon(3)*(rdiff3+xred(3,iatom)))
1367 : 664960548 : prod=rmet(idir,1)*rdiff1+rmet(idir,2)*rdiff2+rmet(idir,3)*rdiff3
1368 : :
1369 : 664960548 : term=-func*rangem1/difmag*prod
1370 : 664960548 : ifft=i1+n1*(i2-1+n2*(i3-1))
1371 : 664960548 : phr=cos(phase)
1372 [ + + ]: 664960548 : if(cplex==1)then
1373 : 252397494 : xccc3d1(ifft)=xccc3d1(ifft)+term*phr
1374 : : else
1375 : 412563054 : phi=sin(phase)
1376 : 412563054 : xccc3d1(2*ifft-1)=xccc3d1(2*ifft-1)+term*phr
1377 : 412563054 : xccc3d1(2*ifft )=xccc3d1(2*ifft )-term*phi
1378 : : end if
1379 : : else
1380 : : prod=&
1381 : : & (rdiff1*(drmetds(1,1)*rdiff1+drmetds(1,2)*rdiff2+drmetds(1,3)*rdiff3)&
1382 : : & +rdiff2*(drmetds(2,1)*rdiff1+drmetds(2,2)*rdiff2+drmetds(2,3)*rdiff3)&
1383 : 12194384 : & +rdiff3*(drmetds(3,1)*rdiff1+drmetds(3,2)*rdiff2+drmetds(3,3)*rdiff3))
1384 : 12194384 : term=prod*func*rangem1/difmag
1385 : :
1386 : 12194384 : ifft=i1+n1*(i2-1+n2*(i3-1))
1387 : 12194384 : xccc3d1(ifft)=xccc3d1(ifft)+term
1388 : :
1389 : : end if
1390 : :
1391 : : ! End of the condition for the distance not to vanish
1392 : : end if
1393 : :
1394 : : ! End of condition to be inside the range
1395 : : end if
1396 : :
1397 : : ! End loop on ishift1
1398 : : end do
1399 : :
1400 : : ! End loop on ishift2
1401 : : end do
1402 : :
1403 : : ! End loop on ishift3
1404 : : end do
1405 : :
1406 [ - + ]: 13644 : ABI_FREE(ii)
1407 [ - + ]: 56440 : ABI_FREE(rrdiff)
1408 : : ! End loop on atoms
1409 : : end do
1410 : :
1411 : : ! End of the condition ipert corresponds to a phonon type perturbation
1412 : : ! or strain type perturbation
1413 : : end if
1414 : :
1415 : : contains
1416 : :
1417 : 43041 : function cross_mk(xx,yy,zz,aa,bb,cc)
1418 : :
1419 : : real(dp) :: cross_mk
1420 : : real(dp),intent(in) :: xx,yy,zz,aa,bb,cc
1421 : 43041 : cross_mk=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2)
1422 : 43041 : end function cross_mk
1423 : :
1424 : : end subroutine dfpt_mkcore
1425 : : !!***
1426 : :
1427 : : end module m_mkcore
1428 : : !!***
|
