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 9509 : subroutine mkcore(corstr,dyfrx2,grxc,mpi_enreg,natom,nfft,nspden,ntypat,n1,n1xccc,&
101 9509 : & 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 9509 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
127 9509 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
128 9509 : 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 9509 : real(dp),allocatable :: rrdiff(:,:),work(:,:,:)
132 : !************************************************************************
133 :
134 9509 : call timab(12,1,tsec)
135 :
136 : !Make sure option is acceptable
137 9509 : 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 9509 : if (option==1) then
148 : ! Zero out array to permit accumulation over atom types below:
149 49390449 : xccc3d(:)=zero
150 6148 : else if (option==2) then
151 : ! Zero out gradient of Exc array
152 63261 : grxc(:,:)=zero
153 2055 : else if (option==3) then
154 : ! Zero out locally defined stress array
155 9509 : corfra(:,:)=zero
156 9509 : 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 9509 : & rprimd(1,3),rprimd(2,3),rprimd(3,3))
168 : lencp(2)=cross_mkcore(rprimd(1,3),rprimd(2,3),rprimd(3,3),&
169 9509 : & rprimd(1,1),rprimd(2,1),rprimd(3,1))
170 : lencp(3)=cross_mkcore(rprimd(1,1),rprimd(2,1),rprimd(3,1),&
171 9509 : & 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 38036 : scale(:)=ucvol/lencp(:)
176 :
177 : !Compute metric tensor in real space rmet
178 38036 : do nu=1,3
179 123617 : rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+rprimd(3,:)*rprimd(3,nu)
180 : end do
181 :
182 9509 : ngfft(1)=n1
183 9509 : ngfft(2)=n2
184 9509 : ngfft(3)=n3
185 9509 : nfftot=n1*n2*n3
186 9509 : me_fft = mpi_enreg%me_fft
187 :
188 : ! Get the distrib associated with this fft_grid
189 9509 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
190 :
191 9509 : delta=one/(n1xccc-1)
192 9509 : deltam1=n1xccc-1
193 9509 : delta2div6=delta**2/6.0d0
194 :
195 9509 : if (option>=2) then
196 30740 : 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 6148 : if (nspden>=2) then
200 : ifft=1
201 24804 : do i3=1,n3
202 24804 : if(me_fft==fftn3_distrib(i3)) then
203 573400 : do i2=1,n2
204 14211550 : do i1=1,n1
205 13638150 : work(i1,i2,i3)=half*(vxc(ifft,1)+vxc(ifft,2))
206 14187781 : ifft=ifft+1
207 : end do
208 : end do
209 : end if
210 : end do
211 : else
212 : ifft=1
213 130122 : do i3=1,n3
214 130122 : if(me_fft==fftn3_distrib(i3)) then
215 3440853 : do i2=1,n2
216 126830547 : do i1=1,n1
217 123389694 : work(i1,i2,i3)=vxc(ifft,1)
218 126706338 : 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 43644 : do iatom=1,natom
229 :
230 34135 : if(option==2)then
231 14792 : grxc1=zero
232 14792 : grxc2=zero
233 14792 : grxc3=zero
234 : end if
235 :
236 : ! Set search range (density cuts off perfectly beyond range)
237 34135 : itypat=typat(iatom)
238 34135 : range=xcccrc(itypat)
239 :
240 : ! Skip loop if this atom has no core charge
241 34135 : if (abs(range)<1.d-16) cycle
242 :
243 31789 : range2=range**2
244 31789 : rangem1=one/range
245 :
246 : ! Consider each component in turn : compute range
247 127156 : do mu=1,3
248 :
249 : ! Convert reduced coord of given atom to [0,1)
250 95367 : 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 95367 : 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 127156 : irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu)))
259 :
260 : end do
261 :
262 : ! Allocate arrays that depends on the range
263 127156 : mrange=maxval(irange(1:3))
264 127156 : ABI_MALLOC(ii,(2*mrange+1,3))
265 95367 : ABI_MALLOC(rrdiff,(2*mrange+1,3))
266 :
267 : ! Set up counters that explore the relevant range
268 : ! of points around the atom
269 127156 : do mu=1,3
270 95367 : ishift=0
271 2919937 : do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
272 2792781 : ishift=ishift+1
273 2792781 : ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu))
274 2888148 : 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 953024 : do ishift3=1,1+2*irange(3)
280 : ! map back to [1,ngfft(3)] for usual fortran index in unit cell
281 921235 : i3=ii(ishift3,3)
282 921235 : if(fftn3_distrib(i3)/=mpi_enreg%me_fft) cycle
283 : ! find vector from atom location to grid point (reduced)
284 912745 : rdiff3=rrdiff(ishift3,3)
285 :
286 32795547 : do ishift2=1,1+2*irange(2)
287 31851013 : i2=ii(ishift2,2)
288 31851013 : rdiff2=rrdiff(ishift2,2)
289 : ! Prepare the computation of difmag2
290 : difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2&
291 31851013 : & +2.0d0*rmet(3,2)*rdiff3*rdiff2
292 31851013 : difmag2_fact=2.0d0*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2)
293 :
294 1402719161 : do ishift1=1,1+2*irange(1)
295 1369946913 : rdiff1=rrdiff(ishift1,1)
296 :
297 : ! Compute (rgrid-tau-Rprim)**2
298 1369946913 : difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1)
299 :
300 : ! Only accept contribution inside defined range
301 1401797926 : if (difmag2<range2-tol12) then
302 :
303 : ! Prepare computation of core charge function and derivative,
304 : ! using splines
305 498513981 : i1=ii(ishift1,1)
306 498513981 : difmag=sqrt(difmag2)
307 498513981 : yy=difmag*rangem1
308 :
309 : ! Compute index of yy over 1 to n1xccc scale
310 498513981 : jj=1+int(yy*(n1xccc-1))
311 498513981 : 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 498513981 : bb = diff*deltam1
317 498513981 : aa = one-bb
318 498513981 : cc = aa*(aa**2-one)*delta2div6
319 498513981 : dd = bb*(bb**2-one)*delta2div6
320 :
321 :
322 : ! Test first for option 2, the most frequently used
323 498513981 : if (option==2) then
324 :
325 : ! Accumulate contributions to Exc gradients
326 :
327 167786929 : 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 167784009 : & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat)
333 167784009 : term=work(i1,i2,i3)*func/difmag
334 167784009 : grxc1=grxc1+rdiff1*term
335 167784009 : grxc2=grxc2+rdiff2*term
336 167784009 : grxc3=grxc3+rdiff3*term
337 : end if
338 :
339 330727052 : 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 216671601 : & cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat)
345 :
346 : ! Accumulate contributions to core electron density
347 : ! throughout unit cell
348 216671601 : ifft=i1+n1*(i2-1+n2*(ffti3_local(i3)-1))
349 216671601 : xccc3d(ifft)=xccc3d(ifft)+func
350 :
351 114055451 : else if (option==3) then
352 :
353 : ! Accumulate contributions to stress tensor
354 : ! in reduced coordinates
355 :
356 92318265 : 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 92316670 : & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat)
362 92316670 : term=work(i1,i2,i3)*func*rangem1/difmag/dble(n1*n2*n3)
363 : ! Write out the 6 symmetric components
364 92316670 : corfra(1,1)=corfra(1,1)+term*rdiff1**2
365 92316670 : corfra(2,2)=corfra(2,2)+term*rdiff2**2
366 92316670 : corfra(3,3)=corfra(3,3)+term*rdiff3**2
367 92316670 : corfra(3,2)=corfra(3,2)+term*rdiff3*rdiff2
368 92316670 : corfra(3,1)=corfra(3,1)+term*rdiff3*rdiff1
369 92316670 : 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 92318265 : & cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat)
379 92318265 : 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 31789 : ABI_FREE(ii)
440 31789 : ABI_FREE(rrdiff)
441 :
442 41298 : if(option==2)then
443 13381 : fact=-(ucvol/dble(nfftot))/range
444 13381 : grxc(1,iatom)=grxc1*fact
445 13381 : grxc(2,iatom)=grxc2*fact
446 13381 : grxc(3,iatom)=grxc3*fact
447 : end if
448 :
449 : end do ! End big loop on atoms
450 :
451 9509 : if (option==2) then
452 :
453 : ! Apply rmet as needed to get reduced coordinate gradients
454 18885 : do iatom=1,natom
455 14792 : t1=grxc(1,iatom)
456 14792 : t2=grxc(2,iatom)
457 14792 : t3=grxc(3,iatom)
458 63261 : grxc(:,iatom)=rmet(:,1)*t1+rmet(:,2)*t2+rmet(:,3)*t3
459 :
460 : end do
461 : end if
462 :
463 9509 : if (option==3) then
464 :
465 : ! Transform stress tensor from full storage mode to symmetric storage mode
466 1809 : corstr(1)=corfra(1,1)
467 1809 : corstr(2)=corfra(2,2)
468 1809 : corstr(3)=corfra(3,3)
469 1809 : corstr(4)=corfra(3,2)
470 1809 : corstr(5)=corfra(3,1)
471 1809 : corstr(6)=corfra(2,1)
472 :
473 : ! Transform stress tensor from reduced coordinates to cartesian coordinates
474 1809 : 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 1809 : ifft=0 ; strdia=zero
479 41324 : do i3=1,n3
480 41324 : if(me_fft==fftn3_distrib(i3)) then
481 944282 : do i2=1,n2
482 27526212 : do i1=1,n1
483 26581930 : ifft=ifft+1
484 27487097 : strdia=strdia+work(i1,i2,i3)*xccc3d(ifft)
485 : end do
486 : end do
487 : end if
488 : end do
489 1809 : strdia=strdia/dble(nfftot)
490 : ! strdia=DDOT(nfft,work,1,xccc3d,1)/dble(nfftot)
491 :
492 : ! Add diagonal term to stress tensor
493 1809 : corstr(1)=corstr(1)+strdia
494 1809 : corstr(2)=corstr(2)+strdia
495 1809 : corstr(3)=corstr(3)+strdia
496 : end if
497 :
498 9509 : if(option>=2) then
499 6148 : ABI_FREE(work)
500 : end if
501 :
502 9509 : 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 19018 : 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 28527 : 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 28527 : cross_mkcore=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2)
536 28527 : 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 : !!***
|
