Branch data Line data Source code
1 : : !!****m* ABINIT/m_dfpt_mkvxcstr
2 : : !! NAME
3 : : !! m_dfpt_mkvxcstr
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2001-2026 ABINIT group (DRH,XG)
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_dfpt_mkvxcstr
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_abicore
27 : :
28 : : use defs_abitypes, only : MPI_type
29 : : use m_time, only : timab
30 : : use m_matrix, only : matr3inv
31 : : use m_xctk, only : xcden, xcpot
32 : :
33 : : implicit none
34 : :
35 : : private
36 : : !!***
37 : :
38 : : public :: dfpt_mkvxcstr
39 : : !!***
40 : :
41 : : contains
42 : : !!***
43 : :
44 : : !!****f* ABINIT/dfpt_mkvxcstr
45 : : !! NAME
46 : : !! dfpt_mkvxcstr
47 : : !!
48 : : !! FUNCTION
49 : : !! Compute the first-order change of exchange-correlation potential
50 : : !! due to strain: assemble the first-order density change with the
51 : : !! frozen-core density change, then use the exchange-correlation kernel.
52 : : !!
53 : : !! INPUTS
54 : : !! cplex= if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
55 : : !! idir=direction of the current perturbation
56 : : !! ipert=type of the perturbation
57 : : !! kxc(nfft,nkxc)=exchange and correlation kernel (see rhotoxc.f)
58 : : !! mpi_enreg=information about MPI parallelization
59 : : !! natom=number of atoms in cell.
60 : : !! nfft=(effective) number of FFT grid points (for this processor)
61 : : !! ngfft(18)=contain all needed information about 3D FFT,
62 : : !! see ~abinit/doc/variables/vargs.htm#ngfft
63 : : !! nhat(nfft,nspden*nhatdim)= -PAW only- compensation density
64 : : !! nhat1(cplex*nfft,2nspden*usepaw)= -PAW only- 1st-order compensation density
65 : : !! nkxc=second dimension of the kxc array
66 : : !! non_magnetic_xc= if true, handle density/potential as non-magnetic (even if it is)
67 : : !! nspden=number of spin-density components
68 : : !! n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used, otherwise, nfft
69 : : !! option=if 0, work only with strain-derivative frozen-wavefunction
70 : : !! charge and the XC core-correction,
71 : : !! if 1, treat both density change and XC core correction
72 : : !! if 2, like 0 but multiply gradient strain derivative term by 2.0 for GGA.
73 : : !! qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
74 : : !! rhor(nfft,nspden)=array for GS electron density in electrons/bohr**3.
75 : : !! rhor1(cplex*nfft,nspden)=array for electron density in electrons/bohr**3.
76 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
77 : : !! usepaw= 0 for non paw calculation; =1 for paw calculation
78 : : !! usexcnhat= -PAW only- flag controling use of compensation density in Vxc
79 : : !! xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc
80 : : !!
81 : : !! OUTPUT
82 : : !! vxc1(cplex*nfft,nspden)=change in exchange-correlation potential (including
83 : : !! core-correction, if applicable)
84 : : !!
85 : : !! SOURCE
86 : :
87 : 3109 : subroutine dfpt_mkvxcstr(cplex,idir,ipert,kxc,mpi_enreg,natom,nfft,ngfft,nhat,nhat1,&
88 : : & nkxc,non_magnetic_xc,nspden,n3xccc,option,qphon,&
89 : 3109 : & rhor,rhor1,rprimd,usepaw,usexcnhat,vxc1,xccc3d1)
90 : :
91 : : !Arguments ------------------------------------
92 : : !scalars
93 : : integer,intent(in) :: cplex,idir,ipert,n3xccc,natom,nfft,nkxc,nspden,option
94 : : integer,intent(in) :: usepaw,usexcnhat
95 : : logical,intent(in) :: non_magnetic_xc
96 : : type(MPI_type),intent(in) :: mpi_enreg
97 : : !arrays
98 : : integer,intent(in) :: ngfft(18)
99 : : real(dp),target,intent(in) :: nhat(nfft,nspden)
100 : : real(dp),target,intent(in) :: nhat1(cplex*nfft,nspden)
101 : : real(dp),intent(in) :: kxc(nfft,nkxc),qphon(3)
102 : : real(dp),target,intent(in) :: rhor(nfft,nspden),rhor1(cplex*nfft,nspden)
103 : : real(dp),intent(in) :: rprimd(3,3)
104 : : real(dp),intent(in) :: xccc3d1(cplex*n3xccc)
105 : : real(dp),intent(out) :: vxc1(cplex*nfft,nspden)
106 : :
107 : : !Local variables-------------------------------
108 : : !scalars
109 : : integer :: ii,ir,istr
110 : : real(dp) :: rho1_dn,rho1_up,spin_scale,str_scale
111 : : character(len=500) :: message
112 : : !arrays
113 : : real(dp) :: gprimd(3,3),tsec(2)
114 : 3109 : real(dp),allocatable :: rhor1tmp(:,:),rhowk1(:,:)
115 : 3109 : real(dp),pointer :: rhor_(:,:),rhor1_(:,:)
116 : :
117 : : ! *************************************************************************
118 : :
119 : 3109 : call timab(181,1,tsec)
120 : :
121 [ - + ]: 3109 : if(nspden/=1 .and. nspden/=2) then
122 : 0 : message = ' dfpt_mkvxc, Only for nspden==1 and 2.'
123 : 0 : ABI_BUG(message)
124 : : end if
125 : :
126 [ + + + + ]: 3109 : if (usepaw==1.and.usexcnhat==0) then
127 [ - + + - : 1032 : ABI_MALLOC(rhor_,(nfft,nspden))
- + + - +
- + - - +
- + ]
128 [ + - + + : 5064048 : rhor_(:,:)=rhor(:,:)-nhat(:,:)
+ + + + +
+ ]
129 : : else
130 : 2851 : rhor_ => rhor
131 : : end if
132 : :
133 [ + + + + : 3109 : if (usepaw==1.and.usexcnhat==0.and.option==1) then
+ + ]
134 [ - + + - : 48 : ABI_MALLOC(rhor1_,(nfft,nspden))
- + + - +
- + - - +
- + ]
135 [ + - + + : 186672 : rhor1_(:,:)=rhor1(:,:)-nhat1(:,:)
+ + + + +
+ ]
136 : : else
137 : 3097 : rhor1_ => rhor1
138 : : end if
139 : :
140 : : !Inhomogeneous term for diagonal strain
141 [ - + + - : 12436 : ABI_MALLOC(rhowk1,(nfft,nspden))
- + + - +
- + - - +
- + ]
142 [ + + ]: 3109 : if(option==0 .or. option==2) then
143 [ + + ]: 2941 : if(ipert==natom+3) then
144 [ + + + + ]: 10661975 : rhowk1(:,:)=-rhor_(:,:)
145 : : else
146 [ + + + + ]: 9771286 : rhowk1(:,:)=zero
147 : : end if
148 [ + - ]: 168 : else if(option==1) then
149 [ + + ]: 168 : if(ipert==natom+3) then
150 [ + + + + ]: 535475 : rhowk1(:,:)=rhor1_(:,:)-rhor_(:,:)
151 : : else
152 [ + + + + ]: 520016 : rhowk1(:,:)=rhor1_(:,:)
153 : : end if
154 : : end if
155 : :
156 [ - + ]: 3109 : if (non_magnetic_xc) then
157 [ # # # # ]: 0 : if(nspden==2) rhowk1(:,2)=rhowk1(:,1)*half
158 [ # # # # : 0 : if(nspden==4) rhowk1(:,2:4)=zero
# # ]
159 : : end if
160 : :
161 : : !Treat first LDA
162 [ + + ]: 3109 : if(nkxc==1.or.nkxc==3)then
163 : :
164 : : ! Case without non-linear core correction
165 [ + + ]: 2829 : if(n3xccc==0)then
166 : :
167 : : ! Non-spin-polarized
168 [ - + ]: 417 : if(nspden==1)then
169 [ + + ]: 2094701 : do ir=1,nfft
170 : 2094701 : vxc1(ir,1)=kxc(ir,1)*rhowk1(ir,1)
171 : : end do
172 : :
173 : : ! Spin-polarized
174 : : else
175 [ # # ]: 0 : do ir=1,nfft
176 : 0 : rho1_dn=rhowk1(ir,1)-rhowk1(ir,2)
177 : 0 : vxc1(ir,1)=kxc(ir,1)*rhowk1(ir,2)+kxc(ir,2)*rho1_dn
178 : 0 : vxc1(ir,2)=kxc(ir,2)*rhowk1(ir,2)+kxc(ir,3)*rho1_dn
179 : : end do
180 : : end if ! nspden==1
181 : :
182 : : ! Treat case with non-linear core correction
183 : : else
184 [ + + ]: 2412 : if(nspden==1)then
185 [ + + ]: 17308565 : do ir=1,nfft
186 : 17308565 : vxc1(ir,1)=kxc(ir,1)*(rhowk1(ir,1)+xccc3d1(ir))
187 : : end do
188 : : else
189 [ + + ]: 425809 : do ir=1,nfft
190 : 425736 : rho1_dn=rhowk1(ir,1)-rhowk1(ir,2) + xccc3d1(ir)*half
191 : 425736 : rho1_up=rhowk1(ir,2) + xccc3d1(ir)*half
192 : 425736 : vxc1(ir,1)=kxc(ir,1)*rho1_up+kxc(ir,2)*rho1_dn
193 : 425809 : vxc1(ir,2)=kxc(ir,2)*rho1_up+kxc(ir,3)*rho1_dn
194 : : end do
195 : : end if ! nspden==1
196 : :
197 : : end if ! n3xccc==0
198 : :
199 : : ! Treat GGA
200 [ + - ]: 280 : else if (nkxc==7.or.nkxc==19) then
201 : :
202 : : ! Generates gprimd and its strain derivative
203 : : ! Note that unlike the implicitly symmetric metric tensor strain
204 : : ! derivatives, we must explicltly symmetrize the strain derivative
205 : : ! here.
206 : 280 : call matr3inv(rprimd,gprimd)
207 : 280 : istr=idir + 3*(ipert-natom-3)
208 [ - + ]: 280 : if(istr<1 .or. istr>6)then
209 : : write(message, '(a,i10,a,a,a)' )&
210 : 0 : & 'Input dir gives istr=',istr,' not allowed.',ch10,&
211 : 0 : & 'Possible values are 1,2,3,4,5,6 only.'
212 : 0 : ABI_BUG(message)
213 : : end if
214 : :
215 : : ! Rescalling needed for use in dfpt_eltfrxc for elastic tensor (not internal strain tensor).
216 [ + + ]: 280 : str_scale=one;if(option==2) str_scale=two
217 : :
218 : : ! Transfer the data to spin-polarized storage
219 [ - + + - : 1120 : ABI_MALLOC(rhor1tmp,(cplex*nfft,nspden))
- + + - +
- + - - +
- + ]
220 [ + + ]: 280 : if(nspden==1)then
221 [ + + ]: 892445 : do ir=1,cplex*nfft
222 : 892445 : rhor1tmp(ir,1)=rhowk1(ir,1)
223 : : end do
224 : : else
225 [ + + ]: 169157 : do ir=1,cplex*nfft
226 : 169128 : rho1_dn=rhowk1(ir,1)-rhowk1(ir,2)
227 : 169128 : rhor1tmp(ir,1)=rhowk1(ir,2)
228 : 169157 : rhor1tmp(ir,2)=rho1_dn
229 : : end do
230 : : end if ! nspden==1
231 [ + + ]: 280 : if(n3xccc/=0)then
232 [ + + ]: 77 : spin_scale=one;if (nspden==2) spin_scale=half
233 [ + + ]: 183 : do ii=1,nspden
234 [ + + ]: 757383 : do ir=1,cplex*nfft
235 : 757306 : rhor1tmp(ir,ii)=rhor1tmp(ir,ii)+xccc3d1(ir)*spin_scale
236 : : end do
237 : : end do
238 : : end if
239 : :
240 : : call dfpt_mkvxcstrgga(cplex,gprimd,istr,kxc,mpi_enreg,nfft,ngfft,nkxc,&
241 : 280 : & nspden,qphon,rhor1tmp,str_scale,vxc1)
242 [ - + ]: 280 : ABI_FREE(rhor1tmp)
243 : :
244 : : else
245 : 0 : ABI_BUG('Invalid nkxc!')
246 : :
247 : : end if ! LDA or GGA
248 : :
249 [ + + + + ]: 3109 : if (usepaw==1.and.usexcnhat==0) then
250 [ - + ]: 258 : ABI_FREE(rhor_)
251 : : end if
252 [ + + + + : 3109 : if (usepaw==1.and.usexcnhat==0.and.option==1) then
+ + ]
253 [ - + ]: 12 : ABI_FREE(rhor1_)
254 : : end if
255 : :
256 [ - + ]: 3109 : ABI_FREE(rhowk1)
257 : :
258 : 3109 : call timab(181,2,tsec)
259 : :
260 : 3109 : end subroutine dfpt_mkvxcstr
261 : : !!***
262 : :
263 : : !!****f* ABINIT/dfpt_mkvxcstrgga
264 : : !! NAME
265 : : !! dfpt_mkvxcstrgga
266 : : !!
267 : : !! FUNCTION
268 : : !! Compute the first-order change of exchange-correlation potential
269 : : !! for the strain perturbation in case of GGA functionals
270 : : !! Use the exchange-correlation kernel.
271 : : !!
272 : : !! INPUTS
273 : : !! cplex= if 1, real space 1-order functions on FFT grid are REAL,
274 : : !! if 2, COMPLEX
275 : : !! gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1)
276 : : !! istr=index of the strain perturbation (1..6)
277 : : !! kxc(nfft,nkxc)=exchange and correlation kernel (see rhotoxc.f)
278 : : !! mpi_enreg=information about MPI parallelization
279 : : !! nfft=(effective) number of FFT grid points (for this processor)
280 : : !! ngfft(18)=contain all needed information about 3D FFT,
281 : : !! see ~abinit/doc/variables/vargs.htm#ngfft
282 : : !! nkxc=second dimension of the kxc array
283 : : !! nspden=number of spin-density components
284 : : !! qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
285 : : !! rhor1tmp(cplex*nfft,2)=array for first-order electron spin-density
286 : : !! in electrons/bohr**3 (second index corresponds to spin-up and spin-down)
287 : : !! str_scale=scaling factor for gradient operator strain-derivative (1. or 2.)
288 : : !!
289 : : !! OUTPUT
290 : : !! vxc1(cplex*nfft,nspden)=change in exchange-correlation potential
291 : : !!
292 : : !! NOTES
293 : : !! Closely related to dfpt_mkvxcgga.
294 : : !! Content of Kxc array:
295 : : !! ===== if GGA
296 : : !! if nspden==1:
297 : : !! kxc(:,1)= d2Exc/drho2
298 : : !! kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)|
299 : : !! kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho
300 : : !! kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| )
301 : : !! kxc(:,5)= gradx(rho)
302 : : !! kxc(:,6)= grady(rho)
303 : : !! kxc(:,7)= gradz(rho)
304 : : !! if nspden>=2:
305 : : !! kxc(:,1)= d2Exc/drho_up drho_up
306 : : !! kxc(:,2)= d2Exc/drho_up drho_dn
307 : : !! kxc(:,3)= d2Exc/drho_dn drho_dn
308 : : !! kxc(:,4)= 1/|grad(rho_up)| dEx/d|grad(rho_up)|
309 : : !! kxc(:,5)= 1/|grad(rho_dn)| dEx/d|grad(rho_dn)|
310 : : !! kxc(:,6)= 1/|grad(rho_up)| d2Ex/d|grad(rho_up)| drho_up
311 : : !! kxc(:,7)= 1/|grad(rho_dn)| d2Ex/d|grad(rho_dn)| drho_dn
312 : : !! kxc(:,8)= 1/|grad(rho_up)| * d/d|grad(rho_up)| ( 1/|grad(rho_up)| dEx/d|grad(rho_up)| )
313 : : !! kxc(:,9)= 1/|grad(rho_dn)| * d/d|grad(rho_dn)| ( 1/|grad(rho_dn)| dEx/d|grad(rho_dn)| )
314 : : !! kxc(:,10)=1/|grad(rho)| dEc/d|grad(rho)|
315 : : !! kxc(:,11)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_up
316 : : !! kxc(:,12)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_dn
317 : : !! kxc(:,13)=1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dEc/d|grad(rho)| )
318 : : !! kxc(:,14)=gradx(rho_up)
319 : : !! kxc(:,15)=gradx(rho_dn)
320 : : !! kxc(:,16)=grady(rho_up)
321 : : !! kxc(:,17)=grady(rho_dn)
322 : : !! kxc(:,18)=gradz(rho_up)
323 : : !! kxc(:,19)=gradz(rho_dn)
324 : : !!
325 : : !! SOURCE
326 : :
327 : 280 : subroutine dfpt_mkvxcstrgga(cplex,gprimd,istr,kxc,mpi_enreg,nfft,ngfft,&
328 : 280 : & nkxc,nspden,qphon,rhor1tmp,str_scale,vxc1)
329 : :
330 : : !Arguments ------------------------------------
331 : : !scalars
332 : : integer,intent(in) :: cplex,istr,nfft,nkxc,nspden
333 : : real(dp) :: str_scale
334 : : type(MPI_type),intent(in) :: mpi_enreg
335 : : !arrays
336 : : integer,intent(in) :: ngfft(18)
337 : : real(dp),intent(in) :: gprimd(3,3),kxc(nfft,nkxc)
338 : : real(dp),intent(in) :: qphon(3),rhor1tmp(cplex*nfft,2)
339 : : real(dp),intent(out) :: vxc1(cplex*nfft,nspden)
340 : : !Local variables-------------------------------
341 : : !scalars
342 : : integer :: ii,ir,ishift,ispden,mgga,ngrad,nspgrad
343 : : real(dp) :: coeff_grho,coeff_grho_corr,coeff_grho_dn,coeff_grho_up
344 : : real(dp) :: gradrho_gradrho1,gradrho_gradrho1_dn,gradrho_gradrho1_up
345 : : character(len=500) :: msg
346 : : !arrays
347 : : real(dp) :: r0(3),r0_dn(3),r0_up(3),r1(3),r1_dn(3),r1_up(3)
348 : 280 : real(dp),allocatable :: dnexcdn(:,:),rho1now(:,:,:),rhodgnow(:,:,:)
349 : :
350 : : ! *************************************************************************
351 : :
352 : : DBG_ENTER("COLL")
353 : :
354 [ - + ]: 280 : if (nkxc/=12*min(nspden,2)-5) then
355 : 0 : msg='Wrong nkxc value for GGA!'
356 : 0 : ABI_BUG(msg)
357 : : end if
358 [ - + ]: 280 : if (nspden>2) then
359 : 0 : msg='Not compatible with non-collinear magnetism!'
360 : 0 : ABI_ERROR(msg)
361 : : end if
362 : :
363 : : !metaGGA contributions are not taken into account here
364 : 280 : mgga=0
365 : :
366 : : !if you uncomment the following line, you will have to modify
367 : : !the original function call to pass in gmet and gsqcut
368 : : !call filterpot(cplex,gmet,gsqcut,nfft,ngfft,2,qphon,rhor1tmp)
369 : :
370 : : !Compute the gradients of the first-order density
371 : : !rho1now(:,:,1) contains the first-order density, and
372 : : !rho1now(:,:,2:4) contains the gradients of the first-order density
373 : 280 : ishift=0 ; ngrad=2
374 [ - + + - : 1400 : ABI_MALLOC(rho1now,(cplex*nfft,nspden,ngrad*ngrad))
- + + - +
- + - + -
- + - + ]
375 : 280 : call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden,qphon,rhor1tmp,rho1now)
376 : :
377 : : !Calculate the 1st-order contribution to grad(n) from the strain derivative
378 : : ! acting on the gradient operator acting on the GS charge density,
379 : : !Simply use the following formula:
380 : : ! (dGprim/ds_alpha_beta)(i,j) = -half.( delta_alpha,i Gprim(beta,j) + delta_beta,i Gprim(alpha,j) )
381 : : !To finally get:
382 : : ! (nabla)^(alpha,beta)_i[n] = -half ( delta_alpha,i nabla_beta[n] + delta_beta,i nabla_alpha[n] )
383 [ - + + - : 1400 : ABI_MALLOC(rhodgnow,(cplex*nfft,nspden,3))
- + + - +
- + - + -
- + + - ]
384 [ + + + + : 3693397 : rhodgnow(1:nfft,1:nspden,1:3)=zero
+ + ]
385 [ + + ]: 280 : if (nspden==1) then
386 [ + + + + ]: 197200 : if (istr==1) rhodgnow(1:nfft,1,1)=- kxc(1:nfft,5)
387 [ + + + + ]: 126700 : if (istr==2) rhodgnow(1:nfft,1,2)=- kxc(1:nfft,6)
388 [ + + + + ]: 125200 : if (istr==3) rhodgnow(1:nfft,1,3)=- kxc(1:nfft,7)
389 [ + + + + ]: 197200 : if (istr==4) rhodgnow(1:nfft,1,2)=-half*kxc(1:nfft,7)
390 [ + + + + ]: 197200 : if (istr==4) rhodgnow(1:nfft,1,3)=-half*kxc(1:nfft,6)
391 [ + + + + ]: 123700 : if (istr==5) rhodgnow(1:nfft,1,1)=-half*kxc(1:nfft,7)
392 [ + + + + ]: 123700 : if (istr==5) rhodgnow(1:nfft,1,3)=-half*kxc(1:nfft,5)
393 [ + + + + ]: 123700 : if (istr==6) rhodgnow(1:nfft,1,1)=-half*kxc(1:nfft,6)
394 [ + + + + ]: 123700 : if (istr==6) rhodgnow(1:nfft,1,2)=-half*kxc(1:nfft,5)
395 : : else
396 [ + + + + ]: 17525 : if (istr==1) rhodgnow(1:nfft,1,1)=- kxc(1:nfft,15)
397 [ + + + + ]: 17525 : if (istr==2) rhodgnow(1:nfft,1,2)=- kxc(1:nfft,17)
398 [ + + + + ]: 17525 : if (istr==3) rhodgnow(1:nfft,1,3)=- kxc(1:nfft,19)
399 [ + + + + ]: 81677 : if (istr==4) rhodgnow(1:nfft,1,2)=-half*kxc(1:nfft,19)
400 [ + + + + ]: 81677 : if (istr==4) rhodgnow(1:nfft,1,3)=-half*kxc(1:nfft,17)
401 [ + + + + ]: 17525 : if (istr==5) rhodgnow(1:nfft,1,1)=-half*kxc(1:nfft,19)
402 [ + + + + ]: 17525 : if (istr==5) rhodgnow(1:nfft,1,3)=-half*kxc(1:nfft,15)
403 [ + + + + ]: 17525 : if (istr==6) rhodgnow(1:nfft,1,1)=-half*kxc(1:nfft,17)
404 [ + + + + ]: 17525 : if (istr==6) rhodgnow(1:nfft,1,2)=-half*kxc(1:nfft,15)
405 [ + + + + ]: 17525 : if (istr==1) rhodgnow(1:nfft,2,1)=- (kxc(1:nfft,14)-kxc(1:nfft,15))
406 [ + + + + ]: 17525 : if (istr==2) rhodgnow(1:nfft,2,2)=- (kxc(1:nfft,16)-kxc(1:nfft,17))
407 [ + + + + ]: 17525 : if (istr==3) rhodgnow(1:nfft,2,3)=- (kxc(1:nfft,18)-kxc(1:nfft,19))
408 [ + + + + ]: 81677 : if (istr==4) rhodgnow(1:nfft,2,2)=-half*(kxc(1:nfft,18)-kxc(1:nfft,19))
409 [ + + + + ]: 81677 : if (istr==4) rhodgnow(1:nfft,2,3)=-half*(kxc(1:nfft,16)-kxc(1:nfft,17))
410 [ + + + + ]: 17525 : if (istr==5) rhodgnow(1:nfft,2,1)=-half*(kxc(1:nfft,18)-kxc(1:nfft,19))
411 [ + + + + ]: 17525 : if (istr==5) rhodgnow(1:nfft,2,3)=-half*(kxc(1:nfft,14)-kxc(1:nfft,15))
412 [ + + + + ]: 17525 : if (istr==6) rhodgnow(1:nfft,2,1)=-half*(kxc(1:nfft,16)-kxc(1:nfft,17))
413 [ + + + + ]: 17525 : if (istr==6) rhodgnow(1:nfft,2,2)=-half*(kxc(1:nfft,14)-kxc(1:nfft,15))
414 : : end if
415 : :
416 : : !Add to the gradients of the first-order density
417 [ + + ]: 1120 : do ii=1,3
418 [ + + ]: 2047 : do ispden=1,nspden
419 [ + + ]: 3692277 : rhodgnow(1:nfft,ispden,ii)=str_scale*rhodgnow(1:nfft,ispden,ii)
420 [ + + ]: 3693117 : rho1now(1:nfft,ispden,1+ii)=rho1now(1:nfft,ispden,1+ii)+rhodgnow(1:nfft,ispden,ii)
421 : : end do
422 : : end do
423 : :
424 : : !rho1now(:,:,1) contains the 1st-order density, and rho1now(:,:,2:4) contains the grads of the 1st-order density
425 : :
426 : : !Apply the XC kernel
427 [ + + ]: 280 : nspgrad=2; if (nspden==2) nspgrad=5
428 [ - + + - : 1120 : ABI_MALLOC(dnexcdn,(cplex*nfft,nspgrad))
- + + - +
- + - - +
- + ]
429 : :
430 : : !== Non polarized
431 [ + + ]: 280 : if (nspden==1) then
432 [ + + ]: 892445 : do ir=1,nfft
433 [ + + + + ]: 6245358 : r0(:)=kxc(ir,5:7) ; r1(:)=rho1now(ir,1,2:4)
434 [ + + ]: 3568776 : gradrho_gradrho1=dot_product(r0,r1)
435 : 892194 : dnexcdn(ir,1)=kxc(ir,1)*rho1now(ir,1,1) + kxc(ir,3)*gradrho_gradrho1
436 : 892194 : coeff_grho=kxc(ir,3)*rho1now(ir,1,1) + kxc(ir,4)*gradrho_gradrho1
437 : : ! Grad strain derivative contribution enters the following term with a
438 : : ! factor of two compared to above terms, so add it again.
439 [ + + ]: 3568776 : r1(:)=r1(:)+rhodgnow(ir,1,1:3)
440 : : ! Reuse the storage in rho1now
441 [ + + ]: 3569027 : rho1now(ir,1,2:4)=r1(:)*kxc(ir,2)+r0(:)*coeff_grho
442 : : end do
443 : :
444 : : !== Spin-polarized
445 : : else ! nspden==2
446 [ + + ]: 169157 : do ir=1,nfft
447 [ + + ]: 676512 : do ii=1,3 ! grad of spin-up ans spin_dn GS rho
448 : 676512 : r0_up(ii)=kxc(ir,13+2*ii);r0_dn(ii)=kxc(ir,12+2*ii)-kxc(ir,13+2*ii)
449 : : end do
450 [ + + ]: 676512 : r0(:)=r0_up(:)+r0_dn(:) ! grad of GS rho
451 [ + + ]: 676512 : r1_up(:)=rho1now(ir,1,2:4) ! grad of spin-up rho1
452 [ + + ]: 676512 : r1_dn(:)=rho1now(ir,2,2:4) ! grad of spin-down rho1
453 [ + + ]: 676512 : r1(:)=r1_up(:)+r1_dn(:) ! grad of GS rho1
454 [ + + ]: 676512 : gradrho_gradrho1_up=dot_product(r0_up,r1_up)
455 [ + + ]: 676512 : gradrho_gradrho1_dn=dot_product(r0_dn,r1_dn)
456 [ + + ]: 676512 : gradrho_gradrho1 =dot_product(r0,r1)
457 : :
458 : : dnexcdn(ir,1)=kxc(ir, 1)*rho1now(ir,1,1) &
459 : : & +kxc(ir, 2)*rho1now(ir,2,1) &
460 : : & +kxc(ir, 6)*gradrho_gradrho1_up &
461 : 169128 : & +kxc(ir,11)*gradrho_gradrho1
462 : : dnexcdn(ir,2)=kxc(ir, 3)*rho1now(ir,2,1) &
463 : : & +kxc(ir, 2)*rho1now(ir,1,1) &
464 : : & +kxc(ir, 7)*gradrho_gradrho1_dn &
465 : 169128 : & +kxc(ir,12)*gradrho_gradrho1
466 : : coeff_grho_corr=kxc(ir,11)*rho1now(ir,1,1) &
467 : : & +kxc(ir,12)*rho1now(ir,2,1) &
468 : 169128 : & +kxc(ir,13)*gradrho_gradrho1
469 : 169128 : coeff_grho_up=kxc(ir,6)*rho1now(ir,1,1)+kxc(ir,8)*gradrho_gradrho1_up
470 : 169128 : coeff_grho_dn=kxc(ir,7)*rho1now(ir,2,1)+kxc(ir,9)*gradrho_gradrho1_dn
471 : :
472 : : ! Grad strain derivative contribution enters the following term with a
473 : : ! factor of two compared to above terms, so add it again.
474 [ + + ]: 676512 : r1_up(:)=r1_up(:)+rhodgnow(ir,1,1:3)
475 [ + + ]: 676512 : r1_dn(:)=r1_dn(:)+rhodgnow(ir,2,1:3)
476 : :
477 : : ! Reuse the storage in rho1now
478 : : rho1now(ir,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1_up(:) &
479 : : & +kxc(ir,10) *r1_dn(:) &
480 : : & +coeff_grho_up *r0_up(:) &
481 [ + + ]: 676512 : & +coeff_grho_corr *r0(:)
482 : : rho1now(ir,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1_dn(:) &
483 : : & +kxc(ir,10) *r1_up(:) &
484 : : & +coeff_grho_dn *r0_dn(:) &
485 [ + + ]: 676541 : & +coeff_grho_corr *r0(:)
486 : : end do
487 : :
488 : : end if ! nspden
489 [ - + ]: 280 : ABI_FREE(rhodgnow)
490 : :
491 [ + + + + ]: 1231039 : vxc1(:,:)=zero
492 : : call xcpot(cplex,gprimd,ishift,mgga,mpi_enreg,nfft,ngfft,ngrad,nspden,&
493 : 280 : & nspgrad,qphon,depsxc=dnexcdn,rhonow=rho1now,vxc=vxc1)
494 : :
495 : : !if you uncomment the following line, you will have to modify
496 : : !the original function call to pass in gmet and gsqcut
497 : : !call filterpot(cplex,gmet,gsqcut,nfft,ngfft,nspden,qphon,vxc1)
498 : :
499 [ - + ]: 280 : ABI_FREE(dnexcdn)
500 [ - + ]: 280 : ABI_FREE(rho1now)
501 : :
502 : 280 : end subroutine dfpt_mkvxcstrgga
503 : : !!***
504 : :
505 : : end module m_dfpt_mkvxcstr
506 : : !!***
|
