Branch data 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 : 56279 : subroutine dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1,nhat1dim,nhat1gr,nhat1grdim,&
130 : 56279 : & 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 : 56279 : real(dp), contiguous, pointer :: nhat1_(:,:),nhat1gr_(:,:,:),rhor1_(:,:)
155 : : ! *************************************************************************
156 : :
157 : : DBG_ENTER("COLL")
158 : :
159 : 56279 : call timab(181,1,tsec)
160 : :
161 [ - + ]: 56279 : 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 [ + + ]: 56279 : 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 [ + + ]: 49780 : if (option/=0) then
177 [ + + + + : 25478 : 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 : 24264 : rhor1_ => rhor1
190 : : end if
191 : : end if
192 : :
193 : : ! Case without non-linear core correction
194 [ + + + + ]: 49780 : if(n3xccc==0 .or. option==2)then
195 : :
196 [ + + ]: 28794 : 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 [ + + ]: 23753 : if(nspden==1)then
204 [ + + ]: 21539 : if(cplex==1)then
205 [ + + ]: 417586412 : do ir=1,nfft
206 : 417586412 : 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 [ + + ]: 2214 : if(cplex==1)then
218 [ + + ]: 37582331 : do ir=1,nfft
219 : 37580282 : rho1_dn=rhor1_(ir,1)-rhor1_(ir,2)
220 : 37580282 : vxc1(ir,1)=kxc(ir,1)*rhor1_(ir,2)+kxc(ir,2)*rho1_dn
221 : 37582331 : 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 [ + + ]: 20986 : if(option==0)then
241 : :
242 [ + + ]: 19261 : 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 [ + + ]: 2221 : if(cplex==1)then
255 [ + + ]: 19912724 : do ir=1,nfft
256 : 19911064 : vxc1(ir,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(ir)*half
257 : 19912724 : 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 [ + + ]: 1725 : 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 [ + + ]: 503 : if(cplex==1)then
284 [ + + ]: 3079255 : do ir=1,nfft
285 : 3078771 : rho1_dn=rhor1_(ir,1)-rhor1_(ir,2) + xccc3d1(ir)*half
286 : 3078771 : rho1_up=rhor1_(ir,2) + xccc3d1(ir)*half
287 : 3078771 : vxc1(ir,1)=kxc(ir,1)*rho1_up+kxc(ir,2)*rho1_dn
288 : 3079255 : 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 [ + + + + : 49780 : 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 [ - + + - : 22148 : ABI_MALLOC(rhor1_,(cplex*nfft,nspden))
- + + - +
- + - - +
- + ]
319 [ + + ]: 5537 : if (option==1 .or. option==2) then
320 [ + + ]: 2662 : if (nspden==1) then
321 [ + + ]: 158101320 : do ir=1,cplex*nfft
322 : 158101320 : 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 [ + + ]: 5834 : do ispden=1,nspden
341 [ + + ]: 189382562 : do ir=1,cplex*nfft
342 : 189379687 : rhor1_(ir,ispden)=zero
343 : : end do
344 : : end do
345 : : end if
346 : :
347 [ + + + + ]: 5537 : 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 : 5537 : nhat1dim_=nhat1dim ; nhat1rgdim_=nhat1grdim
358 [ + + + + : 5537 : 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 [ + + ]: 5482 : else if (option==0) then
373 [ - + ]: 2875 : ABI_MALLOC(nhat1_,(0,0))
374 : 2875 : nhat1dim_=0
375 : : else
376 : 2607 : nhat1_ => nhat1
377 : : end if
378 [ + + + + : 5537 : 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 [ + + ]: 5529 : else if (option==0) then
397 [ - + ]: 2875 : ABI_MALLOC(nhat1gr_,(0,0,0))
398 : 2875 : nhat1rgdim_=0
399 : : else
400 : 2654 : nhat1gr_ => nhat1gr
401 : : end if
402 : :
403 : 5537 : call matr3inv(rprimd,gprimd)
404 : :
405 : : call dfpt_mkvxcgga(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat1_,nhat1dim_,&
406 : 5537 : & nhat1gr_,nhat1rgdim_,nkxc,nspden,qphon,rhor1_,usexcnhat,vxc1)
407 : :
408 [ - + ]: 5537 : ABI_FREE(rhor1_)
409 [ + + + + : 5537 : if ((option==0).or.(nhat1dim==1.and.nspden==2)) then
+ + ]
410 [ - + ]: 2930 : ABI_FREE(nhat1_)
411 : : end if
412 [ + + + + : 5537 : if ((option==0).or.(nhat1grdim==1.and.nspden==2)) then
+ - ]
413 [ - + ]: 2883 : ABI_FREE(nhat1gr_)
414 : : end if
415 : :
416 : : else
417 : 0 : ABI_BUG('Invalid nkxc!')
418 : :
419 : : end if ! LDA or GGA
420 : :
421 : 55317 : call timab(181,2,tsec)
422 : :
423 : : DBG_EXIT("COLL")
424 : :
425 : 56279 : 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 : 5537 : subroutine dfpt_mkvxcgga(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,&
497 : 5537 : & nhat1,nhat1dim,nhat1gr,nhat1grdim,nkxc,&
498 : 5537 : & 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 : 5537 : real(dp),allocatable :: dnexcdn(:,:),rho1now(:,:,:)
527 : 5537 : real(dp),contiguous, pointer :: rhor1_ptr(:,:)
528 : : ! *************************************************************************
529 : :
530 : : DBG_ENTER("COLL")
531 : :
532 : 5537 : 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 : 5537 : use_laplacian=0
539 : :
540 : : !PAW: substract 1st-order compensation density from 1st-order density
541 [ + + + + : 5537 : test_nhat=((nhat1dim==1).and.(usexcnhat==0.or.nhat1grdim==1))
+ + ]
542 : : if (test_nhat) then
543 [ - + + - : 2728 : ABI_MALLOC(rhor1_ptr,(cplex*nfft,nspden))
- + + - +
- + - - +
- + ]
544 [ + - + + : 240081162 : rhor1_ptr(:,:)=rhor1(:,:)-nhat1(:,:)
+ + + + +
+ ]
545 : : else
546 : 4855 : 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 : 5537 : ishift=0 ; ngrad=2
555 [ - + + - : 27685 : ABI_MALLOC(rho1now,(cplex*nfft,nspden,ngrad*ngrad))
- + + - +
- + - + -
- + - + ]
556 : 5537 : 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 [ + + + + ]: 5537 : if (test_nhat.and.usexcnhat==1) then
560 [ + + + + ]: 186648 : rho1now(:,1:nspden,1)=rho1now(:,1:nspden,1)+nhat1(:,1:nspden)
561 : : end if
562 [ + + ]: 5537 : 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 [ + + ]: 5537 : if (test_nhat) then
568 [ - + ]: 682 : ABI_FREE(rhor1_ptr)
569 : : end if
570 : :
571 : : !Apply the XC kernel
572 [ + + ]: 5537 : nspgrad=2; if (nspden==2) nspgrad=5
573 [ - + + - : 22148 : ABI_MALLOC(dnexcdn,(cplex*nfft,nspgrad))
- + + - +
- + - - +
- + ]
574 : :
575 [ + + ]: 5537 : if (cplex==1) then ! Treat real case first
576 [ + + ]: 4140 : if (nspden==1) then
577 [ + + ]: 317986490 : do ir=1,nfft
578 [ + + + + ]: 2225877731 : r0(:)=kxc(ir,5:7) ; r1(:)=rho1now(ir,1,2:4)
579 [ + + ]: 1271930132 : gradrho_gradrho1=dot_product(r0,r1)
580 : 317982533 : dnexcdn(ir,1)=kxc(ir,1)*rho1now(ir,1,1) + kxc(ir,3)*gradrho_gradrho1
581 : 317982533 : coeff_grho=kxc(ir,3)*rho1now(ir,1,1) + kxc(ir,4)*gradrho_gradrho1
582 : : ! Reuse the storage in rho1now
583 [ + + ]: 1271934089 : 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 [ + + + + ]: 349933208 : vxc1(:,:)=zero
705 : : call xcpot(cplex,gprimd,ishift,use_laplacian,mpi_enreg,nfft,ngfft,ngrad,nspden,&
706 : 5537 : & nspgrad,qphon,depsxc=dnexcdn,rhonow=rho1now,vxc=vxc1)
707 : :
708 : : !call filterpot(paral_kgb,cplex,gmet,gsqcut,nfft,ngfft,nspden,qphon,vxc1)
709 : :
710 [ - + ]: 5537 : ABI_FREE(dnexcdn)
711 [ - + ]: 5537 : ABI_FREE(rho1now)
712 : :
713 : : DBG_EXIT("COLL")
714 : :
715 [ - + ]: 5537 : 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 : : !!***
|
