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 : !!***
|
