Line data Source code
1 : !!****m* ABINIT/m_dfpt_mkvxc
2 : !! NAME
3 : !! m_dfpt_mkvxc
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2001-2026 ABINIT group (XG, DRH, FR, EB, SPr)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_dfpt_mkvxc
22 :
23 : use defs_basis
24 : use m_errors
25 : use m_abicore
26 : use m_xc_noncoll
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, xcpotdq
32 :
33 : implicit none
34 :
35 : private
36 : !!***
37 :
38 : public :: dfpt_mkvxc
39 : public :: dfpt_mkvxc_noncoll
40 : public :: dfpt_mkvxcggadq
41 : public :: dfpt_mkvxcgga_n0met
42 : public :: dfpt_mkvxcccdq
43 : !!***
44 :
45 : contains
46 : !!***
47 :
48 : !!****f* ABINIT/dfpt_mkvxc
49 : !! NAME
50 : !! dfpt_mkvxc
51 : !!
52 : !! FUNCTION
53 : !! Compute the first-order change of exchange-correlation potential
54 : !! due to atomic displacement: assemble the first-order
55 : !! density change with the frozen-core density change, then use
56 : !! the exchange-correlation kernel.
57 : !!
58 : !! INPUTS
59 : !! cplex= if 1, real space 1-order functions on FFT grid are REAL,
60 : !! if 2, COMPLEX
61 : !! ixc= choice of exchange-correlation scheme
62 : !! kxc(nfft,nkxc)=exchange and correlation kernel (see below)
63 : !! mpi_enreg=information about MPI parallelization
64 : !! nfft=(effective) number of FFT grid points (for this processor)
65 : !! ngfft(18)=contain all needed information about 3D FFT,
66 : !! see ~abinit/doc/variables/vargs.htm#ngfft
67 : !! nhat1(cplex*nfft,2nspden*nhat1dim)= -PAW only- 1st-order compensation density
68 : !! nhat1dim= -PAW only- 1 if nhat1 array is used ; 0 otherwise
69 : !! nhat1gr(cplex*nfft,nspden,3*nhat1grdim)= -PAW only- gradients of 1st-order compensation density
70 : !! nhat1grdim= -PAW only- 1 if nhat1gr array is used ; 0 otherwise
71 : !! nkxc=second dimension of the kxc array
72 : !! non_magnetic_xc= if true, handle density/potential as non-magnetic (even if it is)
73 : !! nspden=number of spin-density components
74 : !! n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used, otherwise, nfft
75 : !! option=if 0, work only with the XC core-correction,
76 : !! if 1, treat both density change and XC core correction
77 : !! if 2, treat only density change
78 : !! qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
79 : !! rhor1(cplex*nfft,nspden)=array for electron density in electrons/bohr**3.
80 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
81 : !! usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
82 : !! xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc
83 : !!
84 : !! OUTPUT
85 : !! vxc1(cplex*nfft,nspden)=change in exchange-correlation potential (including
86 : !! core-correction, if applicable)
87 : !!
88 : !! NOTES
89 : !! Content of Kxc array:
90 : !! ===== if LDA
91 : !! if nspden==1: kxc(:,1)= d2Exc/drho2
92 : !! (kxc(:,2)= d2Exc/drho_up drho_dn)
93 : !! if nspden>=2: kxc(:,1)= d2Exc/drho_up drho_up
94 : !! kxc(:,2)= d2Exc/drho_up drho_dn
95 : !! kxc(:,3)= d2Exc/drho_dn drho_dn
96 : !! ===== if GGA (or mGGA)
97 : !! if nspden==1:
98 : !! kxc(:,1)= d2Exc/drho2
99 : !! kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)|
100 : !! kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho
101 : !! kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| )
102 : !! kxc(:,5)= gradx(rho)
103 : !! kxc(:,6)= grady(rho)
104 : !! kxc(:,7)= gradz(rho)
105 : !! if nspden>=2:
106 : !! kxc(:,1)= d2Exc/drho_up drho_up
107 : !! kxc(:,2)= d2Exc/drho_up drho_dn
108 : !! kxc(:,3)= d2Exc/drho_dn drho_dn
109 : !! kxc(:,4)= 1/|grad(rho_up)| dEx/d|grad(rho_up)|
110 : !! kxc(:,5)= 1/|grad(rho_dn)| dEx/d|grad(rho_dn)|
111 : !! kxc(:,6)= 1/|grad(rho_up)| d2Ex/d|grad(rho_up)| drho_up
112 : !! kxc(:,7)= 1/|grad(rho_dn)| d2Ex/d|grad(rho_dn)| drho_dn
113 : !! kxc(:,8)= 1/|grad(rho_up)| * d/d|grad(rho_up)| ( 1/|grad(rho_up)| dEx/d|grad(rho_up)| )
114 : !! kxc(:,9)= 1/|grad(rho_dn)| * d/d|grad(rho_dn)| ( 1/|grad(rho_dn)| dEx/d|grad(rho_dn)| )
115 : !! kxc(:,10)=1/|grad(rho)| dEc/d|grad(rho)|
116 : !! kxc(:,11)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_up
117 : !! kxc(:,12)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_dn
118 : !! kxc(:,13)=1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dEc/d|grad(rho)| )
119 : !! kxc(:,14)=gradx(rho_up)
120 : !! kxc(:,15)=gradx(rho_dn)
121 : !! kxc(:,16)=grady(rho_up)
122 : !! kxc(:,17)=grady(rho_dn)
123 : !! kxc(:,18)=gradz(rho_up)
124 : !! kxc(:,19)=gradz(rho_dn)
125 : !! Note about mGGA: 2nd derivatives involving Tau or Laplacian are not taken into account (yet)
126 : !!
127 : !! SOURCE
128 :
129 55756 : subroutine dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1,nhat1dim,nhat1gr,nhat1grdim,&
130 55756 : & nkxc,non_magnetic_xc,nspden,n3xccc,option,qphon,rhor1,rprimd,usexcnhat,vxc1,xccc3d1)
131 :
132 : !Arguments ------------------------------------
133 : !scalars
134 : integer,intent(in) :: cplex,ixc,n3xccc,nfft,nhat1dim,nhat1grdim
135 : integer,intent(in) :: nkxc,nspden,option,usexcnhat
136 : logical,intent(in) :: non_magnetic_xc
137 : type(MPI_type),intent(in) :: mpi_enreg
138 : !arrays
139 : integer,intent(in) :: ngfft(18)
140 : real(dp),intent(in),target :: nhat1(cplex*nfft,nspden*nhat1dim)
141 : real(dp),intent(in),target :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim)
142 : real(dp),intent(in) :: kxc(nfft,nkxc),qphon(3)
143 : real(dp),intent(in),target :: rhor1(cplex*nfft,nspden)
144 : real(dp),intent(in) :: rprimd(3,3),xccc3d1(cplex*n3xccc)
145 : real(dp),intent(out) :: vxc1(cplex*nfft,nspden)
146 :
147 : !Local variables-------------------------------
148 : !scalars
149 : integer :: ii,ir,ispden,nhat1dim_,nhat1rgdim_
150 : real(dp) :: rho1_dn,rho1_up,rho1im_dn,rho1im_up,rho1re_dn,rho1re_up
151 : real(dp) :: spin_scale
152 : !arrays
153 : real(dp) :: gprimd(3,3),tsec(2)
154 55756 : real(dp), contiguous, pointer :: nhat1_(:,:),nhat1gr_(:,:,:),rhor1_(:,:)
155 : ! *************************************************************************
156 :
157 : DBG_ENTER("COLL")
158 :
159 55756 : call timab(181,1,tsec)
160 :
161 55756 : if(nspden/=1 .and. nspden/=2) then
162 0 : ABI_BUG('For nspden==4 please use dfpt_mkvxc_noncoll!')
163 : end if
164 :
165 : !Special case: no XC applied
166 55756 : if (ixc==0.or.nkxc==0) then
167 962 : ABI_WARNING('Note that no xc is applied (ixc=0)')
168 11332295 : vxc1=zero
169 962 : return
170 : end if
171 :
172 : !Treat first LDA
173 : if(nkxc==1.or.nkxc==3)then
174 :
175 : ! PAW: eventually substract compensation density
176 49413 : if (option/=0) then
177 25187 : if ((usexcnhat==0.and.nhat1dim==1).or.(non_magnetic_xc)) then
178 4856 : ABI_MALLOC(rhor1_,(cplex*nfft,nspden))
179 1214 : if (usexcnhat==0.and.nhat1dim==1) then
180 63056689 : rhor1_(:,:)=rhor1(:,:)-nhat1(:,:)
181 : else
182 1576591 : rhor1_(:,:)=rhor1(:,:)
183 : end if
184 1214 : if (non_magnetic_xc) then
185 394142 : if(nspden==2) rhor1_(:,2)=rhor1_(:,1)*half
186 23 : if(nspden==4) rhor1_(:,2:4)=zero
187 : end if
188 : else
189 23973 : rhor1_ => rhor1
190 : end if
191 : end if
192 :
193 : ! Case without non-linear core correction
194 49413 : if(n3xccc==0 .or. option==2)then
195 :
196 28505 : if(option==0)then ! No straight XC to compute
197 :
198 49015236 : vxc1(:,:)=zero
199 :
200 : else ! XC, without non-linear XC correction
201 :
202 : ! Non-spin-polarized
203 23464 : if(nspden==1)then
204 21291 : if(cplex==1)then
205 414164092 : do ir=1,nfft
206 414164092 : vxc1(ir,1)=kxc(ir,1)*rhor1_(ir,1)
207 : end do
208 : else
209 32743204 : do ir=1,nfft
210 32738092 : vxc1(2*ir-1,1)=kxc(ir,1)*rhor1_(2*ir-1,1)
211 32743204 : vxc1(2*ir ,1)=kxc(ir,1)*rhor1_(2*ir ,1)
212 : end do
213 : end if ! cplex==1
214 :
215 : ! Spin-polarized
216 : else
217 2173 : if(cplex==1)then
218 33048018 : do ir=1,nfft
219 33046010 : rho1_dn=rhor1_(ir,1)-rhor1_(ir,2)
220 33046010 : vxc1(ir,1)=kxc(ir,1)*rhor1_(ir,2)+kxc(ir,2)*rho1_dn
221 33048018 : vxc1(ir,2)=kxc(ir,2)*rhor1_(ir,2)+kxc(ir,3)*rho1_dn
222 : end do
223 : else
224 923917 : do ir=1,nfft
225 923752 : rho1re_dn=rhor1_(2*ir-1,1)-rhor1_(2*ir-1,2)
226 923752 : rho1im_dn=rhor1_(2*ir ,1)-rhor1_(2*ir ,2)
227 923752 : vxc1(2*ir-1,1)=kxc(ir,1)*rhor1_(2*ir-1,2)+kxc(ir,2)*rho1re_dn
228 923752 : vxc1(2*ir ,1)=kxc(ir,1)*rhor1_(2*ir ,2)+kxc(ir,2)*rho1im_dn
229 923752 : vxc1(2*ir-1,2)=kxc(ir,2)*rhor1_(2*ir-1,2)+kxc(ir,3)*rho1re_dn
230 923917 : vxc1(2*ir ,2)=kxc(ir,2)*rhor1_(2*ir ,2)+kxc(ir,3)*rho1im_dn
231 : end do
232 : end if ! cplex==1
233 : end if ! nspden==1
234 :
235 : end if ! option==0
236 :
237 : ! Treat case with non-linear core correction
238 : else
239 :
240 20908 : if(option==0)then
241 :
242 19185 : if(nspden==1)then
243 17040 : if(cplex==1)then
244 127401478 : do ir=1,nfft
245 127401478 : vxc1(ir,1)=kxc(ir,1)*xccc3d1(ir)
246 : end do
247 : else
248 63774383 : do ir=1,nfft
249 63766943 : vxc1(2*ir-1,1)=kxc(ir,1)*xccc3d1(2*ir-1)
250 63774383 : vxc1(2*ir ,1)=kxc(ir,1)*xccc3d1(2*ir )
251 : end do
252 : end if ! cplex==1
253 : else
254 2145 : if(cplex==1)then
255 11507656 : do ir=1,nfft
256 11506072 : vxc1(ir,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(ir)*half
257 11507656 : vxc1(ir,2)=(kxc(ir,2)+kxc(ir,3))*xccc3d1(ir)*half
258 : end do
259 : else
260 2229223 : do ir=1,nfft
261 2228662 : vxc1(2*ir-1,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(2*ir-1)*half
262 2228662 : vxc1(2*ir ,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(2*ir )*half
263 2228662 : vxc1(2*ir-1,2)=(kxc(ir,2)+kxc(ir,3))*xccc3d1(2*ir-1)*half
264 2229223 : vxc1(2*ir ,2)=(kxc(ir,2)+kxc(ir,3))*xccc3d1(2*ir )*half
265 : end do
266 : end if ! cplex==1
267 : end if ! nspden==1
268 :
269 : else ! option/=0
270 :
271 1723 : if(nspden==1)then
272 1222 : if(cplex==1)then
273 7442937 : do ir=1,nfft
274 7442937 : vxc1(ir,1)=kxc(ir,1)*(rhor1_(ir,1)+xccc3d1(ir))
275 : end do
276 : else
277 3280034 : do ir=1,nfft
278 3279618 : vxc1(2*ir-1,1)=kxc(ir,1)*(rhor1_(2*ir-1,1)+xccc3d1(2*ir-1))
279 3280034 : vxc1(2*ir ,1)=kxc(ir,1)*(rhor1_(2*ir ,1)+xccc3d1(2*ir ))
280 : end do
281 : end if ! cplex==1
282 : else
283 501 : if(cplex==1)then
284 2858069 : do ir=1,nfft
285 2857587 : rho1_dn=rhor1_(ir,1)-rhor1_(ir,2) + xccc3d1(ir)*half
286 2857587 : rho1_up=rhor1_(ir,2) + xccc3d1(ir)*half
287 2857587 : vxc1(ir,1)=kxc(ir,1)*rho1_up+kxc(ir,2)*rho1_dn
288 2858069 : vxc1(ir,2)=kxc(ir,2)*rho1_up+kxc(ir,3)*rho1_dn
289 : end do
290 : else
291 93809 : do ir=1,nfft
292 93790 : rho1re_dn=rhor1_(2*ir-1,1)-rhor1_(2*ir-1,2) + xccc3d1(2*ir-1)*half
293 93790 : rho1im_dn=rhor1_(2*ir ,1)-rhor1_(2*ir ,2) + xccc3d1(2*ir )*half
294 93790 : rho1re_up=rhor1_(2*ir-1,2) + xccc3d1(2*ir-1)*half
295 93790 : rho1im_up=rhor1_(2*ir ,2) + xccc3d1(2*ir )*half
296 93790 : vxc1(2*ir-1,1)=kxc(ir,1)*rho1re_up+kxc(ir,2)*rho1re_dn
297 93790 : vxc1(2*ir ,1)=kxc(ir,1)*rho1im_up+kxc(ir,2)*rho1im_dn
298 93790 : vxc1(2*ir-1,2)=kxc(ir,2)*rho1re_up+kxc(ir,3)*rho1re_dn
299 93809 : vxc1(2*ir ,2)=kxc(ir,2)*rho1im_up+kxc(ir,3)*rho1im_dn
300 : end do
301 : end if ! cplex==1
302 : end if ! nspden==1
303 :
304 : end if ! option==0
305 :
306 : end if ! n3xccc==0
307 :
308 49413 : if (option/=0.and.((usexcnhat==0.and.nhat1dim==1).or.(non_magnetic_xc))) then
309 1214 : ABI_FREE(rhor1_)
310 : end if
311 :
312 : ! Treat GGA
313 : else if (nkxc==7.or.nkxc==19) then
314 :
315 : ! Transfer the data to spin-polarized storage
316 :
317 : ! Treat the density change
318 21524 : ABI_MALLOC(rhor1_,(cplex*nfft,nspden))
319 5381 : if (option==1 .or. option==2) then
320 2578 : if (nspden==1) then
321 74101236 : do ir=1,cplex*nfft
322 74101236 : rhor1_(ir,1)=rhor1(ir,1)
323 : end do
324 : else
325 139 : if(non_magnetic_xc) then
326 0 : do ir=1,cplex*nfft
327 0 : rho1_dn=rhor1(ir,1)*half
328 0 : rhor1_(ir,1)=rho1_dn
329 0 : rhor1_(ir,2)=rho1_dn
330 : end do
331 : else
332 1223332 : do ir=1,cplex*nfft
333 1223193 : rho1_dn=rhor1(ir,1)-rhor1(ir,2)
334 1223193 : rhor1_(ir,1)=rhor1(ir,2)
335 1223332 : rhor1_(ir,2)=rho1_dn
336 : end do
337 : end if
338 : end if
339 : else
340 5690 : do ispden=1,nspden
341 117382418 : do ir=1,cplex*nfft
342 117379615 : rhor1_(ir,ispden)=zero
343 : end do
344 : end do
345 : end if
346 :
347 5381 : if( (option==0 .or. option==1) .and. n3xccc/=0)then
348 2573 : spin_scale=one;if (nspden==2) spin_scale=half
349 5213 : do ispden=1,nspden
350 96862533 : do ir=1,cplex*nfft
351 96859960 : rhor1_(ir,ispden)=rhor1_(ir,ispden)+xccc3d1(ir)*spin_scale
352 : end do
353 : end do
354 : end if
355 :
356 : ! PAW: treat also compensation density (and gradients)
357 5381 : nhat1dim_=nhat1dim ; nhat1rgdim_=nhat1grdim
358 5381 : if (option/=0.and.nhat1dim==1.and.nspden==2) then
359 165 : ABI_MALLOC(nhat1_,(cplex*nfft,nspden))
360 55 : if (non_magnetic_xc) then
361 0 : do ir=1,cplex*nfft
362 0 : rho1_dn=nhat1(ir,1)*half
363 0 : nhat1_(ir,1:2)=rho1_dn
364 : end do
365 : else
366 864439 : do ir=1,cplex*nfft
367 864384 : rho1_dn=nhat1(ir,1)-nhat1(ir,2)
368 864384 : nhat1_(ir,1)=nhat1(ir,2)
369 864439 : nhat1_(ir,2)=rho1_dn
370 : end do
371 : end if
372 5326 : else if (option==0) then
373 2803 : ABI_MALLOC(nhat1_,(0,0))
374 2803 : nhat1dim_=0
375 : else
376 2523 : nhat1_ => nhat1
377 : end if
378 5381 : if (option/=0.and.nhat1grdim==1.and.nspden==2) then
379 40 : ABI_MALLOC(nhat1gr_,(cplex*nfft,nspden,3))
380 8 : if (non_magnetic_xc) then
381 0 : do ii=1,3
382 0 : do ir=1,cplex*nfft
383 0 : rho1_dn=nhat1(ir,1)*half
384 0 : nhat1gr_(ir,1:2,ii)=rho1_dn
385 : end do
386 : end do
387 : else
388 32 : do ii=1,3
389 279968 : do ir=1,cplex*nfft
390 279936 : rho1_dn=nhat1gr(ir,1,ii)-nhat1gr(ir,2,ii)
391 279936 : nhat1gr_(ir,1,ii)=nhat1gr(ir,2,ii)
392 279960 : nhat1gr_(ir,2,ii)=rho1_dn
393 : end do
394 : end do
395 : end if
396 5373 : else if (option==0) then
397 2803 : ABI_MALLOC(nhat1gr_,(0,0,0))
398 2803 : nhat1rgdim_=0
399 : else
400 2570 : nhat1gr_ => nhat1gr
401 : end if
402 :
403 5381 : call matr3inv(rprimd,gprimd)
404 :
405 : call dfpt_mkvxcgga(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat1_,nhat1dim_,&
406 5381 : & nhat1gr_,nhat1rgdim_,nkxc,nspden,qphon,rhor1_,usexcnhat,vxc1)
407 :
408 5381 : ABI_FREE(rhor1_)
409 5381 : if ((option==0).or.(nhat1dim==1.and.nspden==2)) then
410 2858 : ABI_FREE(nhat1_)
411 : end if
412 5381 : if ((option==0).or.(nhat1grdim==1.and.nspden==2)) then
413 2811 : ABI_FREE(nhat1gr_)
414 : end if
415 :
416 : else
417 0 : ABI_BUG('Invalid nkxc!')
418 :
419 : end if ! LDA or GGA
420 :
421 54794 : call timab(181,2,tsec)
422 :
423 : DBG_EXIT("COLL")
424 :
425 55756 : end subroutine dfpt_mkvxc
426 : !!***
427 :
428 : !!****f* ABINIT/dfpt_mkvxcgga
429 : !! NAME
430 : !! dfpt_mkvxcgga
431 : !!
432 : !! FUNCTION
433 : !! Compute the first-order change of exchange-correlation potential
434 : !! in case of GGA functionals
435 : !! Use the exchange-correlation kernel.
436 : !!
437 : !! INPUTS
438 : !! cplex= if 1, real space 1-order functions on FFT grid are REAL,
439 : !! if 2, COMPLEX
440 : !! gmet(3,3)=metrix tensor in G space in Bohr**-2.
441 : !! gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1)
442 : !! gsqcut=cutoff value on G**2 for sphere inside fft box.
443 : !! kxc(nfft,nkxc)=exchange and correlation kernel (see below)
444 : !! mpi_enreg=information about MPI parallelization
445 : !! nfft=(effective) number of FFT grid points (for this processor)
446 : !! ngfft(18)=contain all needed information about 3D FFT
447 : !! nhat1(cplex*nfft,2*nhat1dim)= -PAW only- 1st-order compensation density
448 : !! nhat1dim= -PAW only- 1 if nhat1 array is used ; 0 otherwise
449 : !! nhat1gr(cplex*nfft,2,3*nhat1grdim)= -PAW only- gradients of 1st-order compensation density
450 : !! nhat1grdim= -PAW only- 1 if nhat1gr array is used ; 0 otherwise
451 : !! nkxc=second dimension of the kxc array
452 : !! nspden=number of spin-density components
453 : !! qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
454 : !! rhor1tmp(cplex*nfft,2)=array for first-order electron spin-density
455 : !! in electrons/bohr**3 (second index corresponds to spin-up and spin-down)
456 : !! usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
457 : !!
458 : !! OUTPUT
459 : !! vxc1(cplex*nfft,nspden)=change in exchange-correlation potential
460 : !!
461 : !! NOTES
462 : !! For the time being, a rather crude coding, to be optimized ...
463 : !! Content of Kxc array:
464 : !! ===== if GGA
465 : !! if nspden==1:
466 : !! kxc(:,1)= d2Exc/drho2
467 : !! kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)|
468 : !! kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho
469 : !! kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| )
470 : !! kxc(:,5)= gradx(rho)
471 : !! kxc(:,6)= grady(rho)
472 : !! kxc(:,7)= gradz(rho)
473 : !! if nspden>=2:
474 : !! kxc(:,1)= d2Exc/drho_up drho_up
475 : !! kxc(:,2)= d2Exc/drho_up drho_dn
476 : !! kxc(:,3)= d2Exc/drho_dn drho_dn
477 : !! kxc(:,4)= 1/|grad(rho_up)| dEx/d|grad(rho_up)|
478 : !! kxc(:,5)= 1/|grad(rho_dn)| dEx/d|grad(rho_dn)|
479 : !! kxc(:,6)= 1/|grad(rho_up)| d2Ex/d|grad(rho_up)| drho_up
480 : !! kxc(:,7)= 1/|grad(rho_dn)| d2Ex/d|grad(rho_dn)| drho_dn
481 : !! kxc(:,8)= 1/|grad(rho_up)| * d/d|grad(rho_up)| ( 1/|grad(rho_up)| dEx/d|grad(rho_up)| )
482 : !! kxc(:,9)= 1/|grad(rho_dn)| * d/d|grad(rho_dn)| ( 1/|grad(rho_dn)| dEx/d|grad(rho_dn)| )
483 : !! kxc(:,10)=1/|grad(rho)| dEc/d|grad(rho)|
484 : !! kxc(:,11)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_up
485 : !! kxc(:,12)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_dn
486 : !! kxc(:,13)=1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dEc/d|grad(rho)| )
487 : !! kxc(:,14)=gradx(rho_up)
488 : !! kxc(:,15)=gradx(rho_dn)
489 : !! kxc(:,16)=grady(rho_up)
490 : !! kxc(:,17)=grady(rho_dn)
491 : !! kxc(:,18)=gradz(rho_up)
492 : !! kxc(:,19)=gradz(rho_dn)
493 : !!
494 : !! SOURCE
495 :
496 5381 : subroutine dfpt_mkvxcgga(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,&
497 5381 : & nhat1,nhat1dim,nhat1gr,nhat1grdim,nkxc,&
498 5381 : & nspden,qphon,rhor1,usexcnhat,vxc1)
499 :
500 : !Arguments ------------------------------------
501 : !scalars
502 : integer,intent(in) :: cplex,nfft,nhat1dim,nhat1grdim,nkxc,nspden,usexcnhat
503 : type(MPI_type),intent(in) :: mpi_enreg
504 : !arrays
505 : integer,intent(in) :: ngfft(18)
506 : real(dp),intent(in) :: gprimd(3,3)
507 : real(dp),intent(in) :: kxc(nfft,nkxc)
508 : real(dp),intent(in) :: nhat1(cplex*nfft,nspden*nhat1dim)
509 : real(dp),intent(in) :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim)
510 : real(dp),intent(in) :: qphon(3)
511 : real(dp),intent(in),target :: rhor1(cplex*nfft,nspden)
512 : real(dp),intent(out) :: vxc1(cplex*nfft,nspden)
513 :
514 : !Local variables-------------------------------
515 : !scalars
516 : integer :: ii,ir,ishift,ngrad,nspgrad,use_laplacian
517 : logical :: test_nhat
518 : real(dp) :: coeff_grho,coeff_grho_corr,coeff_grho_dn,coeff_grho_up
519 : real(dp) :: coeffim_grho,coeffim_grho_corr,coeffim_grho_dn,coeffim_grho_up
520 : real(dp) :: gradrho_gradrho1,gradrho_gradrho1_dn,gradrho_gradrho1_up
521 : real(dp) :: gradrho_gradrho1im,gradrho_gradrho1im_dn,gradrho_gradrho1im_up
522 : character(len=500) :: msg
523 : !arrays
524 : real(dp) :: r0(3),r0_dn(3),r0_up(3),r1(3),r1_dn(3),r1_up(3)
525 : real(dp) :: r1im(3),r1im_dn(3),r1im_up(3)
526 5381 : real(dp),allocatable :: dnexcdn(:,:),rho1now(:,:,:)
527 5381 : real(dp),contiguous, pointer :: rhor1_ptr(:,:)
528 : ! *************************************************************************
529 :
530 : DBG_ENTER("COLL")
531 :
532 5381 : if (nkxc/=12*min(nspden,2)-5) then
533 0 : msg='Wrong nkxc value for GGA!'
534 0 : ABI_BUG(msg)
535 : end if
536 :
537 : !metaGGA contributions are not taken into account here
538 5381 : use_laplacian=0
539 :
540 : !PAW: substract 1st-order compensation density from 1st-order density
541 5381 : test_nhat=((nhat1dim==1).and.(usexcnhat==0.or.nhat1grdim==1))
542 : if (test_nhat) then
543 2392 : ABI_MALLOC(rhor1_ptr,(cplex*nfft,nspden))
544 72080826 : rhor1_ptr(:,:)=rhor1(:,:)-nhat1(:,:)
545 : else
546 4783 : rhor1_ptr => rhor1
547 : end if
548 :
549 : !call filterpot(paral_kgb,cplex,gmet,gsqcut,nfft,ngfft,2,qphon,rhor1_ptr)
550 :
551 : !Compute the gradients of the first-order density
552 : !rho1now(:,:,1) contains the first-order density, and
553 : !rho1now(:,:,2:4) contains the gradients of the first-order density
554 5381 : ishift=0 ; ngrad=2
555 26905 : ABI_MALLOC(rho1now,(cplex*nfft,nspden,ngrad*ngrad))
556 5381 : call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden,qphon,rhor1_ptr,rho1now)
557 :
558 : !PAW: add "exact" gradients of compensation density
559 5381 : if (test_nhat.and.usexcnhat==1) then
560 186648 : rho1now(:,1:nspden,1)=rho1now(:,1:nspden,1)+nhat1(:,1:nspden)
561 : end if
562 5381 : if (nhat1grdim==1) then
563 32 : do ii=1,3
564 559952 : rho1now(:,1:nspden,ii+1)=rho1now(:,1:nspden,ii+1)+nhat1gr(:,1:nspden,ii)
565 : end do
566 : end if
567 5381 : if (test_nhat) then
568 598 : ABI_FREE(rhor1_ptr)
569 : end if
570 :
571 : !Apply the XC kernel
572 5381 : nspgrad=2; if (nspden==2) nspgrad=5
573 21524 : ABI_MALLOC(dnexcdn,(cplex*nfft,nspgrad))
574 :
575 5381 : if (cplex==1) then ! Treat real case first
576 3984 : if (nspden==1) then
577 161986334 : do ir=1,nfft
578 1133877731 : r0(:)=kxc(ir,5:7) ; r1(:)=rho1now(ir,1,2:4)
579 647930132 : gradrho_gradrho1=dot_product(r0,r1)
580 161982533 : dnexcdn(ir,1)=kxc(ir,1)*rho1now(ir,1,1) + kxc(ir,3)*gradrho_gradrho1
581 161982533 : coeff_grho=kxc(ir,3)*rho1now(ir,1,1) + kxc(ir,4)*gradrho_gradrho1
582 : ! Reuse the storage in rho1now
583 647933933 : rho1now(ir,1,2:4)=r1(:)*kxc(ir,2)+r0(:)*coeff_grho
584 : end do
585 : else
586 1703232 : do ir=1,nfft
587 6812196 : do ii=1,3 ! grad of spin-up ans spin_dn GS rho
588 6812196 : r0_up(ii)=kxc(ir,13+2*ii);r0_dn(ii)=kxc(ir,12+2*ii)-kxc(ir,13+2*ii)
589 : end do
590 6812196 : r0(:)=r0_up(:)+r0_dn(:) ! grad of GS rho
591 6812196 : r1_up(:)=rho1now(ir,1,2:4) ! grad of spin-up rho1
592 6812196 : r1_dn(:)=rho1now(ir,2,2:4) ! grad of spin-down rho1
593 6812196 : r1(:)=r1_up(:)+r1_dn(:) ! grad of GS rho1
594 6812196 : gradrho_gradrho1_up=dot_product(r0_up,r1_up)
595 6812196 : gradrho_gradrho1_dn=dot_product(r0_dn,r1_dn)
596 6812196 : gradrho_gradrho1 =dot_product(r0,r1)
597 : dnexcdn(ir,1)=kxc(ir, 1)*rho1now(ir,1,1) &
598 : & +kxc(ir, 2)*rho1now(ir,2,1) &
599 : & +kxc(ir, 6)*gradrho_gradrho1_up &
600 1703049 : & +kxc(ir,11)*gradrho_gradrho1
601 : dnexcdn(ir,2)=kxc(ir, 3)*rho1now(ir,2,1) &
602 : & +kxc(ir, 2)*rho1now(ir,1,1) &
603 : & +kxc(ir, 7)*gradrho_gradrho1_dn &
604 1703049 : & +kxc(ir,12)*gradrho_gradrho1
605 : coeff_grho_corr=kxc(ir,11)*rho1now(ir,1,1) &
606 : & +kxc(ir,12)*rho1now(ir,2,1) &
607 1703049 : & +kxc(ir,13)*gradrho_gradrho1
608 1703049 : coeff_grho_up=kxc(ir,6)*rho1now(ir,1,1)+kxc(ir,8)*gradrho_gradrho1_up
609 1703049 : coeff_grho_dn=kxc(ir,7)*rho1now(ir,2,1)+kxc(ir,9)*gradrho_gradrho1_dn
610 : ! Reuse the storage in rho1now
611 : rho1now(ir,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1_up(:) &
612 : & +kxc(ir,10) *r1_dn(:) &
613 : & +coeff_grho_up *r0_up(:) &
614 6812196 : & +coeff_grho_corr *r0(:)
615 : rho1now(ir,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1_dn(:) &
616 : & +kxc(ir,10) *r1_up(:) &
617 : & +coeff_grho_dn *r0_dn(:) &
618 6812379 : & +coeff_grho_corr *r0(:)
619 : end do
620 : end if ! nspden
621 :
622 : else ! if cplex==2
623 1397 : if (nspden==1) then
624 13801437 : do ir=1,nfft
625 55200320 : r0(:)=kxc(ir,5:7)
626 55200320 : r1(:) =rho1now(2*ir-1,1,2:4)
627 55200320 : r1im(:)=rho1now(2*ir ,1,2:4)
628 55200320 : gradrho_gradrho1 =dot_product(r0,r1)
629 55200320 : gradrho_gradrho1im=dot_product(r0,r1im)
630 13800080 : dnexcdn(2*ir-1,1)=kxc(ir,1)*rho1now(2*ir-1,1,1) + kxc(ir,3)*gradrho_gradrho1
631 13800080 : dnexcdn(2*ir ,1)=kxc(ir,1)*rho1now(2*ir ,1,1) + kxc(ir,3)*gradrho_gradrho1im
632 13800080 : coeff_grho =kxc(ir,3)*rho1now(2*ir-1,1,1) + kxc(ir,4)*gradrho_gradrho1
633 13800080 : coeffim_grho=kxc(ir,3)*rho1now(2*ir ,1,1) + kxc(ir,4)*gradrho_gradrho1im
634 : ! Reuse the storage in rho1now
635 55200320 : rho1now(2*ir-1,1,2:4)=r1(:) *kxc(ir,2)+r0(:)*coeff_grho
636 55201677 : rho1now(2*ir ,1,2:4)=r1im(:)*kxc(ir,2)+r0(:)*coeffim_grho
637 : end do
638 : else
639 233320 : do ir=1,nfft
640 933120 : do ii=1,3 ! grad of spin-up ans spin_dn GS rho
641 933120 : r0_up(ii)=kxc(ir,13+2*ii);r0_dn(ii)=kxc(ir,12+2*ii)-kxc(ir,13+2*ii)
642 : end do
643 933120 : r0(:)=r0_up(:)+r0_dn(:) ! grad of GS rho
644 933120 : r1_up(:)=rho1now(2*ir-1,1,2:4) ! grad of spin-up rho1
645 933120 : r1im_up(:)=rho1now(2*ir,1,2:4) ! grad of spin-up rho1 , im part
646 933120 : r1_dn(:)=rho1now(2*ir-1,2,2:4) ! grad of spin-down rho1
647 933120 : r1im_dn(:)=rho1now(2*ir,2,2:4) ! grad of spin-down rho1 , im part
648 933120 : r1(:)=r1_up(:)+r1_dn(:) ! grad of GS rho1
649 933120 : r1im(:)=r1im_up(:)+r1im_dn(:) ! grad of GS rho1, im part
650 933120 : gradrho_gradrho1_up =dot_product(r0_up,r1_up)
651 933120 : gradrho_gradrho1_dn =dot_product(r0_dn,r1_dn)
652 933120 : gradrho_gradrho1 =dot_product(r0,r1)
653 933120 : gradrho_gradrho1im_up=dot_product(r0_up,r1im_up)
654 933120 : gradrho_gradrho1im_dn=dot_product(r0_dn,r1im_dn)
655 933120 : gradrho_gradrho1im =dot_product(r0,r1im)
656 : dnexcdn(2*ir-1,1)=kxc(ir, 1)*rho1now(2*ir-1,1,1) &
657 : & +kxc(ir, 2)*rho1now(2*ir-1,2,1) &
658 : & +kxc(ir, 6)*gradrho_gradrho1_up &
659 233280 : & +kxc(ir,11)*gradrho_gradrho1
660 : dnexcdn(2*ir-1,2)=kxc(ir, 3)*rho1now(2*ir-1,2,1) &
661 : & +kxc(ir, 2)*rho1now(2*ir-1,1,1) &
662 : & +kxc(ir, 7)*gradrho_gradrho1_dn &
663 233280 : & +kxc(ir,12)*gradrho_gradrho1
664 : dnexcdn(2*ir ,1)=kxc(ir, 1)*rho1now(2*ir ,1,1) &
665 : & +kxc(ir, 2)*rho1now(2*ir ,2,1) &
666 : & +kxc(ir, 6)*gradrho_gradrho1im_up &
667 233280 : & +kxc(ir,11)*gradrho_gradrho1im
668 : dnexcdn(2*ir ,2)=kxc(ir, 3)*rho1now(2*ir ,2,1) &
669 : & +kxc(ir, 2)*rho1now(2*ir ,1,1) &
670 : & +kxc(ir, 7)*gradrho_gradrho1im_dn &
671 233280 : & +kxc(ir,12)*gradrho_gradrho1im
672 : coeff_grho_corr =kxc(ir,11)*rho1now(2*ir-1,1,1) &
673 : & +kxc(ir,12)*rho1now(2*ir-1,2,1) &
674 233280 : & +kxc(ir,13)*gradrho_gradrho1
675 : coeffim_grho_corr=kxc(ir,11)*rho1now(2*ir ,1,1) &
676 : & +kxc(ir,12)*rho1now(2*ir ,2,1) &
677 233280 : & +kxc(ir,13)*gradrho_gradrho1im
678 233280 : coeff_grho_up =kxc(ir,6)*rho1now(2*ir-1,1,1)+kxc(ir,8)*gradrho_gradrho1_up
679 233280 : coeff_grho_dn =kxc(ir,7)*rho1now(2*ir-1,2,1)+kxc(ir,9)*gradrho_gradrho1_dn
680 233280 : coeffim_grho_up=kxc(ir,6)*rho1now(2*ir ,1,1)+kxc(ir,8)*gradrho_gradrho1im_up
681 233280 : coeffim_grho_dn=kxc(ir,7)*rho1now(2*ir ,2,1)+kxc(ir,9)*gradrho_gradrho1im_dn
682 : ! Reuse the storage in rho1now
683 : rho1now(2*ir-1,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1_up(:) &
684 : & +kxc(ir,10) *r1_dn(:) &
685 : & +coeff_grho_up *r0_up(:) &
686 933120 : & +coeff_grho_corr*r0(:)
687 : rho1now(2*ir-1,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1_dn(:) &
688 : & +kxc(ir,10) *r1_up(:) &
689 : & +coeff_grho_dn *r0_dn(:) &
690 933120 : & +coeff_grho_corr*r0(:)
691 : rho1now(2*ir ,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1im_up(:) &
692 : & +kxc(ir,10) *r1im_dn(:) &
693 : & +coeffim_grho_up *r0_up(:) &
694 933120 : & +coeffim_grho_corr *r0(:)
695 : rho1now(2*ir ,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1im_dn(:) &
696 : & +kxc(ir,10) *r1im_up(:) &
697 : & +coeffim_grho_dn *r0_dn(:) &
698 933160 : & +coeffim_grho_corr *r0(:)
699 : end do
700 : end if ! nspden
701 :
702 : end if
703 :
704 193932896 : vxc1(:,:)=zero
705 : call xcpot(cplex,gprimd,ishift,use_laplacian,mpi_enreg,nfft,ngfft,ngrad,nspden,&
706 5381 : & nspgrad,qphon,depsxc=dnexcdn,rhonow=rho1now,vxc=vxc1)
707 :
708 : !call filterpot(paral_kgb,cplex,gmet,gsqcut,nfft,ngfft,nspden,qphon,vxc1)
709 :
710 5381 : ABI_FREE(dnexcdn)
711 5381 : ABI_FREE(rho1now)
712 :
713 : DBG_EXIT("COLL")
714 :
715 5381 : end subroutine dfpt_mkvxcgga
716 : !!***
717 :
718 : !!****f* ABINIT/dfpt_mkvxc_noncoll
719 : !! NAME
720 : !! dfpt_mkvxc_noncoll
721 : !!
722 : !! FUNCTION
723 : !! Compute the first-order change of exchange-correlation potential
724 : !! due to atomic displacement for non-collinear spins: assemble the first-order
725 : !! density change with the frozen-core density change, then use
726 : !! the exchange-correlation kernel.
727 : !!
728 : !! INPUTS
729 : !! cplex= if 1, real space 1-order functions on FFT grid are REAL,
730 : !! if 2, COMPLEX
731 : !! ixc= choice of exchange-correlation scheme
732 : !! ixcrot= option for rotation of collinear spin potential to non collinear full matrix
733 : !! kxc(nfft,nkxc)=exchange and correlation kernel (see rhotoxc.F90)
734 : !! mpi_enreg=information about MPI parallelization
735 : !! nfft=(effective) number of FFT grid points (for this processor)
736 : !! ngfft(18)=contain all needed information about 3D FFT,
737 : !! see ~abinit/doc/variables/vargs.htm#ngfft
738 : !! nhat(nfft,nspden*nhatdim)= -PAW only- GS compensation density
739 : !! nhatdim= -PAW only- 1 if nhat array is used ; 0 otherwise
740 : !! nhat1(cplex*nfft,nspden*nhat1dim)= -PAW only- 1st-order compensation density
741 : !! nhat1dim= -PAW only- 1 if nhat1 array is used ; 0 otherwise
742 : !! nhat1gr(cplex*nfft,nspden,3*nhat1grdim)= -PAW only- gradients of 1st-order compensation density
743 : !! nhat1grdim= -PAW only- 1 if nhat1gr array is used ; 0 otherwise
744 : !! nkxc=second dimension of the kxc array
745 : !! non_magnetic_xc= if true, handle density/potential as non-magnetic (even if it is)
746 : !! nspden=number of spin-density components
747 : !! n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used, otherwise, nfft
748 : !! optnc=option for non-collinear magnetism (nspden=4):
749 : !! 1: the whole 2x2 Vres matrix is computed
750 : !! 2: only Vres^{11} and Vres^{22} are computed
751 : !! option=if 0, work only with the XC core-correction,
752 : !! if 1, treat both density change and XC core correction
753 : !! if 2, treat only density change
754 : !! qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
755 : !! rhor(nfft,nspden)=GS electron density in real space
756 : !! rhor1(cplex*nfft,nspden)=1st-order electron density in real space
757 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
758 : !! usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
759 : !! vxc(nfft,nspden)=GS XC potential
760 : !!
761 : !!
762 : !! OUTPUT
763 : !! vxc1(cplex*nfft,nspden)=change in exchange-correlation potential (including
764 : !! core-correction, if applicable)
765 : !!
766 : !! SOURCE
767 :
768 1086 : subroutine dfpt_mkvxc_noncoll(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat,nhatdim,nhat1,nhat1dim,&
769 1086 : & nhat1gr,nhat1grdim,nkxc,non_magnetic_xc,nspden,n3xccc,optnc,option,qphon,&
770 1086 : & rhor,rhor1,rprimd,usexcnhat,vxc,vxc1,xccc3d1,ixcrot)
771 :
772 : !Arguments ------------------------------------
773 : !scalars
774 : integer,intent(in) :: cplex,ixc,n3xccc,nfft,nhatdim,nhat1dim,nhat1grdim,optnc
775 : integer,intent(in) :: nkxc,nspden,option,usexcnhat
776 : logical,intent(in) :: non_magnetic_xc
777 : type(MPI_type),intent(in) :: mpi_enreg
778 : !arrays
779 : integer,intent(in) :: ngfft(18)
780 : real(dp),intent(in) :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim)
781 : real(dp),intent(in) :: kxc(nfft,nkxc)
782 : real(dp),intent(in) :: vxc(nfft,nspden)
783 : real(dp),intent(in) :: nhat(nfft,nspden*nhatdim),nhat1(cplex*nfft,nspden*nhat1dim)
784 : real(dp),intent(in),target :: rhor(nfft,nspden),rhor1(cplex*nfft,nspden)
785 : real(dp),intent(in) :: qphon(3),rprimd(3,3),xccc3d1(cplex*n3xccc)
786 : real(dp),intent(out) :: vxc1(cplex*nfft,nspden)
787 : integer,optional,intent(in) :: ixcrot
788 : !Local variables-------------------------------
789 : !scalars
790 : !arrays
791 : real(dp) :: nhat1_zero(0,0),nhat1gr_zero(0,0,0),tsec(2)
792 1086 : real(dp),allocatable :: m_norm(:),rhor1_diag(:,:),vxc1_diag(:,:)
793 1086 : real(dp), contiguous, pointer :: mag(:,:),rhor_(:,:),rhor1_(:,:)
794 : ! *************************************************************************
795 :
796 : ! Non-collinear magnetism
797 : ! Has to locally "rotate" rho(r)^(1) (according to magnetization),
798 : ! Compute Vxc(r)^(1) in the spin frame aligned with \vec{m} and rotate it back
799 :
800 : DBG_ENTER("COLL")
801 : ABI_UNUSED(nhat1gr)
802 :
803 1086 : call timab(181,1,tsec)
804 :
805 1086 : if(nspden/=4) then
806 0 : ABI_BUG('only for nspden=4!')
807 : end if
808 :
809 1086 : if(nkxc/=2*min(nspden,2)-1) then
810 0 : ABI_BUG('nspden=4 works only with LSDA.')
811 : end if
812 :
813 : !Special case: no XC applied
814 1086 : if (ixc==0.or.nkxc==0) then
815 0 : ABI_WARNING('Note that no xc is applied (ixc=0)')
816 0 : vxc1(:,:)=zero
817 : return
818 : end if
819 :
820 :
821 :
822 : !Treat first LDA
823 1086 : if(nkxc==1.or.nkxc==3)then
824 :
825 32038638 : vxc1(:,:)=zero
826 :
827 : ! PAW: possibly substract compensation density
828 1086 : if ((usexcnhat==0.and.nhatdim==1).or.(non_magnetic_xc)) then
829 40 : ABI_MALLOC(rhor_,(nfft,nspden))
830 10 : if (usexcnhat==0.and.nhatdim==1) then
831 0 : rhor_(:,:) =rhor(:,:)-nhat(:,:)
832 : else
833 1106010 : rhor_(:,:) =rhor(:,:)
834 : end if
835 10 : if (non_magnetic_xc) then
836 10 : if(nspden==2) rhor_(:,2)=rhor_(:,1)*half
837 414760 : if(nspden==4) rhor_(:,2:4)=zero
838 : end if
839 : else
840 1076 : rhor_ => rhor
841 : end if
842 1086 : if ((usexcnhat==0.and.nhat1dim==1).or.(non_magnetic_xc)) then
843 40 : ABI_MALLOC(rhor1_,(cplex*nfft,nspden))
844 10 : if (usexcnhat==0.and.nhatdim==1) then
845 0 : rhor1_(:,:)=rhor1(:,:)-nhat1(:,:)
846 : else
847 1106010 : rhor1_(:,:)=rhor1(:,:)
848 : end if
849 10 : if (non_magnetic_xc) then
850 10 : if(nspden==2) rhor1_(:,2)=rhor1_(:,1)*half
851 414760 : if(nspden==4) rhor1_(:,2:4)=zero
852 : end if
853 : else
854 1076 : rhor1_ => rhor1
855 : end if
856 :
857 : ! Magnetization
858 1086 : mag => rhor_(:,2:4)
859 3258 : ABI_MALLOC(rhor1_diag,(cplex*nfft,2))
860 2172 : ABI_MALLOC(vxc1_diag,(cplex*nfft,2))
861 3258 : ABI_MALLOC(m_norm,(nfft))
862 :
863 : ! -- Rotate rho(r)^(1)
864 : ! SPr: for option=0 the rhor is not used, only core density xccc3d1
865 : ! rotate_mag is only to compute the m_norm
866 : call rotate_mag(rhor1_,rhor1_diag,mag,nfft,cplex,mag_norm_out=m_norm,&
867 1086 : & rho_out_format=2)
868 :
869 : ! -- Compute Vxc(r)^(1)=Kxc(r).rho(r)^(1)_rotated
870 : ! Note for PAW: nhat has already been substracted; don't use it in dfpt_mkvxc
871 : ! (put all nhat options to zero).
872 : ! The collinear routine dfpt_mkvxc wants a general density built as (tr[rho],rho_upup)
873 : call dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1_zero,0,nhat1gr_zero,0,&
874 1086 : & nkxc,non_magnetic_xc,2,n3xccc,option,qphon,rhor1_diag,rprimd,0,vxc1_diag,xccc3d1)
875 :
876 : !call test_rotations(0,1)
877 :
878 : ! -- Rotate back Vxc(r)^(1)
879 1086 : if (optnc==1) then
880 997 : if(present(ixcrot)) then
881 : call rotate_back_mag_dfpt(option,vxc1_diag,vxc1,vxc,kxc,rhor1_,mag,nfft,cplex,&
882 997 : & mag_norm_in=m_norm,rot_method=ixcrot)
883 : else
884 : call rotate_back_mag_dfpt(option,vxc1_diag,vxc1,vxc,kxc,rhor1_,mag,nfft,cplex,&
885 0 : & mag_norm_in=m_norm)
886 : end if
887 : else
888 89 : call rotate_back_mag(vxc1_diag,vxc1,mag,nfft,mag_norm_in=m_norm)
889 1817961 : vxc1(:,3:4)=zero
890 : end if
891 :
892 1086 : ABI_FREE(rhor1_diag)
893 1086 : ABI_FREE(vxc1_diag)
894 1086 : ABI_FREE(m_norm)
895 1086 : if ((usexcnhat==0.and.nhatdim==1).or.(non_magnetic_xc)) then
896 10 : ABI_FREE(rhor_)
897 : end if
898 1086 : if ((usexcnhat==0.and.nhat1dim==1).or.(non_magnetic_xc)) then
899 10 : ABI_FREE(rhor1_)
900 : end if
901 :
902 : end if ! nkxc=1 or nkxc=3
903 :
904 1086 : call timab(181,2,tsec)
905 :
906 : DBG_EXIT("COLL")
907 :
908 1086 : end subroutine dfpt_mkvxc_noncoll
909 : !!***
910 :
911 : !!****f* ABINIT/dfpt_mkvxcggadq
912 : !! NAME
913 : !! dfpt_mkvxcggadq
914 : !!
915 : !! FUNCTION
916 : !! Compute the first-order change of exchange-correlation potential
917 : !! in case of GGA functionals
918 : !! Use the q-gradient (Cartesian) of the exchange-correlation kernel.
919 : !!
920 : !! INPUTS
921 : !! cplex= if 1, real space 1-order functions on FFT grid are REAL,
922 : !! if 2, COMPLEX
923 : !! gmet(3,3)=metrix tensor in G space in Bohr**-2.
924 : !! gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1)
925 : !! gsqcut=cutoff value on G**2 for sphere inside fft box.
926 : !! kxc(nfft,nkxc)=exchange and correlation kernel (see below)
927 : !! mpi_enreg=information about MPI parallelization
928 : !! nfft=(effective) number of FFT grid points (for this processor)
929 : !! ngfft(18)=contain all needed information about 3D FFT
930 : !! nkxc=second dimension of the kxc array
931 : !! nspden=number of spin-density components
932 : !! qdirc= indicates the Cartesian direction of the q-gradient (1,2 or 3)
933 : !! rhor1tmp(cplex*nfft,2)=array for first-order electron spin-density
934 : !! in electrons/bohr**3 (second index corresponds to spin-up and spin-down)
935 : !!
936 : !! OUTPUT
937 : !! vxc1(2*nfft,nspden)=change in exchange-correlation potential
938 : !!
939 : !! NOTES
940 : !! For the time being, a rather crude coding, to be optimized ...
941 : !! Content of Kxc array:
942 : !! Only works with nspden=1
943 : !! ===== if GGA
944 : !! if nspden==1:
945 : !! kxc(:,1)= d2Exc/drho2
946 : !! kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)|
947 : !! kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho
948 : !! kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| )
949 : !! kxc(:,5)= gradx(rho)
950 : !! kxc(:,6)= grady(rho)
951 : !! kxc(:,7)= gradz(rho)
952 : !!
953 : !! SOURCE
954 :
955 669 : subroutine dfpt_mkvxcggadq(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,&
956 669 : & nkxc,nspden,qdirc,rhor1,vxc1)
957 :
958 : !Arguments ------------------------------------
959 : !scalars
960 : integer,intent(in) :: cplex,nfft,nkxc,nspden,qdirc
961 : type(MPI_type),intent(in) :: mpi_enreg
962 : !arrays
963 : integer,intent(in) :: ngfft(18)
964 : real(dp),intent(in) :: gprimd(3,3)
965 : real(dp),intent(in) :: kxc(nfft,nkxc)
966 : real(dp),intent(in),target :: rhor1(cplex*nfft,nspden)
967 : real(dp),intent(out) :: vxc1(2*nfft,nspden)
968 :
969 : !Local variables-------------------------------
970 : !scalars
971 : integer :: ii,ir,ishift,ngrad,nspgrad
972 : real(dp) :: gradrho_gradrho1
973 : !character(len=500) :: msg
974 : !arrays
975 : real(dp) :: qphon(3)
976 : real(dp) :: r0(3),r1(3)
977 669 : real(dp),allocatable :: ar1(:,:)
978 669 : real(dp),allocatable :: a_gradi_r1(:,:)
979 669 : real(dp),allocatable :: dadgradn_t1(:,:,:),dadgradn_t2(:,:)
980 669 : real(dp),allocatable :: rho1now(:,:,:)
981 669 : real(dp),contiguous, pointer :: rhor1_ptr(:,:)
982 : ! *************************************************************************
983 :
984 : DBG_EXIT("COLL")
985 :
986 669 : if (nkxc/=7) then
987 0 : ABI_BUG('Wrong nkxc value for GGA in the longwave driver (optdriver=10)!')
988 : end if
989 :
990 : !Compute the gradients of the first-order density
991 : !rho1now(:,:,1) contains the first-order density, and
992 : !rho1now(:,:,2:4) contains the gradients of the first-order density
993 669 : ishift=0 ; ngrad=2
994 669 : qphon(:)=zero
995 669 : rhor1_ptr => rhor1
996 3345 : ABI_MALLOC(rho1now,(cplex*nfft,nspden,ngrad*ngrad))
997 669 : call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden,qphon,rhor1_ptr,rho1now)
998 :
999 : !Apply the XC kernel
1000 669 : nspgrad=1
1001 2007 : ABI_MALLOC(ar1,(cplex*nfft,nspgrad))
1002 1338 : ABI_MALLOC(a_gradi_r1,(cplex*nfft,nspgrad))
1003 2676 : ABI_MALLOC(dadgradn_t1,(cplex*nfft,nspgrad,3))
1004 1338 : ABI_MALLOC(dadgradn_t2,(cplex*nfft,nspgrad))
1005 2362368 : do ir=1,nfft
1006 16531893 : r0(:)=kxc(ir,5:7); r1(:)=rho1now(ir,1,2:4)
1007 9446796 : gradrho_gradrho1=dot_product(r0,r1)
1008 2361699 : ar1(ir,1)=kxc(ir,2)*rho1now(ir,1,1)
1009 2361699 : a_gradi_r1(ir,1)=kxc(ir,2)*r1(qdirc)
1010 2361699 : dadgradn_t2(ir,1)=kxc(ir,4)*gradrho_gradrho1*r0(qdirc)
1011 9447465 : dadgradn_t1(ir,1,:)=kxc(ir,4)*r0(:)*r0(qdirc)*rho1now(ir,1,1)
1012 : end do
1013 2676 : do ii=1,3
1014 2364375 : if (ii==qdirc) dadgradn_t1(:,1,ii)=dadgradn_t1(:,1,ii)+ar1(:,1)
1015 : end do
1016 :
1017 : !Incorporate the terms that do not need further treatment
1018 : !(a -i factor is applied here)
1019 2362368 : do ir=1,nfft
1020 2361699 : ii=2*ir
1021 2361699 : vxc1(ii-1,1)=zero
1022 2362368 : vxc1(ii,1)= -a_gradi_r1(ir,1) -dadgradn_t2(ir,1)
1023 : end do
1024 669 : ABI_FREE(rho1now)
1025 669 : ABI_FREE(a_gradi_r1)
1026 669 : ABI_FREE(dadgradn_t2)
1027 669 : ABI_FREE(ar1)
1028 :
1029 : !Now the term whose sum over real-space derivatives has to be computed
1030 : call xcpotdq(dadgradn_t1,cplex,gprimd,ishift,mpi_enreg,nfft, &
1031 669 : & ngfft,ngrad,nspden,nspgrad,vxc1)
1032 :
1033 669 : ABI_FREE(dadgradn_t1)
1034 :
1035 669 : end subroutine dfpt_mkvxcggadq
1036 : !!***
1037 :
1038 : !!****f* ABINIT/dfpt_mkvxcgga_n0met
1039 : !! NAME
1040 : !! dfpt_mkvxcgga_n0met
1041 : !!
1042 : !! FUNCTION
1043 : !! Compute the contribution to the second q-gradient of the metric
1044 : !! perturbation that comes from gga XC potentials and depends only
1045 : !! on ground state rho
1046 : !!
1047 : !! INPUTS
1048 : !! beta= indicates the Cartesian direction of the metric perturbation
1049 : !! cplex= if 1, real space 1-order functions on FFT grid are REAL,
1050 : !! if 2, COMPLEX
1051 : !! delta= indicates the Cartesian direction of the first q-gradient
1052 : !! gamma= indicates the Cartesian direction of the second q-gradient
1053 : !! gmet(3,3)=metrix tensor in G space in Bohr**-2.
1054 : !! gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1)
1055 : !! gsqcut=cutoff value on G**2 for sphere inside fft box.
1056 : !! kxc(nfft,nkxc)=exchange and correlation kernel (see below)
1057 : !! mpi_enreg=information about MPI parallelization
1058 : !! nfft=(effective) number of FFT grid points (for this processor)
1059 : !! ngfft(18)=contain all needed information about 3D FFT
1060 : !! nkxc=second dimension of the kxc array
1061 : !! nspden=number of spin-density components
1062 : !! rho(cplex*nfft,2)=array for ground-state electron spin-density
1063 : !! in electrons/bohr**3 (second index corresponds to spin-up and spin-down)
1064 : !!
1065 : !! OUTPUT
1066 : !! vxc1(2*nfft,nspden)=change in exchange-correlation potential
1067 : !!
1068 : !! NOTES
1069 : !! For the time being, a rather crude coding, to be optimized ...
1070 : !! Content of Kxc array:
1071 : !! Only works with nspden=1
1072 : !! ===== if GGA
1073 : !! if nspden==1:
1074 : !! kxc(:,1)= d2Exc/drho2
1075 : !! kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)|
1076 : !! kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho
1077 : !! kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| )
1078 : !! kxc(:,5)= gradx(rho)
1079 : !! kxc(:,6)= grady(rho)
1080 : !! kxc(:,7)= gradz(rho)
1081 : !!
1082 : !! SOURCE
1083 :
1084 243 : subroutine dfpt_mkvxcgga_n0met(beta,cplex,delta,gamma,gprimd,kxc,mpi_enreg,nfft,ngfft,&
1085 243 : & nkxc,nspden,rhor,vxc1)
1086 :
1087 : !Arguments ------------------------------------
1088 : !scalars
1089 : integer,intent(in) :: beta,cplex,delta,gamma,nfft,nkxc,nspden
1090 : type(MPI_type),intent(in) :: mpi_enreg
1091 : !arrays
1092 : integer,intent(in) :: ngfft(18)
1093 : real(dp),intent(in) :: gprimd(3,3)
1094 : real(dp),intent(in) :: kxc(nfft,nkxc)
1095 : real(dp),intent(in),target :: rhor(cplex*nfft,nspden)
1096 : real(dp),intent(out) :: vxc1(2*nfft,nspden)
1097 :
1098 : !Local variables-------------------------------
1099 : !scalars
1100 : integer :: alpha,ii,ir,ishift,ngrad,nspgrad
1101 : real(dp) :: delag,delad,delbd,delbg,deldg
1102 : real(dp) :: gmodsq
1103 : character(len=500) :: msg
1104 : !arrays
1105 : real(dp) :: r0(3)
1106 243 : real(dp),allocatable :: dadgg(:,:),dadgtgn(:,:),gna(:,:),dadgngn_1(:,:),dadgngn_2(:,:)
1107 243 : real(dp),allocatable :: dadgngn(:,:,:),kro_an(:,:,:),sumgrad(:,:,:)
1108 : ! *************************************************************************
1109 :
1110 : DBG_EXIT("COLL")
1111 :
1112 243 : if (nkxc/=7) then
1113 0 : msg='Wrong nkxc value for GGA in the longwave driver (optdriver=10)!'
1114 0 : ABI_BUG(msg)
1115 : end if
1116 :
1117 : !Kronecker deltas
1118 243 : delbd=0.0_dp; delbg=0.0_dp; deldg=0.0_dp
1119 243 : if (beta==delta) delbd=1.0_dp
1120 243 : if (beta==gamma) delbg=1.0_dp
1121 243 : if (delta==gamma) deldg=1.0_dp
1122 :
1123 : !Apply the XC kernel
1124 243 : nspgrad=1
1125 729 : ABI_MALLOC(dadgg,(cplex*nfft,nspgrad))
1126 486 : ABI_MALLOC(dadgtgn,(cplex*nfft,nspgrad))
1127 486 : ABI_MALLOC(gna,(cplex*nfft,nspgrad))
1128 486 : ABI_MALLOC(dadgngn_1,(cplex*nfft,nspgrad))
1129 486 : ABI_MALLOC(dadgngn_2,(cplex*nfft,nspgrad))
1130 820368 : do ir=1,nfft
1131 3280500 : r0(:)=kxc(ir,5:7)
1132 820125 : gmodsq=r0(1)**2+r0(2)**2+r0(3)**2
1133 820125 : dadgg(ir,1)=kxc(ir,4)*gmodsq*(delbd*r0(gamma)+delbg*r0(delta))
1134 820125 : dadgtgn(ir,1)=two*kxc(ir,4)*r0(beta)*r0(delta)*r0(gamma)
1135 820125 : gna(ir,1)=(delbg*r0(delta)+delbd*r0(gamma)+two*deldg*r0(beta))*kxc(ir,2)
1136 820125 : dadgngn_1(ir,1)=delbd*kxc(ir,4)*rhor(ir,1)*r0(gamma)
1137 820368 : dadgngn_2(ir,1)=delbg*kxc(ir,4)*rhor(ir,1)*r0(delta)
1138 : end do
1139 :
1140 : !Incorporate the terms that do not need further treatment
1141 820368 : do ir=1,nfft
1142 820125 : ii=2*ir
1143 820125 : vxc1(ii-1,1)= -dadgg(ir,1)-dadgtgn(ir,1)-gna(ir,1)
1144 820368 : vxc1(ii,1)= zero
1145 : end do
1146 243 : ABI_FREE(dadgg)
1147 243 : ABI_FREE(dadgtgn)
1148 243 : ABI_FREE(gna)
1149 :
1150 : !Build the last term whose gradient needs to be computed
1151 972 : ABI_MALLOC(dadgngn,(cplex*nfft,nspgrad,3))
1152 486 : ABI_MALLOC(kro_an,(cplex*nfft,nspgrad,3))
1153 486 : ABI_MALLOC(sumgrad,(cplex*nfft,nspgrad,3))
1154 972 : do alpha=1,3
1155 729 : delad=0.0_dp; delag=0.0_dp
1156 729 : if (alpha==delta) delad=1.0_dp
1157 729 : if (alpha==gamma) delag=1.0_dp
1158 2461347 : do ir=1,nfft
1159 9841500 : r0(:)=kxc(ir,5:7)
1160 2460375 : dadgngn(ir,1,alpha)=(dadgngn_1(ir,1)+dadgngn_2(ir,1))*r0(alpha)
1161 2460375 : kro_an(ir,1,alpha)=(delbd*delag+delbg*delad)*rhor(ir,1)*kxc(ir,2)
1162 2461104 : sumgrad(ir,1,alpha)=dadgngn(ir,1,alpha)+kro_an(ir,1,alpha)
1163 : end do
1164 : end do
1165 :
1166 243 : ABI_FREE(dadgngn_1)
1167 243 : ABI_FREE(dadgngn_2)
1168 243 : ABI_FREE(dadgngn)
1169 243 : ABI_FREE(kro_an)
1170 :
1171 : !Now the term whose sum over real-space derivatives has to be computed.
1172 : !(Use the same routine as in the q-gradient of the XC kernel. It saves
1173 : ! the gradient sum in the imaginary part of vxc1 and includes an additional
1174 : ! two_pi factor. Need to fix this after the call.)
1175 243 : ishift=0 ; ngrad=2
1176 : call xcpotdq(sumgrad,cplex,gprimd,ishift,mpi_enreg,nfft, &
1177 243 : & ngfft,ngrad,nspden,nspgrad,vxc1)
1178 :
1179 820368 : do ir=1,nfft
1180 820125 : ii=2*ir
1181 820125 : vxc1(ii-1,1)=vxc1(ii-1,1)+vxc1(ii,1)/two_pi
1182 820368 : vxc1(ii,1)=zero
1183 : end do
1184 :
1185 :
1186 243 : ABI_FREE(sumgrad)
1187 :
1188 243 : end subroutine dfpt_mkvxcgga_n0met
1189 : !!***
1190 :
1191 : !!****f* ABINIT/dfpt_mkvxcccdq
1192 : !! NAME
1193 : !! dfpt_mkvxcccdq
1194 : !!
1195 : !! FUNCTION
1196 : !! Computes the first q-gradient of the first-order exchange-correlation potential
1197 : !! due to the pseudocore density.
1198 : !!
1199 : !! INPUTS
1200 : !! cplex= if 1, real space 1-order functions on FFT grid are REAL,
1201 : !! if 2, COMPLEX
1202 : !! i3dir= reduced direction of the q-gradient
1203 : !! ixc= choice of exchange-correlation scheme
1204 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2
1205 : !! gprimd(3,3)=reciprocal space dimensional primitive translations
1206 : !! kxc(nfft,nkxc)=exchange and correlation kernel
1207 : !! mpi_enreg=information about MPI parallelization
1208 : !! nfft=(effective) number of FFT grid points (for this processor)
1209 : !! ngfft(1:18)=integer array with FFT box dimensions and other
1210 : !! nspden=number of spin-density components
1211 : !! nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed.
1212 : !! qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
1213 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
1214 : !! xccc3d1(cplex*nfft)=3D change in core charge density
1215 : !! xccc3d1dq(2*nfft)=q_i3dir-gradient of 3D change in core charge density
1216 : !!
1217 : !! OUTPUT
1218 : !! vxccc1dq(2*nfft,nspden)= q-gradient of first-order XC potential due to pseudocore charge
1219 : !!
1220 : !! SOURCE
1221 :
1222 27 : subroutine dfpt_mkvxcccdq(cplex,i3dir,ixc,gprimd,kxc,mpi_enreg,nfft, &
1223 27 : & ngfft,nkxc,nspden,qphon,rprimd,vxccc1dq,xccc3d1,xccc3d2dq)
1224 :
1225 : !Arguments ------------------------------------
1226 : !scalars
1227 : integer , intent(in) :: cplex,i3dir,ixc,nfft,nkxc,nspden
1228 : type(MPI_type),intent(inout) :: mpi_enreg
1229 :
1230 : !arrays
1231 : integer,intent(in) :: ngfft(18)
1232 : real(dp), intent(in) :: gprimd(3,3)
1233 : real(dp), intent(in) :: kxc(nfft,nkxc)
1234 : real(dp), intent(in) :: qphon(3),rprimd(3,3)
1235 : real(dp), intent(in) :: xccc3d1(cplex*nfft)
1236 : real(dp), intent(in) :: xccc3d2dq(2*nfft)
1237 : real(dp), intent(out) :: vxccc1dq(2*nfft,nspden)
1238 :
1239 : !Local variables-------------------------------
1240 : !scalars
1241 : integer :: ii,ispden,ir,jj,nhat1grdim,option,qcar,usexcnhat,usepaw
1242 : real(dp) :: spin_scale
1243 : logical :: non_magnetic_xc
1244 : !arrays
1245 27 : real(dp),allocatable :: nhat1(:,:),nhat1gr(:,:,:)
1246 27 : real(dp),allocatable :: rhor1(:,:), rhor1_cplx(:,:)
1247 27 : real(dp),allocatable :: vxc1dq_a(:,:),vxc1dq_b(:,:),vxc1dq_car(:,:,:)
1248 : ! *************************************************************************
1249 :
1250 : DBG_ENTER("COLL")
1251 :
1252 244662 : vxccc1dq= zero
1253 :
1254 : !If GGA xc first calculate the contribution from the q gradient of the xc potential
1255 27 : if (nkxc == 7) then
1256 :
1257 : !Adapt the format of xccc3d1
1258 96 : ABI_MALLOC(rhor1,(cplex*nfft,nspden))
1259 24 : spin_scale=one;if (nspden==2) spin_scale=half
1260 48 : do ispden=1,nspden
1261 98352 : do ir=1,cplex*nfft
1262 98328 : rhor1(ir,ispden)=xccc3d1(ir)*spin_scale
1263 : end do
1264 : end do
1265 :
1266 : !The gradient of the potential is calculated in Cartesian coordinates
1267 120 : ABI_MALLOC(vxc1dq_a,(2*nfft,nspden))
1268 120 : ABI_MALLOC(vxc1dq_car,(2*nfft,nspden,3))
1269 96 : do qcar=1,3
1270 72 : call dfpt_mkvxcggadq(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,qcar,rhor1,vxc1dq_a)
1271 :
1272 : !Here we apply an i factor, to compensate the lake of the -i factor in
1273 : !vxc1dq_b (see notes in dfpt_vlocaldq).
1274 295008 : do ir=1,nfft
1275 294912 : ii=2*ir-1
1276 294912 : jj=2*ir
1277 589824 : vxc1dq_car(ii,:,qcar)=-vxc1dq_a(jj,:)
1278 589896 : vxc1dq_car(jj,:,qcar)= vxc1dq_a(ii,:)
1279 : end do
1280 : end do
1281 24 : ABI_FREE(rhor1)
1282 :
1283 : !Convert to reduced coordinate i3dir
1284 196656 : vxc1dq_a=zero
1285 96 : do qcar=1,3
1286 589992 : vxc1dq_a(:,:)=vxc1dq_a(:,:) + gprimd(qcar,i3dir) * vxc1dq_car(:,:,qcar)
1287 : end do
1288 24 : ABI_FREE(vxc1dq_car)
1289 :
1290 : !Accumulate this term
1291 196656 : vxccc1dq= vxc1dq_a
1292 :
1293 24 : ABI_FREE(vxc1dq_a)
1294 : end if
1295 :
1296 : !Calculate the term with the gradient of the first-order pseudocore density
1297 : !Dummy arguments for mkvxc
1298 135 : ABI_MALLOC(rhor1_cplx,(2*nfft,nspden))
1299 27 : usexcnhat= 0
1300 27 : nhat1grdim= 0
1301 27 : ABI_MALLOC(nhat1gr,(0,0,0))
1302 : nhat1gr(:,:,:)= zero
1303 27 : usepaw= 0
1304 54 : ABI_MALLOC(nhat1,(2*nfft,nspden*usepaw))
1305 : nhat1= zero
1306 27 : non_magnetic_xc= .true.
1307 27 : option= 0
1308 81 : ABI_MALLOC(vxc1dq_b,(2*nfft,nspden))
1309 : call dfpt_mkvxc(2,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,&
1310 27 : & non_magnetic_xc,nspden,nfft,option,qphon,rhor1_cplx,rprimd,usexcnhat,vxc1dq_b,xccc3d2dq)
1311 :
1312 : !Accumulate this term
1313 244662 : vxccc1dq= vxccc1dq + vxc1dq_b
1314 :
1315 : !Deallocations
1316 27 : ABI_FREE(vxc1dq_b)
1317 27 : ABI_FREE(rhor1_cplx)
1318 27 : ABI_FREE(nhat1)
1319 27 : ABI_FREE(nhat1gr)
1320 :
1321 : DBG_EXIT("COLL")
1322 :
1323 27 : end subroutine dfpt_mkvxcccdq
1324 : !!***
1325 :
1326 :
1327 : end module m_dfpt_mkvxc
1328 : !!***
|
