Line data Source code
1 : !!****m* ABINIT/m_kxc
2 : !! NAME
3 : !! m_kxc
4 : !!
5 : !! FUNCTION
6 : !! Helper functions to compute the XC kernel in reciprocal space.
7 : !! WARNING: At present (10/01/14) these routines are not tested
8 : !! since the ACFD code has been disabled.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 1998-2026 ABINIT group (DCA, MF, XG, GMR, LSI, YMN, Rhaltaf, MS)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !!
16 : !! NOTES:
17 : !! libxc_functionals.F90 uses a global structure (funcs) to store the XC parameters.
18 : !! This structure is initialized in driver with the value of ixc specified by the user in the input file.
19 : !! In order to change the value of ixc at run-time, we have to reinitialize the global structure
20 : !! with the new value of ixc before computing XC quantities.
21 : !! Moreover one has to reinstate the old functional before returning so that the other routines
22 : !! will continue to used the previous ixc. This task can be accomplished with the following pseudocode
23 : !!
24 : !! ! Reinitialize the libxc module with the overridden values
25 : !! if (old_ixc<0) call libxc_functionals_end()
26 : !! if (new_ixc<0) call libxc_functionals_init(new_ixc,nspden)
27 : !! ! Compute XC stuff here.
28 : !! ! Revert libxc module to the original settings
29 : !! if (new_ixc<0) call libxc_functionals_end()
30 : !! if (old_ixc<0) call libxc_functionals_init(old_ixc,nspden)
31 : !!
32 : !! SOURCE
33 :
34 : #if defined HAVE_CONFIG_H
35 : #include "config.h"
36 : #endif
37 :
38 : #include "abi_common.h"
39 :
40 : MODULE m_kxc
41 :
42 : use defs_basis
43 : use m_abicore
44 : use m_errors
45 : use m_xmpi
46 : use m_crystal
47 : use m_distribfft
48 : use m_xcdata
49 : use libxc_functionals
50 : use m_dtset
51 :
52 : use defs_abitypes, only : MPI_type
53 : use m_fstrings, only : sjoin, itoa
54 : use m_numeric_tools, only : hermitianize
55 : use m_fft_mesh, only : g2ifft
56 : use m_fft, only : fourdp_6d, fourdp
57 : use m_mpinfo, only : initmpi_seq, destroy_mpi_enreg
58 : use m_spacepar, only : hartre
59 : use m_rhotoxc, only : rhotoxc
60 : use m_dfpt_mkvxc, only : dfpt_mkvxc
61 :
62 : implicit none
63 :
64 : private
65 : !!***
66 :
67 : public :: kxc_rpa ! Hartree kernel
68 : public :: kxc_local ! Compute local xc kernel in G space.
69 : public :: kxc_alda ! AL(S)DA kernel in reciprocal space, on the FFT grid.
70 : public :: kxc_pgg ! Compute the PGG-exchange kernel in reciprocal space (Phys. Rev. Lett. 76, 1212 (1996) [[cite:Petersilka1996]]).
71 : public :: kxc_eok ! linear or non-linear (ixceok = 2) energy optimized kernel of Dobson and Wang.
72 : public :: kxc_driver ! Driver routine (TODO)
73 : public :: kxc_ADA ! Adiabatic density approximation
74 :
75 :
76 : CONTAINS !=========================================================================================================================
77 : !!***
78 :
79 : !!****f* m_kxc/kxc_rpa
80 : !! NAME
81 : !! kxc_rpa
82 : !!
83 : !! FUNCTION
84 : !! Return the Hartree kernel:
85 : !! If option = 0, the bare Hartree kernel:
86 : !! krpa(ipw) = 4.0*pi/gsq(ipw) if gsq(ipw) /= 0.,
87 : !! krpa(ipw) = 0.0 if gsq(ipw) == 0. (1 <= ipw <= npw).
88 : !! If option /= 0, the Hartree kernel with a cut-off in real space beyond rcut_coulomb:
89 : !! krpa(ipw) = (4.0*pi/gsq(ipw))*(1.0-cos(sqrt(gsq(ipw))*rcut_coulomb)) if gsq(ipw) /= 0.,
90 : !! krpa(ipw) = 2.0*pi*rcut_coulomb**2 if gsq(ipw) == 0.
91 : !!
92 : !! INPUTS
93 : !! gsq(npw) = the squared norm of the planewaves.
94 : !! npw = number of planewaves in the gsq array.
95 : !! option = 0 for the bare Hartree kernel, /=0 for the cut-off Hartree kernel.
96 : !! rcut_coulomb = real space cut-off radius for the Coulomb interaction in Bohr.
97 : !!
98 : !! OUTPUT
99 : !! krpa(npw) = the Hartree kernel.
100 : !!
101 : !! SOURCE
102 :
103 0 : subroutine kxc_rpa(gsq,krpa,npw,option,rcut_coulomb)
104 :
105 : !Arguments -------------------------------------------------------------
106 : !scalars
107 : integer,intent(in) :: npw,option
108 : real(dp),intent(in) :: rcut_coulomb
109 : !arrays
110 : real(dp),intent(in) :: gsq(npw)
111 : real(dp),intent(out) :: krpa(npw)
112 :
113 : !Local variables -------------------------------------------------------
114 : !scalars
115 : integer :: ipw
116 : !***********************************************************************
117 :
118 0 : if (option == 0) then
119 : ! Compute the bare Hartree kernel.
120 0 : do ipw = 1,npw
121 0 : if (gsq(ipw) > tol12) then
122 0 : krpa(ipw) = four_pi/gsq(ipw)
123 : else
124 0 : krpa(ipw) = zero
125 : end if
126 : end do
127 :
128 : else
129 :
130 : ! Compute the Hartree kernel with a cut-off in real space beyond rcut_coulomb:
131 0 : do ipw = 1,npw
132 0 : if (gsq(ipw) > tol12) then
133 0 : krpa(ipw) = (four_pi/gsq(ipw))*(1._dp-cos(sqrt(gsq(ipw))*rcut_coulomb))
134 : else
135 0 : krpa(ipw) = two_pi*rcut_coulomb**2
136 : end if
137 : end do
138 :
139 : end if
140 :
141 0 : end subroutine kxc_rpa
142 : !!***
143 :
144 : !----------------------------------------------------------------------
145 :
146 : !!****f* m_kxc/kxc_local
147 : !! NAME
148 : !! kxc_local
149 : !!
150 : !! FUNCTION
151 : !! In a planewave basis set, the matrix of a local xc kernel:
152 : !!
153 : !! $f_{\rm xc}(\vec{r},\vec{r}') = f(\vec{r})\delta(\vec{r}-\vec{r}')$
154 : !!
155 : !! is just:
156 : !!
157 : !! $f_{\rm xc}(\vec{G},\vec{G}') = f(\vec{G}-\vec{G}')$.
158 : !!
159 : !! This subroutine calculates the matrix of such a local xc kernel given $f(\vec{G})$ on the FFT grid.
160 : !!
161 : !! INPUTS
162 : !! ispxc = 1 for the up-up spin channel.
163 : !! = 2 for the up-down (and down-up) spin channels.
164 : !! = 3 for the down-down spin channel.
165 : !! ispxc must be 1 if nspden = 1.
166 : !! kg_diel(3,npwdiel) = reduced planewave coordinates for the kxc matrix.
167 : !! kxcg(2,nfft) = $f(\vec{G})$ on the FFT grid.
168 : !! nfft = number of fft grid points.
169 : !! ngfft(1:3) = integer fft box dimensions, see getng for ngfft(4:8).
170 : !! npwdiel = number of planewaves for the susceptibility matrix.
171 : !! nspden = number of spin-density components.
172 : !! option = 0 do not compute the first row and column of the matrix of the
173 : !! xc kernel (which we assume to the G = 0 row and column).
174 : !! /= 0 compute the full matrix of the xc kernel.
175 : !!
176 : !! OUTPUT
177 : !! kxc(2,npwdiel,nspden,npwdiel,nspden) = the matrix of the xc kernel.
178 : !!
179 : !! SOURCE
180 :
181 0 : subroutine kxc_local(ispxc,kg_diel,kxc,kxcg,nfft,ngfft,npwdiel,nspden,option)
182 :
183 : !Arguments -------------------------------------------------------------
184 : !scalars
185 : integer,intent(in) :: ispxc,nfft,npwdiel,nspden,option
186 : !arrays
187 : integer,intent(in) :: kg_diel(3,npwdiel),ngfft(18)
188 : real(dp),intent(in) :: kxcg(2,nfft)
189 : real(dp),intent(out) :: kxc(2,npwdiel,nspden,npwdiel,nspden)
190 :
191 : !Local variables -------------------------------------------------------
192 : !For debugging purposes:
193 : !real(dp) :: c1,c2,c3
194 : !scalars
195 : integer :: i1,i2,i3,ifft,ipw1,ipw2,ipwstart,isp1,isp2,j1,j2,j3,k1,k2,k3,n1,n2
196 : integer :: n3
197 : logical :: ok
198 : character(len=500) :: msg
199 : !***********************************************************************
200 :
201 : !Check input parameters.
202 0 : if (nspden > 2) then
203 0 : ABI_ERROR('kxc_local does not work yet for nspden > 2.')
204 : end if
205 :
206 0 : isp1 = 1
207 0 : isp2 = 1
208 0 : ok = .true.
209 0 : if (nspden == 1) then
210 0 : select case (ispxc)
211 : case (1)
212 : isp1 = 1
213 : isp2 = 1
214 : case default
215 0 : ok = .false.
216 : end select
217 : else
218 0 : select case (ispxc)
219 : case (1)
220 : isp1 = 1
221 : isp2 = 1
222 : case (2)
223 : isp1 = 1
224 : isp2 = 2
225 : case (3)
226 : isp1 = 2
227 : isp2 = 2
228 : case default
229 0 : ok = .false.
230 : end select
231 : end if
232 :
233 : if (.not.ok) then
234 0 : write (msg,'(2(a,i0))')' The input ispxc = ',ispxc,' is not compatible with nspden = ',nspden
235 0 : ABI_BUG(msg)
236 : end if
237 :
238 0 : if (option == 0) then
239 0 : ipwstart = 2
240 0 : kxc(:,1,isp1,:,isp2) = 0._dp
241 0 : kxc(:,:,isp1,1,isp2) = 0._dp
242 : else
243 : ipwstart = 1
244 : end if
245 :
246 : ! Calculate the xc matrix.
247 0 : n1 = ngfft(1) ; n2 = ngfft(2) ; n3 = ngfft(3)
248 :
249 0 : do ipw2 = ipwstart,npwdiel
250 0 : j1 = kg_diel(1,ipw2) ; j2 = kg_diel(2,ipw2) ; j3 = kg_diel(3,ipw2)
251 :
252 : !Fill the diagonal.
253 :
254 0 : kxc(:,ipw2,isp1,ipw2,isp2) = kxcg(:,1)
255 :
256 : !Fill the off-diagonal elements.
257 :
258 0 : do ipw1 = ipw2+1,npwdiel
259 :
260 0 : i1 = kg_diel(1,ipw1) ; i2 = kg_diel(2,ipw1) ; i3 = kg_diel(3,ipw1)
261 :
262 : ! Compute the difference between G vectors.
263 : ! The use of two mod calls handles both i1-j1 >= n1 AND i1-j1 < 0.
264 :
265 0 : k1 = mod(n1+mod(i1-j1,n1),n1)
266 0 : k2 = mod(n2+mod(i2-j2,n2),n2)
267 0 : k3 = mod(n3+mod(i3-j3,n3),n3)
268 :
269 0 : ifft = k1+n1*(k2+n2*k3)+1
270 :
271 0 : kxc(1,ipw1,isp1,ipw2,isp2) = kxcg(1,ifft)
272 0 : kxc(2,ipw1,isp1,ipw2,isp2) = kxcg(2,ifft)
273 :
274 0 : kxc(1,ipw2,isp1,ipw1,isp2) = kxcg(1,ifft)
275 0 : kxc(2,ipw2,isp1,ipw1,isp2) = -kxcg(2,ifft)
276 : end do
277 : end do
278 :
279 : ! If needed, copy the up-down to the down-up spin channel.
280 0 : if (ispxc == 2) then
281 0 : do ipw2 = 1,npwdiel
282 0 : do ipw1 = 1,npwdiel
283 0 : kxc(1,ipw2,isp2,ipw1,isp1) = kxc(1,ipw1,isp1,ipw2,isp2)
284 0 : kxc(2,ipw2,isp2,ipw1,isp1) = -kxc(2,ipw1,isp1,ipw2,isp2)
285 : end do
286 : end do
287 : end if
288 :
289 : !DEBUG
290 : !See kxc_alda.f, "test kernel" DEBUG section.
291 : !do ipw2 = 1,npwdiel
292 : !j1 = kg_diel(1,ipw2) ; j2 = kg_diel(2,ipw2) ; j3 = kg_diel(3,ipw2)
293 : !do ipw1 = ipw2+1,npwdiel
294 : !i1 = kg_diel(1,ipw1) ; i2 = kg_diel(2,ipw1) ; i3 = kg_diel(3,ipw1)
295 : !k1 = mod(n1+mod(i1-j1,n1),n1)
296 : !k2 = mod(n2+mod(i2-j2,n2),n2)
297 : !k3 = mod(n3+mod(i3-j3,n3),n3)
298 : !ifft = k1+n1*(k2+n2*k3)+1
299 : !c1 = 0._dp ; c2 = 0._dp ; c3 = 0._dp
300 : !if (i1-j1 == 0) c1 = c1+0.0_dp
301 : !if (i2-j2 == 0) c2 = c2+0.0_dp
302 : !if (i3-j3 == 0) c3 = c3+0.0_dp
303 : !if (i1-j1 == 1) c1 = c1+0.5_dp
304 : !if (i2-j2 == 2) c2 = c2+0.5_dp
305 : !if (i3-j3 == 3) c3 = c3+0.5_dp
306 : !if (i1-j1 == -1) c1 = c1+0.5_dp
307 : !if (i2-j2 == -2) c2 = c2+0.5_dp
308 : !if (i3-j3 == -3) c3 = c3+0.5_dp
309 : !if ((abs(kxcg(1,ifft)-c1*c2*c3) > tol10).or.(abs(kxcg(2,ifft)) > tol10)) then
310 : !write (std_out,*) ' i1 i2 i3 ifft: ',i1,i2,i3,ifft
311 : !write (std_out,*) ' exp.: ',c1*c2*c3,' got: ',kxcg(:,ifft)
312 : !end if
313 : !end do
314 : !end do
315 : !ENDDEBUG
316 :
317 0 : end subroutine kxc_local
318 : !!***
319 :
320 : !----------------------------------------------------------------------
321 :
322 : !!****f* m_kxc/kxc_alda
323 : !! NAME
324 : !! kxc_alda
325 : !!
326 : !! FUNCTION
327 : !! If option = 1:
328 : !! Compute the AL(S)DA kernel in reciprocal space, on the FFT grid.
329 : !! If option = 2:
330 : !! Only computes the up-down channel of the AL(S)DA kernel, on the
331 : !! FFT grid, for use in the BPG kernel.
332 : !!
333 : !! INPUTS
334 : !! dtset <type(dataset_type)>=all input variables in this dataset
335 : !! ixc = choice of exchange-correlation functional.
336 : !! mpi_enreg=information about MPI parallelization
337 : !! nfft = number of fft grid points.
338 : !! ngfft(1:3) = integer fft box dimensions, see getng for ngfft(4:8).
339 : !! nspden = number of spin-density components.
340 : !! option = 1 compute the AL(S)DA kernel in reciprocal space.
341 : !! = 2 only computes the up-down channel of the AL(S)DA kernel,
342 : !! for use in the BPG kernel.
343 : !! rhor(nfft,nspden) = electron density in real space in electrons/bohr**3
344 : !! (total in first half and spin-up in second half if nspden = 2).
345 : !! rhocut = cut-off density for the local kernels (ALDA, EOK),
346 : !! relative to max(rhor(:,:)).
347 : !! rprimd(3,3) = dimensional primitive translations for real space in Bohr.
348 : !!
349 : !! OUTPUT
350 : !! kxcg(2,nfft,*) = the AL(S)DA kernel in reciprocal space, on the FFT grid
351 : !! (the third dimension is 2*nspden-1 if option = 1, and 1 if option = 2).
352 : !!
353 : !! WARNINGS
354 : !! Spin-polarized case not tested.
355 : !!
356 : !! SOURCE
357 :
358 0 : subroutine kxc_alda(dtset,ixc,kxcg,mpi_enreg,nfft,ngfft,nspden,option,rhor,rhocut,rprimd)
359 :
360 : !Arguments -------------------------------------------------------------
361 : !scalars
362 : integer,intent(in) :: ixc,nfft,nspden,option
363 : real(dp),intent(in) :: rhocut
364 : type(MPI_type),intent(in) :: mpi_enreg
365 : type(dataset_type),intent(in) :: dtset
366 : !arrays
367 : integer,intent(in) :: ngfft(18)
368 : real(dp),intent(in) :: rhor(nfft,2*nspden-1),rprimd(3,3)
369 : real(dp),intent(out) :: kxcg(2,nfft,*)
370 :
371 : !Local variables -------------------------------------------------------
372 : !No improved xc quadrature.
373 : !No core correction.
374 : !Dummy here.
375 : !For debugging purposes (see tests below):
376 : !integer :: i1,i2,i3,k1,n1,n2,n3
377 : !real(dp) :: kx,rho,rhomax,ftest
378 : !scalars
379 : integer :: ifft,ikxc,isp,n3xccc,ncut,nk3xc,nkxc,optionrhoxc,tim_fourdp
380 : logical :: non_magnetic_xc
381 : real(dp),parameter :: gsqcut=1._dp
382 : real(dp) :: el_temp,bigexc,bigsxc,rhocuttot,rhomin,vxcavg
383 : character(len=500) :: msg
384 : type(xcdata_type) :: xcdata
385 : !arrays
386 : real(dp) :: dum(0)
387 : real(dp),parameter :: dummyvgeo(3)=zero
388 0 : real(dp),allocatable :: kxcr(:,:),rhog(:,:),rhorcut(:,:),vhartree(:)
389 0 : real(dp),allocatable :: vxc(:,:),xccc3d(:)
390 : !***********************************************************************
391 :
392 : !For debugging purposes (see tests below):
393 : !ftest(i1,n1,k1) = 0._dp+1._dp*cos(k1*two_pi*float(i1)/float(n1))
394 :
395 : ! Check input parameters.
396 0 : if (nspden > 2) then
397 0 : ABI_ERROR('kxc_alda does not work yet for nspden > 2.')
398 : end if
399 :
400 : ! Allocate memory.
401 0 : ABI_MALLOC(rhorcut, (nfft,nspden))
402 0 : ABI_MALLOC(rhog, (2,nfft))
403 0 : ABI_MALLOC(vhartree, (nfft))
404 0 : ABI_MALLOC(vxc, (nfft,nspden))
405 :
406 0 : call xcdata_init(xcdata,dtset=dtset,intxc=0,ixc=ixc,nspden=nspden)
407 :
408 0 : non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
409 :
410 : ! Reinitialize the libxc module with the overridden values
411 0 : if (dtset%ixc<0) then
412 0 : call libxc_functionals_end()
413 : end if
414 :
415 0 : if (ixc<0) then
416 0 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
417 0 : call libxc_functionals_init(ixc,dtset%nspden,el_temp=el_temp,xc_tb09_c=dtset%xc_tb09_c)
418 : end if
419 :
420 : !to be adjusted for the call to rhotoxc
421 0 : nk3xc=1
422 :
423 : ! Cut-off the density.
424 0 : rhorcut(:,:) = rhor(:,:)
425 :
426 0 : do isp = 1,nspden
427 0 : rhomin = maxval(rhorcut(:,isp))*rhocut
428 :
429 0 : ncut = 0
430 0 : rhocuttot = 0._dp
431 :
432 0 : do ifft = 1,nfft
433 0 : if (rhorcut(ifft,isp) < rhomin) then
434 0 : ncut = ncut+1
435 0 : rhocuttot = rhocuttot+rhorcut(ifft,isp)
436 0 : rhorcut(ifft,isp) = rhomin
437 : end if
438 : end do
439 :
440 0 : if (ncut > 0) then
441 : write (msg,'(a,es10.3,3a,i1,a,i6,a,f6.3,3a,f6.3,a)') &
442 0 : 'rhocut = ',rhocut,'.',ch10,&
443 0 : 'For isp = ',isp,' the density was cut-off at ',ncut,' (',100._dp*float(ncut)/float(ifft),'%) grid points.',ch10,&
444 0 : 'These points account for ',100._dp*rhocuttot/sum(rhor(:,isp)),'% of the total density.'
445 0 : ABI_WARNING(msg)
446 : end if
447 :
448 : end do
449 :
450 : ! Calculate the AL(S)DA kernel in real space.
451 0 : rhog(:,:) = zero !We do not need the Hartree potential.
452 0 : tim_fourdp=0
453 :
454 0 : if ((option == 1).or.((option == 2).and.(nspden == 2))) then
455 :
456 0 : nkxc = 2*nspden-1
457 0 : n3xccc=0
458 0 : ABI_MALLOC(kxcr,(nfft,nkxc))
459 0 : ABI_MALLOC(xccc3d,(n3xccc))
460 :
461 0 : optionrhoxc = 2 !See rhotoxc.f
462 :
463 0 : call hartre(1,gsqcut,3,0,mpi_enreg,nfft,ngfft,1,zero,rhog,rprimd,dummyvgeo,vhartree)
464 : call rhotoxc(bigexc,bigsxc,kxcr,mpi_enreg,nfft,ngfft,dum,0,dum,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
465 0 : optionrhoxc,rhorcut,rprimd,1,vxc,vxcavg,xccc3d,xcdata,vhartr=vhartree)
466 :
467 : ! DEBUG
468 : ! fx for tests.
469 : ! write (std_out,'(a)') ' kxc_alda: Using exchange-only kernel for tests.'
470 : ! rhomin = minval(rhor(:,1))
471 : ! rhomax = maxval(rhor(:,1))
472 : ! write (std_out,'(a,es12.5,a,es12.5)') ' kxc_alda: rhomin = ',rhomin,' rhomax = ',rhomax
473 : ! write (std_out,'(a)') ' kxc_alda: loping below 0.2*rhomax.'
474 : ! kx = (3._dp/4._dp)*((3._dp/pi)**(1._dp/3._dp))
475 : ! do ifft = 1,nfft
476 : ! rho = rhor(ifft,1)
477 : ! rho = max(rho,0.2_dp*rhomax)
478 : ! kxcr(ifft,1) = -(4._dp/9._dp)*kx*(rho**(-2._dp/3._dp))
479 : ! write (std_out,'(i4,a,es12.5)') ifft,': ',kxcr(ifft,1)
480 : ! end do
481 : ! write (std_out,'(a,es12.5)') 'kxcrmin: ',minval(kxcr(:,1))
482 : ! write (std_out,'(a,es12.5)') 'kxcrmax: ',maxval(kxcr(:,1))
483 : ! ENDDEBUG
484 :
485 : ! DEBUG
486 : ! test kernel.
487 : ! write(std_out,'(a)') ' kxc_alda: Using test kernel for tests.'
488 : ! n1 = ngfft(1) ; n2 = ngfft(2) ; n3 = ngfft(3)
489 : ! do i3 = 0,n3-1
490 : ! do i2 = 0,n2-1
491 : ! do i1 = 0,n1-1
492 : ! ifft = i1+n1*(i2+n2*i3)+1
493 : ! kxcr(ifft,1) = ftest(i1,n1,1)*ftest(i2,n2,2)*ftest(i3,n3,3)
494 : ! end do
495 : ! end do
496 : ! end do
497 : ! ENDDEBUG
498 :
499 : ! Calculate the Fourier transform of the AL(S)DA kernel.
500 :
501 0 : if (option == 1) then
502 0 : do ikxc = 1,nkxc
503 0 : call fourdp(1,kxcg(:,:,ikxc),kxcr(:,ikxc),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
504 : end do
505 : else
506 0 : call fourdp(1,kxcg(:,:,1),kxcr(:,2),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
507 : end if
508 :
509 0 : else if ((option == 2).and.(nspden == 1)) then
510 :
511 0 : nkxc = 2
512 0 : n3xccc=0
513 0 : ABI_MALLOC(kxcr,(nfft,nkxc))
514 0 : ABI_MALLOC(xccc3d,(n3xccc))
515 :
516 0 : optionrhoxc = -2 !See rhotoxc.f
517 :
518 0 : call hartre(1,gsqcut,3,0,mpi_enreg,nfft,ngfft,1,zero,rhog,rprimd,dummyvgeo,vhartree)
519 : call rhotoxc(bigexc,bigsxc,kxcr,mpi_enreg,nfft,ngfft,dum,0,dum,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
520 0 : & optionrhoxc,rhorcut,rprimd,1,vxc,vxcavg,xccc3d,xcdata,vhartr=vhartree)
521 :
522 0 : kxcr(:,2) = 0.5_dp*kxcr(:,2)
523 :
524 : ! Calculate the Fourier transform of the up-down channel of the AL(S)DA kernel.
525 0 : call fourdp(1,kxcg(:,:,1),kxcr(:,2),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
526 :
527 : else
528 0 : ABI_ERROR(sjoin("Invalid option:", itoa(option)))
529 : end if
530 :
531 : !write(std_out,*)' kxc_alda: Exc = ',bigexc
532 : !write(std_out,*)' kxc_alda: Sxc = ',bigsxc
533 : !write(std_out,*)' kxc_alda: <Vxc> = ',vxcavg
534 :
535 : ! Revert libxc module to the original settings
536 0 : if (ixc<0) then
537 0 : call libxc_functionals_end()
538 : end if
539 0 : if (dtset%ixc<0) then
540 0 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
541 0 : call libxc_functionals_init(dtset%ixc,dtset%nspden,el_temp=el_temp,xc_tb09_c=dtset%xc_tb09_c)
542 : end if
543 :
544 : ! Free memory.
545 0 : ABI_FREE(rhorcut)
546 0 : ABI_FREE(rhog)
547 0 : ABI_FREE(vhartree)
548 0 : ABI_FREE(vxc)
549 0 : ABI_FREE(kxcr)
550 0 : ABI_FREE(xccc3d)
551 :
552 0 : end subroutine kxc_alda
553 : !!***
554 :
555 : !----------------------------------------------------------------------
556 :
557 : !!****f* m_kxc/kxc_pgg
558 : !! NAME
559 : !! kxc_pgg
560 : !!
561 : !! FUNCTION
562 : !! Compute the PGG-exchange kernel in reciprocal space
563 : !! (Phys. Rev. Lett. 76, 1212 (1996) [[cite:Petersilka1996]]).
564 : !!
565 : !! INPUTS
566 : !! gmet=reciprocal space metric (bohr**-2)
567 : !! npw=number of plane waves
568 : !! rcut_coulomb=real space cutoff radius for Coulomb interaction (bohr)
569 : !! susmat(2,npw,npw)=density weighted squared density matrix in reciprocal space
570 : !! ucvol=unit cell volume (bohr**3)
571 : !!
572 : !! OUTPUT
573 : !! khxcg(2,npwdiel,nspden,npwdiel,nspden)=PGG-exhange kernel in G space, at
574 : !! full interaction strength
575 : !!
576 : !! NOTES
577 : !! The density weighted squared density matrix (actually the reduced=summed-over-spin
578 : !! density matrix) is convolved with the spherically cutoff Coulomb interaction.
579 : !!
580 : !! WARNINGS
581 : !! a - 'rcut_coulomb' should be chosen consistently with cutoffs elsewhere,
582 : !! for instance dieltcel8.f
583 : !! b - applicable for spin-polarized case as well, through input 'susmat',
584 : !! but this has not been checked
585 : !!
586 : !! TODO
587 : !! If simply the squared density matrix is input through 'susmat' the
588 : !! exchange energy is recovered as the zero-G component of the resulting 'khxcg'
589 : !! (then not the kernel of course). This could help to check convergence
590 : !! with respect to 'npw'. See +ex_pgg comment.
591 : !!
592 : !! SOURCE
593 :
594 0 : subroutine kxc_pgg(gmet,kg,khxcg,npw,rcut_coulomb,susmat,ucvol)
595 :
596 : !Arguments ------------------------------------
597 : !scalars
598 : integer,intent(in) :: npw
599 : real(dp),intent(in) :: rcut_coulomb,ucvol
600 : !arrays
601 : integer,intent(in) :: kg(3,npw)
602 : real(dp),intent(in) :: gmet(3,3),susmat(2,npw,npw)
603 : real(dp),intent(out) :: khxcg(2,npw,npw)
604 :
605 : !Local variables-------------------------------
606 : !scalars
607 : integer :: i1,i2,i3,ig,ii,ikg11,ikg12,ikg13,ikg21,ikg22,ikg23,ipw1,ipw2
608 : integer :: j1,j2,j3,jg,jj
609 : real(dp),parameter :: diffgsq=1.d-2
610 : real(dp) :: kg_red1,kg_red2,kg_red3,gsquar,tpisq
611 : !arrays
612 : integer :: kgmax(3)
613 0 : integer,allocatable :: index_g_inv(:,:,:),jgarr(:)
614 0 : real(dp),allocatable :: gsq(:),sumg(:),vcoul(:)
615 : ! *************************************************************************
616 :
617 : !write(std_out,*) '%kxc_pgg: enter', 'npw=',npw
618 :
619 : !tpisq is (2 Pi) **2:
620 0 : tpisq=(two_pi)**2
621 :
622 0 : kgmax(:)=0
623 0 : do ipw1=1,npw
624 0 : do jj=1,3
625 0 : kgmax(jj)=max( kg(jj,ipw1), kgmax(jj) )
626 : end do
627 : end do
628 : !write(std_out,*) 'kgmax:',kgmax(1:3)
629 :
630 : ! Perform allocations
631 0 : ABI_MALLOC(index_g_inv,(-2*kgmax(1):2*kgmax(1),-2*kgmax(2):2*kgmax(2),-2*kgmax(3):2*kgmax(3)))
632 0 : ABI_MALLOC(jgarr,(npw))
633 0 : ABI_MALLOC(gsq,(npw))
634 0 : ABI_MALLOC(sumg,(2))
635 0 : ABI_MALLOC(vcoul,(npw))
636 :
637 : !write(std_out,*) '%kxc_pg: creating plane wave index and coulomb potential'
638 0 : index_g_inv(:,:,:)=0
639 0 : do ipw1=1,npw
640 0 : index_g_inv(kg(1,ipw1),kg(2,ipw1),kg(3,ipw1))=ipw1
641 :
642 : !write(std_out,'(i5,2x,3i3,2x,i4)') ipw1,kg(1,ipw1),kg(2,ipw1),kg(3,ipw1)
643 :
644 0 : kg_red1=dble(kg(1,ipw1))
645 0 : kg_red2=dble(kg(2,ipw1))
646 0 : kg_red3=dble(kg(3,ipw1))
647 : gsquar=tpisq*(gmet(1,1)*kg_red1**2+gmet(2,2)*kg_red2**2+gmet(3,3)*kg_red3**2 &
648 : & +2.0_dp*( (gmet(1,2)*kg_red2+gmet(1,3)*kg_red3)* kg_red1 + &
649 0 : & gmet(2,3)*kg_red2*kg_red3) )
650 : ! Distinguish G=0 from other elements
651 0 : if(gsquar > 1.0d-12)then
652 0 : vcoul(ipw1)=four_pi/gsquar*(1._dp-cos(sqrt(gsquar)*rcut_coulomb))
653 : else
654 0 : vcoul(ipw1)=four_pi*0.5_dp*rcut_coulomb**2
655 : end if
656 :
657 : end do
658 :
659 : !write(std_out,*) '%kxc_pg: starting convolution integral'
660 : !loop over G1,G2 components of the density matrix
661 0 : do ipw2=1,npw
662 0 : ikg21=kg(1,ipw2)
663 0 : ikg22=kg(2,ipw2)
664 0 : ikg23=kg(3,ipw2)
665 :
666 0 : do ii=1,npw
667 0 : j1=ikg21-kg(1,ii)
668 0 : j2=ikg22-kg(2,ii)
669 0 : j3=ikg23-kg(3,ii)
670 0 : jgarr(ii)=index_g_inv(j1,j2,j3)
671 : end do
672 :
673 0 : do ipw1=1,ipw2
674 0 : ikg11=kg(1,ipw1)
675 0 : ikg12=kg(2,ipw1)
676 0 : ikg13=kg(3,ipw1)
677 :
678 : ! do the convolution integral
679 0 : sumg(:)=0._dp
680 0 : do ii=1,npw
681 :
682 0 : if( jgarr(ii) /= 0 ) then
683 0 : i1=ikg11-kg(1,ii)
684 0 : i2=ikg12-kg(2,ii)
685 0 : i3=ikg13-kg(3,ii)
686 :
687 : ! j1=ikg21-kg(1,ii)
688 : ! j2=ikg22-kg(2,ii)
689 : ! j3=ikg23-kg(3,ii)
690 :
691 0 : ig=index_g_inv(i1,i2,i3)
692 : ! jg=index_g_inv(j1,j2,j3)
693 :
694 0 : if( ig /= 0 ) then
695 0 : jg=jgarr(ii)
696 :
697 : ! write(std_out,'(i5,2x,3i3,1x,3i3,2x,2i4)') ii,i1,i2,i3,&
698 : ! & kg(1,jg),kg(2,jg),kg(3,jg),&
699 : ! & ig,jg
700 :
701 0 : sumg(1)=sumg(1)+susmat(1,ig,jg)*vcoul(ii)
702 0 : sumg(2)=sumg(2)+susmat(2,ig,jg)*vcoul(ii)
703 :
704 : end if
705 : end if
706 :
707 : end do
708 0 : khxcg(:,ipw1,ipw2)=-sumg(:)*ucvol
709 :
710 : ! if(ipw1==ipw2) write(std_out,'(2i4,2(1x,es14.6))') ipw1,ipw2,khxcg(1,ipw1,ipw1),vcoul(ipw1)
711 : ! write(std_out,'(2i4,3(1x,es14.6))') ipw1,ipw2,khxcg(1:2,ipw1,ipw2),vcoul(ipw1)
712 :
713 : end do
714 : end do
715 :
716 : !verify hermiticity, note: ipw1 loop above must end at npw
717 : !write(std_out,*) '%kxc_pgg: check hermiticity of pgg kernel'
718 : !do ipw2=1,npw,max(2,npw/10)
719 : !do ipw1=ipw2,npw,max(2,npw/10)
720 : !write(std_out,'(2i4,2(1x,es14.6))') ipw1,ipw2,&
721 : !& khxcg(1,ipw1,ipw2)-khxcg(1,ipw2,ipw1),&
722 : !& khxcg(2,ipw1,ipw2)+khxcg(2,ipw2,ipw1)
723 : !end do
724 : !end do
725 :
726 : ! Impose hermiticity
727 0 : write(std_out,*) '%kxc_pg: imposing hermiticity'
728 0 : do ipw2=1,npw
729 0 : do ipw1=ipw2+1,npw
730 0 : khxcg(1,ipw1,ipw2)= khxcg(1,ipw2,ipw1)
731 0 : khxcg(2,ipw1,ipw2)=-khxcg(2,ipw2,ipw1)
732 : end do
733 : end do
734 :
735 : !write(std_out,'(a10,2(1x,es20.12))') '+ex_pgg? ', 0.5_dp*khxcg(1,1,1)/ucvol
736 :
737 0 : ABI_FREE(index_g_inv)
738 0 : ABI_FREE(jgarr)
739 0 : ABI_FREE(gsq)
740 0 : ABI_FREE(sumg)
741 0 : ABI_FREE(vcoul)
742 :
743 : !write(std_out,*) '%kxc_pgg: done'
744 :
745 0 : end subroutine kxc_pgg
746 : !!***
747 :
748 : !----------------------------------------------------------------------
749 :
750 : !!****f* m_kxc/kxc_eok
751 : !! NAME
752 : !! kxc_eok
753 : !!
754 : !! FUNCTION
755 : !! Compute the linear (ixceok = 1) or non-linear (ixceok = 2)
756 : !! energy optimized kernel of Dobson and Wang, in reciprocal space, on the FFT grid.
757 : !! See J. Dobson and J. Wang, Phys. Rev. B 62, 10038 (2000) [[cite:Dobson2000]].
758 : !!
759 : !! INPUTS
760 : !! ixceok = 1 linear energy optimized kernel.
761 : !! = 2 non-linear energy optimized kernel.
762 : !! mpi_enreg=information about MPI parallelization
763 : !! nfft = number of fft grid points.
764 : !! ngfft(1:3) = integer fft box dimensions, see getng for ngfft(4:8).
765 : !! nspden = number of spin-density components.
766 : !! rhor(nfft,nspden) = electron density in real space in electrons/bohr**3
767 : !! (total in first half and spin-up in second half if nspden = 2).
768 : !! rhocut = cut-off density for the local kernels (ALDA, EOK),
769 : !! relative to max(rhor(:,:)).
770 : !! OUTPUT
771 : !! kxcg(2,nfft,2*nspden-1) = the EOK kernel in reciprocal space, on the FFT grid.
772 : !!
773 : !! SOURCE
774 :
775 0 : subroutine kxc_eok(ixceok,kxcg,mpi_enreg,nfft,ngfft,nspden,rhor,rhocut)
776 :
777 : !Arguments -------------------------------------------------------------
778 : !scalars
779 : integer,intent(in) :: ixceok,nfft,nspden
780 : real(dp),intent(in) :: rhocut
781 : type(MPI_type),intent(in) :: mpi_enreg
782 : !arrays
783 : integer,intent(in) :: ngfft(18)
784 : real(dp),intent(in) :: rhor(nfft,2*nspden-1)
785 : real(dp),intent(out) :: kxcg(2,nfft,2*nspden-1)
786 :
787 : !Local variables -------------------------------------------------------
788 : !Maximum value allowed for rs.
789 : !scalars
790 : integer :: ifft,ikxc,ncut,nkxc,nlop,tim_fourdp
791 : real(dp),parameter :: rslim=50._dp,dummyvgeo(3)=zero
792 : real(dp) :: a2,a3,a4,rho,rhocuttot,rhomin,rs
793 : character(len=500) :: msg
794 : !arrays
795 0 : real(dp),allocatable :: kxcr(:,:)
796 : !***********************************************************************
797 :
798 : ! Check input parameters.
799 0 : if (nspden > 1) then
800 0 : ABI_ERROR('kxc_eok does not work yet for nspden > 1.')
801 : end if
802 :
803 : ! Values of a2, a3 and a4 for case 1
804 0 : a2 = 0.0_dp
805 0 : a3 = 0.0_dp
806 0 : a4 = 0.0_dp
807 :
808 0 : select case (ixceok)
809 : case (1)
810 : a2 = -0.51887_dp
811 : a3 = 4.9359d-03
812 0 : a4 = -5.9603d-05
813 : case (2)
814 0 : a2 = -0.50044_dp
815 0 : a3 = 4.9653d-03
816 0 : a4 = -3.3660d-05
817 : case default
818 0 : ABI_ERROR(' kxc_eok: ixceok /= 1 (linear EOK) or 2 (non-linear EOK).')
819 : end select
820 :
821 : !Allocate memory.
822 0 : nkxc = 2*nspden-1
823 :
824 0 : ABI_MALLOC(kxcr,(nfft,nkxc))
825 :
826 : ! Calculate the energy optimized kernel in real space.
827 0 : nlop = 0
828 :
829 0 : rhomin = rhocut*maxval(rhor(:,:))
830 :
831 0 : ncut = 0
832 0 : rhocuttot = 0._dp
833 :
834 0 : do ifft = 1,nfft
835 0 : rho = rhor(ifft,1)
836 :
837 0 : if (rho < rhomin) then
838 0 : ncut = ncut+1
839 0 : rhocuttot = rhocuttot+rho
840 0 : rho = rhomin
841 : end if
842 :
843 0 : rs = (3._dp/(4._dp*pi*rho))**(1._dp/3._dp)
844 :
845 0 : if (rs > rslim) then
846 0 : rs = rslim
847 0 : nlop = nlop+1
848 : end if
849 :
850 0 : kxcr(ifft,1) = a2*rs**2+a3*rs**3+a4*rs**4
851 : end do
852 :
853 0 : if (ncut > 0) then
854 : write (msg,'(a,es10.3,3a,i1,a,i6,a,f6.3,3a,f6.3,a)') &
855 0 : 'rhocut = ',rhocut,'.',ch10,&
856 0 : 'For isp = ',1,' the density was cut-off at ',ncut,' (',100._dp*float(ncut)/float(ifft),'%) grid points.',ch10,&
857 0 : 'These points account for ',100._dp*rhocuttot/sum(rhor(:,1)),'% of the total density.'
858 0 : ABI_WARNING(msg)
859 : end if
860 :
861 0 : if (nlop > 0) then
862 : write (msg,'(a,f6.2,a,i6,a,f6.3,a)') &
863 0 : 'rs still exceeds ',rslim,' Bohr at ',nlop,' (',100._dp*float(nlop)/float(ifft),'%) grid points (after cut-off).'
864 0 : ABI_WARNING(msg)
865 : end if
866 :
867 : ! Calculate the Fourier transform of the energy optimized kernel.
868 0 : tim_fourdp=0
869 0 : do ikxc = 1,nkxc
870 0 : call fourdp(1,kxcg(:,:,ikxc),kxcr(:,ikxc),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
871 : end do
872 :
873 0 : ABI_FREE(kxcr)
874 :
875 0 : end subroutine kxc_eok
876 : !!***
877 :
878 : !----------------------------------------------------------------------
879 :
880 : !!****f* m_kxc/kxc_driver
881 : !! NAME
882 : !! kxc_driver
883 : !!
884 : !! FUNCTION
885 : !! Calculate the exchange-correlation kernel in reciprocal space.
886 : !! Require density in real space on the FFT mesh. MPI-FFT is not supported.
887 : !!
888 : !! INPUTS
889 : !! Dtset<dataset_type>=all input variables in this dataset
890 : !! Cryst<crystal_t>=Info on the crystal structure.
891 : !! ixc = choice for the exchange-correlation potential.
892 : !! ngfft(18)=contain all needed information about 3D FFT,
893 : !! see ~abinit/doc/variables/vargs.htm#ngfft
894 : !! nfft_tot = Total number of points on the FFT grid.
895 : !! nspden=Number of independent spin densities.
896 : !! rhor(nfft_tot,nspden) = the charge density on the full FFT grid.
897 : !! (total in first half and spin-up in second half if nspden=2)
898 : !! npw: the size of kernel matrix
899 : !! dim_kxcg=dimension of the kernel.
900 : !! comm=MPI communicator.
901 : !! [dbg_mode]=Optional flag used to switch on the debug mode.
902 : !!
903 : !! OUTPUT
904 : !! FIXME: Why are we using nfft_tot instead of the G-sphere
905 : !! kxcg(nfft_tot,dim_kxcg) = the exchange-correlation potential on the FFT grid.
906 : !! warning: the kernel is not divided by the unit cell volume
907 : !!
908 : !! NOTES
909 : !! No xc quadrature
910 : !! No nl core correction
911 : !!
912 : !! SOURCE
913 :
914 11 : subroutine kxc_driver(Dtset,Cryst,ixc,ngfft,nfft_tot,nspden,rhor,npw,dim_kxcg,kxcg,gvec,comm,dbg_mode)
915 :
916 : !Arguments ------------------------------------
917 : !scalars
918 : integer,intent(in) :: ixc,npw,nfft_tot,nspden,dim_kxcg,comm
919 : logical,optional,intent(in) :: dbg_mode
920 : type(crystal_t),intent(in) :: Cryst
921 : type(dataset_type),intent(in) :: Dtset
922 : !arrays
923 : integer,intent(in) :: gvec(3,npw),ngfft(18)
924 : real(dp),intent(in) :: rhor(nfft_tot,nspden)
925 : complex(gwp),intent(out) :: kxcg(nfft_tot,dim_kxcg)
926 :
927 : !Local variables ------------------------------
928 : !scalars
929 : integer :: cplex,i1,i2,i3,ig,igp,iq,ir,n3xccc,ngfft1,ngfft2,izero
930 : integer :: ngfft3,nkxc,option,ikxc,nk3xc,my_rank,master
931 : logical :: non_magnetic_xc
932 : real(dp) :: el_temp,bigexc,bigsxc,expo,gpqx,gpqy,gpqz,gsqcut,vxcavg
933 : character(len=500) :: fname ! msg,
934 : type(xcdata_type) :: xcdata
935 11 : type(MPI_type) :: MPI_enreg_seq
936 : !arrays
937 : real(dp) :: qphon(3),dum(0)
938 : real(dp),parameter :: dummyvgeo(3)=zero
939 11 : real(dp),allocatable :: kxcpw_g(:,:),kxcr(:,:),phas(:,:,:)
940 11 : real(dp),allocatable :: rhog(:,:),vhartr(:),kxcpw_r(:,:),vxclda(:,:)
941 11 : real(dp),allocatable :: xccc3d(:),xx(:,:), my_kxcg(:,:)
942 : !************************************************************************
943 :
944 11 : ABI_CHECK_IEQ(Dtset%nsppol, 1, 'nsppol/=1 not coded')
945 11 : ABI_CHECK_IEQ(Dtset%nspden, 1, 'nspden /= 1 not coded')
946 44 : ABI_CHECK_IEQ(nfft_tot, PRODUCT(ngfft(1:3)), "mismatch in nfftot")
947 :
948 : ! Fake MPI_type for the sequential part.
949 11 : call initmpi_seq(MPI_enreg_seq)
950 11 : call MPI_enreg_seq%distribfft%init_seq('c',ngfft(2),ngfft(3),'all')
951 11 : my_rank = xmpi_comm_rank(comm)
952 11 : master =0
953 :
954 11 : call wrtout(std_out,sjoin(' kxc_driver: calculating exchange-correlation kernel using ixc: ', itoa(ixc)))
955 11 : call xcdata_init(xcdata,dtset=Dtset,intxc=0,ixc=ixc,nspden=nspden)
956 :
957 11 : if (all(xcdata%xclevel /= [1,2])) then
958 0 : ABI_ERROR(sjoin("Unsupported xclevel: ", itoa(xcdata%xclevel)))
959 : end if
960 :
961 11 : ngfft1=ngfft(1)
962 11 : ngfft2=ngfft(2)
963 11 : ngfft3=ngfft(3)
964 :
965 11 : non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
966 :
967 11 : if (ixc>=1.and.ixc<11) then ! LDA case
968 11 : nkxc= 2*min(nspden,2)-1 ! 1 or 3
969 0 : elseif (ixc==51) then ! TLDA case, same as LDA above
970 0 : nkxc= 2*min(nspden,2)-1
971 : else ! GGA case
972 0 : nkxc=12*min(nspden,2)-5 ! 7 or 19
973 0 : ABI_CHECK_IEQ(dtset%xclevel, 2,"Functional should be GGA")
974 0 : ABI_ERROR("GGA functional not tested")
975 : end if
976 :
977 44 : ABI_MALLOC(kxcr, (nfft_tot,nkxc))
978 :
979 : ! gsqcut and rhog are zeroed because they are not used by rhotoxc if 1<=ixc<=16 and option=0
980 11 : gsqcut=zero
981 :
982 33 : ABI_MALLOC(rhog,(2,nfft_tot))
983 33 : ABI_MALLOC(vhartr,(nfft_tot))
984 236774 : rhog(:,:)=zero
985 : !MG FIXME this is the 3D core electron density for XC core correction (bohr^-3)
986 : !should implement the non linear core correction
987 11 : n3xccc=0
988 11 : ABI_MALLOC(xccc3d,(n3xccc))
989 44 : ABI_MALLOC(vxclda,(nfft_tot,nspden))
990 :
991 11 : option=2 ! 2 for Hxc and kxcr (no paramagnetic part if nspden=1)
992 11 : qphon =zero
993 :
994 : ! to be adjusted for the call to rhotoxc
995 11 : nk3xc=1
996 11 : izero=0
997 :
998 : ! Reinitialize the libxc module with the overridden values
999 11 : if (dtset%ixc<0) then
1000 0 : call libxc_functionals_end()
1001 : end if
1002 11 : if (ixc<0) then
1003 0 : el_temp=merge(Dtset%tphysel,Dtset%tsmear,Dtset%tphysel>tol8.and.Dtset%occopt/=3.and.Dtset%occopt/=9)
1004 0 : call libxc_functionals_init(ixc,Dtset%nspden,el_temp=el_temp,xc_tb09_c=Dtset%xc_tb09_c)
1005 : end if
1006 :
1007 11 : call hartre(1,gsqcut,3,izero,MPI_enreg_seq,nfft_tot,ngfft,1,zero,rhog,Cryst%rprimd,dummyvgeo,vhartr)
1008 :
1009 : ! Compute the XC kernel.
1010 : call rhotoxc(bigexc,bigsxc,kxcr,MPI_enreg_seq,nfft_tot,ngfft,&
1011 : dum,0,dum,0,nkxc,nk3xc,non_magnetic_xc,&
1012 : n3xccc,option,rhor,Cryst%rprimd,&
1013 11 : 1,vxclda,vxcavg,xccc3d,xcdata,vhartr=vhartr)
1014 :
1015 11 : ABI_FREE(rhog)
1016 11 : ABI_FREE(vhartr)
1017 :
1018 : ! print Kxc
1019 11 : if (present(dbg_mode)) then
1020 1 : if (dbg_mode .and. my_rank==master) then
1021 0 : fname = 'xc_Kxc.xsf'
1022 0 : call cryst%write_xsf_data(fname, ngfft1, ngfft2, ngfft3, kxcr(:,1))
1023 : end if
1024 :
1025 : end if
1026 :
1027 11 : ABI_FREE(xccc3d)
1028 11 : ABI_FREE(vxclda)
1029 :
1030 22 : ABI_MALLOC(my_kxcg,(2,nfft_tot))
1031 22 : do ikxc=1,nkxc
1032 11 : call fourdp(1,my_kxcg,kxcr(:,ikxc),-1,MPI_enreg_seq,nfft_tot,1,ngfft,0)
1033 78943 : kxcg(:,ikxc)=CMPLX(my_kxcg(1,:),my_kxcg(2,:))
1034 : end do
1035 :
1036 : !write(std_out,*)"kxcr(r=0)",kxcr(1,1)
1037 : !write(std_out,*)"my_kxg(G=0)",my_kxcg(:,1)
1038 : !write(std_out,*)"SUM kxcr/nfft_tot ",SUM(kxcr(:,1))/nfft_tot
1039 : !write(std_out,*)"SUM my_kxg ",SUM(kxcg(:,1))
1040 :
1041 11 : ABI_FREE(my_kxcg)
1042 :
1043 : !MG this part is never executed, but one should use dfpt_mkvxc for the GGA kernel.
1044 11 : if (xcdata%xclevel==2) then
1045 0 : ABI_ERROR("check GGA implementation")
1046 0 : cplex=2
1047 0 : ABI_MALLOC(phas,(cplex*nfft_tot,npw,nspden))
1048 0 : ABI_MALLOC(kxcpw_r,(cplex*nfft_tot,nspden))
1049 0 : ABI_MALLOC(xx,(3,nfft_tot))
1050 0 : ABI_MALLOC(kxcpw_g,(2,nfft_tot))
1051 :
1052 0 : kxcg = czero
1053 :
1054 : ! find the coordinates for all r in the FFT grid
1055 : ir=0
1056 0 : do i3=1,ngfft3
1057 0 : do i2=1,ngfft2
1058 0 : do i1=1,ngfft1
1059 0 : ir=ir+1
1060 0 : xx(1,ir)=dble((i1-1))/ngfft1
1061 0 : xx(2,ir)=dble((i2-1))/ngfft2
1062 0 : xx(3,ir)=dble((i3-1))/ngfft3
1063 : end do
1064 : end do
1065 : end do
1066 :
1067 0 : do iq=1,1
1068 :
1069 : ! Calculate at once exp(i(G+q).r), for all possible q,G,r
1070 0 : do ig=1,npw
1071 0 : gpqx=dble(gvec(1,ig))
1072 0 : gpqy=dble(gvec(2,ig))
1073 0 : gpqz=dble(gvec(3,ig))
1074 0 : do ir=1,nfft_tot
1075 0 : expo=gpqx*xx(1,ir)+gpqy*xx(2,ir)+gpqz*xx(3,ir)
1076 0 : phas(2*ir-1,ig,1)= cos(two_pi*expo)
1077 0 : phas(2*ir,ig,1) = sin(two_pi*expo)
1078 : end do
1079 : end do
1080 :
1081 : ! Calculate $K(G,G'',q)=\frac{1}{nfft_tot}\sum_{r} exp(-i(q+G_{2}).r_{2} kxcr(r_{1}r_{2}) exp(i(q+G_{1}).r_{1} $
1082 0 : do igp=1,npw
1083 0 : kxcpw_r(:,:)=zero
1084 :
1085 : call dfpt_mkvxc(cplex,ixc,kxcr,MPI_enreg_seq,nfft_tot,ngfft,dum,0,dum,0,nkxc,non_magnetic_xc,&
1086 0 : nspden,n3xccc,option,qphon(:),phas(:,igp,:),Cryst%rprimd,1,kxcpw_r,xccc3d)
1087 :
1088 : ! FFT the first index to --> to G space
1089 0 : call fourdp(cplex,kxcpw_g(:,:),kxcpw_r(:,1),-1,MPI_enreg_seq,nfft_tot,1,ngfft,0)
1090 :
1091 : !kxcg(:,igp,iq)=CMPLX(kxcpw_g(1,igfft(:)),kxcpw_g(2,igfft(:)))
1092 : !kxcg(:,igp)=CMPLX(kxcpw_g(1,igfft(:)),kxcpw_g(2,igfft(:)))
1093 : end do ! igp
1094 : end do ! iq
1095 :
1096 0 : ABI_FREE(phas)
1097 0 : ABI_FREE(kxcpw_r)
1098 0 : ABI_FREE(xx)
1099 0 : ABI_FREE(kxcpw_g)
1100 : end if !xclevel==2
1101 :
1102 : ! Revert libxc module to the original settings
1103 11 : if (ixc<0) then
1104 0 : call libxc_functionals_end()
1105 : end if
1106 11 : if (dtset%ixc<0) then
1107 0 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
1108 0 : call libxc_functionals_init(dtset%ixc,dtset%nspden,el_temp=el_temp,xc_tb09_c=dtset%xc_tb09_c)
1109 : end if
1110 :
1111 11 : call destroy_mpi_enreg(MPI_enreg_seq)
1112 11 : ABI_FREE(kxcr)
1113 :
1114 11 : end subroutine kxc_driver
1115 : !!***
1116 :
1117 : !----------------------------------------------------------------------
1118 :
1119 : !!****f* m_kxc/kxc_ADA
1120 : !! NAME
1121 : !! kxc_ADA
1122 : !!
1123 : !! FUNCTION
1124 : !! Calculate exchange-correlation kernel in reciprocal space
1125 : !!
1126 : !! INPUTS
1127 : !! Dtset <type(dataset_type)>=all input variables in this dataset
1128 : !! Cryst<crystal_t>=Info on the unit cell.
1129 : !! ixc = choice for the exchange-correlation potential.
1130 : !! ngfft(18)=contain all needed information about 3D FFT,
1131 : !! see ~abinit/doc/variables/vargs.htm#ngfft
1132 : !! nfft = total number of points on the FFT grid.
1133 : !! rhor(nfft,nspden) = the charge density on the FFT grid.
1134 : !! (total in first half and spin-up in second half if nsppol=2)
1135 : !! npw: the size of kernel matrix
1136 : !! dim_kxcg=dimension of the kernel.
1137 : !! comm=MPI communicator.
1138 : !! [dbg_mode]=Set it to .TRUE. to switch on the debug mode.
1139 : !!
1140 : !! OUTPUT
1141 : !! kxcg(nfft,dim_kxcg) = the exchange-correlation potential on the FFT grid.
1142 : !! warning: the kernel is not divided by unit cell volume
1143 : !!
1144 : !! NOTES
1145 : !! No xc quadrature
1146 : !! No nl core correction
1147 : !!
1148 : !! SOURCE
1149 :
1150 0 : subroutine kxc_ADA(Dtset,Cryst,ixc,ngfft,nfft,nspden,rhor,&
1151 0 : npw,nqibz,qibz,fxc_ADA,gvec,comm,kappa_init,dbg_mode)
1152 :
1153 : !Arguments ------------------------------------
1154 : !scalars
1155 : integer,intent(in) :: ixc,nfft,nspden,npw,comm
1156 : real(dp),intent(in),optional :: kappa_init
1157 : logical,optional,intent(in) :: dbg_mode
1158 : type(crystal_t),intent(in) :: Cryst
1159 : type(dataset_type),intent(in) :: Dtset
1160 : !arrays
1161 : integer,intent(in) :: gvec(3,npw),ngfft(18)
1162 : integer,intent(in) :: nqibz
1163 : real(dp),intent(in) :: rhor(nfft,nspden)
1164 : real(dp),intent(in) :: qibz(3,nqibz)
1165 : complex(gwp),intent(out) :: fxc_ADA(npw,npw,nqibz)
1166 :
1167 : !Local variables ------------------------------
1168 : !scalars
1169 : integer :: i1,i2,i3,ig,igp,ir,irp,n3xccc,ngfft1,ngfft2,izero !,isp
1170 : integer :: ngfft3,nkxc,option,ikxc,ierr,nproc
1171 : integer :: nk3xc,igrid,iqbz,my_rank,master,gmgp_idx
1172 : logical :: non_magnetic_xc
1173 : real(dp) :: el_temp,bigexc,bigsxc,gsqcut,ucvol !,rs,Kx,Kc
1174 : real(dp) :: vxcavg,kappa,abs_qpg_sq,abs_qpgp_sq
1175 : real(dp) :: difx,dify,difz,inv_kappa_sq
1176 : character(len=500) :: msg,fname
1177 0 : type(MPI_type) :: MPI_enreg_seq
1178 : type(xcdata_type) :: xcdata
1179 : !arrays
1180 : real(dp) :: qpg(3),qpgp(3),qphon(3),q_point(3),dum(0)
1181 : real(dp),parameter :: dummyvgeo(3)=zero
1182 0 : real(dp),allocatable :: kxcr(:,:)
1183 0 : real(dp),allocatable :: rhog(:,:),vhartr(:),vxclda(:,:)
1184 0 : real(dp),allocatable :: xccc3d(:),my_rhor(:,:)
1185 0 : real(dp),allocatable :: my_kxcg(:,:)
1186 0 : real(dp),allocatable :: rhotilder(:,:)
1187 0 : complex(gwp),allocatable :: my_fxc_ADA_ggpq(:,:,:)
1188 0 : complex(gwp),allocatable :: FT_fxc_ADA_ggpq(:,:,:),dummy(:,:)
1189 0 : real(dp),allocatable :: rvec(:,:),my_fxc_ADA_rrp(:,:)
1190 : real(dp) :: rmrp(3),abs_rmrp
1191 0 : integer :: n1,n2,n3,ig_idx_fft(npw)
1192 : ! ************************************************************************
1193 :
1194 0 : ABI_CHECK_IEQ(Dtset%nsppol, 1,'nsppol/=1 not coded')
1195 0 : ABI_CHECK_IEQ(nspden, 1, 'nspden /=1 not coded')
1196 0 : ABI_CHECK_IEQ(nfft, PRODUCT(ngfft(1:3)), "mismatch in nfftot")
1197 :
1198 : ! Fake MPI_type for the sequential part.
1199 0 : call initmpi_seq(MPI_enreg_seq)
1200 :
1201 0 : my_rank = xmpi_comm_rank(comm)
1202 0 : nproc = xmpi_comm_size(comm)
1203 0 : master=0
1204 :
1205 0 : write(msg,'(a,i3)') ' kxc_ADA: calculating exchange-correlation kernel using ixc = ',ixc
1206 0 : call wrtout(std_out,msg)
1207 0 : call wrtout(std_out,' kxc_ADA: using smeared density')
1208 :
1209 0 : if (.not.present(kappa_init)) then
1210 0 : kappa = 2.1_dp
1211 : else
1212 0 : kappa = kappa_init
1213 : end if
1214 :
1215 0 : write(msg,'(a,F10.3)') ' kxc_ADA: inverse smearing length, kappa = ',kappa
1216 0 : call wrtout(std_out,msg)
1217 0 : inv_kappa_sq = one/(kappa*kappa)
1218 :
1219 0 : call xcdata_init(xcdata,dtset=dtset,intxc=0,ixc=ixc,nspden=nspden)
1220 :
1221 0 : if (ALL(xcdata%xclevel/=(/1,2/))) then
1222 0 : ABI_ERROR(sjoin("Unsupported xclevel: ", itoa(xcdata%xclevel)))
1223 : end if
1224 :
1225 0 : ngfft1=ngfft(1)
1226 0 : ngfft2=ngfft(2)
1227 0 : ngfft3=ngfft(3)
1228 :
1229 0 : non_magnetic_xc=(abs(dtset%usepawu)==4.or.dtset%usepawu==14)
1230 :
1231 0 : if (ixc>=1.and.ixc<11) then ! LDA case
1232 0 : nkxc= 2*min(Dtset%nspden,2)-1 ! 1 or 3
1233 0 : elseif (ixc==51) then ! TLDA case, same as LDA above
1234 0 : nkxc= 2*min(Dtset%nspden,2)-1
1235 : else ! GGA case
1236 0 : nkxc=12*min(Dtset%nspden,2)-5 ! 7 or 19
1237 0 : ABI_CHECK_IEQ(dtset%xclevel, 2,"Functional should be GGA")
1238 0 : ABI_ERROR("GGA functional not implemented for ADA vertex")
1239 : end if
1240 :
1241 0 : ABI_MALLOC(kxcr,(nfft,nkxc))
1242 :
1243 : !gsqcut and rhog are zeroed because they are not used by rhotoxc if 1<=ixc<=16 and option=0
1244 0 : gsqcut=zero
1245 :
1246 0 : ABI_MALLOC(rhog,(2,nfft))
1247 0 : ABI_MALLOC(vhartr,(nfft))
1248 0 : rhog(:,:)=zero
1249 : !MG FIXME this is the 3D core electron density for XC core correction (bohr^-3)
1250 : !should implement the non linear core correction
1251 0 : n3xccc=0
1252 0 : ABI_MALLOC(xccc3d,(n3xccc))
1253 0 : ABI_MALLOC(vxclda,(nfft,nspden))
1254 :
1255 0 : option=2 ! 2 for Hxc and kxcr (no paramagnetic part if nspden=1)
1256 : qphon(:)=zero
1257 :
1258 : !to be adjusted for the call to rhotoxc
1259 0 : nk3xc=1
1260 :
1261 : !Compute the kernel.
1262 0 : izero=0
1263 :
1264 : ! print density
1265 0 : if (present(dbg_mode)) then
1266 0 : if (dbg_mode.and.my_rank==master) then
1267 0 : fname = 'xc_ADA_den.xsf'
1268 0 : call cryst%write_xsf_data(fname, ngfft1, ngfft2, ngfft3, rhor(:,1))
1269 : end if
1270 : end if
1271 :
1272 : !Calculate the smeared density
1273 0 : ABI_MALLOC(my_rhor,(nfft,nspden))
1274 0 : ABI_MALLOC(rhotilder,(nfft,nspden))
1275 0 : ucvol = Cryst%ucvol
1276 0 : my_rhor = rhor
1277 : !do isp = 1,nsppol
1278 : !call calc_smeared_density(my_rhor(:,isp),1,rhotilder(:,isp),nfft,ngfft,npw,&
1279 : !& gvec,Cryst%gprimd,Cryst%ucvol,MPI_enreg_seq,paral_kgb0,kappa_in=kappa)
1280 : !my_rhor(:,isp) = rhotilder(:,isp)
1281 : !end do
1282 :
1283 : ! print smeared density
1284 0 : if (present(dbg_mode)) then
1285 0 : if (dbg_mode.and.my_rank==master) then
1286 0 : fname = 'xc_ADA_smeared_den.xsf'
1287 0 : call cryst%write_xsf_data(fname, ngfft1, ngfft2, ngfft3, my_rhor(:,1))
1288 : end if
1289 : end if
1290 :
1291 : ! Reinitialize the libxc module with the overridden values
1292 0 : if (dtset%ixc<0) then
1293 0 : call libxc_functionals_end()
1294 : end if
1295 0 : if (ixc<0) then
1296 0 : el_temp=merge(Dtset%tphysel,Dtset%tsmear,Dtset%tphysel>tol8.and.Dtset%occopt/=3.and.Dtset%occopt/=9)
1297 0 : call libxc_functionals_init(ixc,Dtset%nspden,el_temp=el_temp,xc_tb09_c=Dtset%xc_tb09_c)
1298 : end if
1299 :
1300 0 : call hartre(1,gsqcut,3,izero,MPI_enreg_seq,nfft,ngfft,1,zero,rhog,Cryst%rprimd,dummyvgeo,vhartr)
1301 : call rhotoxc(bigexc,bigsxc,kxcr,MPI_enreg_seq,nfft,ngfft,&
1302 : & dum,0,dum,0,nkxc,nk3xc,non_magnetic_xc,&
1303 : & n3xccc,option,my_rhor,Cryst%rprimd,&
1304 0 : & 1,vxclda,vxcavg,xccc3d,xcdata,vhartr=vhartr)
1305 :
1306 : !Check for extreme (NaN) values
1307 : !do ir=1,nfft
1308 : !if (isnan(kxcr(ir,1))) kxcr(ir,1) = HUGE(kxcr(ir,1))
1309 : !end do
1310 :
1311 : !DEBUG test with direct way of calculating Kxc
1312 : !do i1=1,nfft
1313 : !rs = (three/(four_pi*my_rhor(i1,1)))**third
1314 : !Kx = 16._dp/27._dp*0.3141592653589793e1_dp*(rs**2)*(-0.4581652_dp)
1315 : !
1316 : !Kc = -0.4e1_dp / 0.9e1_dp * 0.3141592654e1_dp * rs ** 4 &
1317 : !* (0.207271333333333333333333333333e-1_dp * &
1318 : !(-0.177442658629204480000000e3_dp * rs - 0.17565190511219200000000e2_dp &
1319 : !* sqrt(rs) - 0.1332650665120000e2_dp * rs ** 2 &
1320 : !- 0.51031691247948928000000e2_dp * rs ** (0.3e1_dp / 0.2e1_dp)) &
1321 : !* rs ** (-0.3e1_dp / 0.2e1_dp) / (rs + 0.37274400e1_dp * sqrt(rs) &
1322 : !+ 0.129352000e2_dp) ** 2 / (-sqrt(rs) - 0.1049800_dp) &
1323 : !+ 0.518178333333333333333333333333e-2_dp * rs ** (-0.3e1_dp / 0.2e1_dp) &
1324 : !* (0.617071835390850041282140897280e3_dp * sqrt(rs) &
1325 : !+ 0.659369347307557491857191871552e5_dp * rs ** 2 + &
1326 : !0.700403648491298930017835369562e5_dp * rs ** (0.3e1_dp / 0.2e1_dp) &
1327 : !+ 0.398437532951539263722720308167e5_dp * rs ** (0.5e1_dp / 0.2e1_dp) &
1328 : !+ 0.368852071032531998953472000000e4_dp * rs ** (0.7e1_dp / 0.2e1_dp) &
1329 : !+ 0.5330602660480000e2_dp * rs ** (0.9e1_dp / 0.2e1_dp) &
1330 : !+ 0.143783940386264738593799346176e5_dp * rs ** 3 &
1331 : !+ 0.124672564145568409213848436081e5_dp * rs &
1332 : !+ 0.557398029956167136000000e3_dp * rs ** 4) &
1333 : !/ (rs + 0.37274400e1_dp * sqrt(rs) + 0.129352000e2_dp) ** 4 &
1334 : !/ (-sqrt(rs) - 0.1049800_dp) ** 2)
1335 : !kxcr(i1,1) = Kx + Kc
1336 : !end do
1337 : !END DEBUG
1338 :
1339 : ! print Kxc
1340 0 : if (present(dbg_mode)) then
1341 0 : if (dbg_mode.and.my_rank==master) then
1342 0 : fname = 'xc_ADA_Kxc.xsf'
1343 0 : call cryst%write_xsf_data(fname, ngfft1, ngfft2, ngfft3, kxcr(:,1))
1344 : end if
1345 : end if
1346 :
1347 0 : ABI_FREE(xccc3d)
1348 0 : ABI_FREE(vxclda)
1349 0 : ABI_FREE(vhartr)
1350 :
1351 0 : ABI_MALLOC(my_kxcg,(2,nfft))
1352 :
1353 0 : do ikxc=1,nkxc
1354 0 : call fourdp(1,my_kxcg,kxcr(:,ikxc),-1,MPI_enreg_seq,nfft,1,ngfft,0)
1355 : !kxcg(:,ikxc)=CMPLX(my_kxcg(1,:),my_kxcg(2,:))
1356 : end do
1357 : !TODO Check symmetry of kxcg
1358 :
1359 : !set up ADA vertex
1360 0 : ABI_MALLOC(my_fxc_ADA_ggpq,(npw,npw,nqibz))
1361 0 : my_fxc_ADA_ggpq = czero
1362 :
1363 : !Calculate f_xc(R,R')=(kappa^2/2)K_xc[\tilde{n}](G-G')
1364 : !x(1/(kappa^2+|q+G|^2) + 1/(kappa^2+|q+G'|^2
1365 : !First get G vectors and indices
1366 :
1367 0 : ierr=0
1368 0 : do iqbz=1,nqibz
1369 0 : q_point(:) = qibz(:,iqbz)
1370 0 : do ig=1,npw
1371 0 : do igp=1,npw
1372 : ! Calculate |q+G| and |q+G'|
1373 0 : qpg(:) = two_pi*MATMUL(Cryst%gprimd,q_point(:)+gvec(:,ig))
1374 0 : qpgp(:) = two_pi*MATMUL(Cryst%gprimd,q_point(:)+gvec(:,igp))
1375 0 : abs_qpg_sq = 1.0_dp/(1.0_dp+dot_product(qpg,qpg)*inv_kappa_sq)
1376 0 : abs_qpgp_sq = 1.0_dp/(1.0_dp+dot_product(qpgp,qpgp)*inv_kappa_sq)
1377 :
1378 0 : gmgp_idx = g2ifft(gvec(:,ig)-gvec(:,igp),ngfft)
1379 0 : if (gmgp_idx>0) then
1380 0 : my_fxc_ADA_ggpq(ig,igp,iqbz) = half*CMPLX(my_kxcg(1,gmgp_idx), my_kxcg(2,gmgp_idx))*(abs_qpg_sq+abs_qpgp_sq)
1381 : else
1382 0 : ierr=ierr+1
1383 0 : my_fxc_ADA_ggpq(ig,igp,iqbz) = czero
1384 : end if
1385 : end do
1386 : end do
1387 0 : if (ierr/=0) then
1388 : write(msg,'(a,i4,3a)')&
1389 0 : ' Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1390 0 : ' Enlarge the FFT mesh to get rid of this problem. '
1391 0 : ABI_WARNING(msg)
1392 : end if
1393 : end do
1394 :
1395 0 : fxc_ADA = my_fxc_ADA_ggpq
1396 :
1397 : !do iqbz=1,nqibz
1398 : !call hermitianize(my_fxc_ADA_ggpq(:,:,iqbz),"All")
1399 : !end do
1400 :
1401 : !DEBUG check symmetry
1402 : if (.FALSE.) then
1403 : ! do iqbz=1,nkptgw
1404 : ! do ig=1,npw
1405 : ! do igp=ig,npw
1406 : ! if (ABS(REAL(fxc_ADA(ig,igp,iqbz))-REAL(fxc_ADA(igp,ig,iqbz)))>tol15.OR.&
1407 : ! ABS(AIMAG(fxc_ADA(ig,igp,iqbz))-AIMAG(-fxc_ADA(igp,ig,iqbz)))>tol15) then
1408 : ! write(std_out,*) 'Elements:'
1409 : ! write(std_out,*) 'fxc_ADA(ig,igp,iqbz):',ig,igp,iqbz,fxc_ADA(ig,igp,iqbz)
1410 : ! write(std_out,*) 'fxc_ADA(igp,ig,iqbz):',igp,ig,iqbz,fxc_ADA(igp,ig,iqbz)
1411 : ! ABI_ERROR('fxc_ADA not symmetric')
1412 : ! end if
1413 : ! end do
1414 : ! end do
1415 : ! end do
1416 :
1417 : ! write(std_out,*)"kxcr(r=0)",kxcr(1,1)
1418 : ! write(std_out,*)"my_kxg(G=0)",my_kxcg(:,1)
1419 : ! write(std_out,*)"SUM kxcr/nfft ",SUM(kxcr(:,1))/nfft
1420 : ! write(std_out,*)"SUM my_kxg ",SUM(kxcg(:,1))
1421 :
1422 : ! DEBUG Check FT to real space
1423 : ! The real-space expression is:
1424 : ! f_xc(R,R')=(1/2)(kappa^2/(4*Pi))
1425 : ! \{K_xc[\tilde{n(R)}]+K_xc[\tilde{n(R')}]\}
1426 : ! x exp(-kappa|R-R'|)/|R-R'|
1427 : ABI_MALLOC(my_fxc_ADA_rrp,(nfft,nfft))
1428 : ABI_MALLOC(FT_fxc_ADA_ggpq,(npw,npw,nqibz))
1429 : ABI_MALLOC(rvec,(3,nfft))
1430 : ABI_MALLOC(dummy,(nfft,nfft))
1431 : my_fxc_ADA_rrp=zero; FT_fxc_ADA_ggpq=czero; dummy=czero; rvec=zero
1432 :
1433 : ! First find coordinates of real-space fft points
1434 : igrid = 0
1435 : ngfft1 = ngfft(1)
1436 : ngfft2 = ngfft(2)
1437 : ngfft3 = ngfft(3)
1438 : do i3=0,ngfft3-1
1439 : difz=dble(i3)/dble(ngfft3)
1440 : do i2=0,ngfft2-1
1441 : dify=dble(i2)/dble(ngfft2)
1442 : do i1=0,ngfft1-1
1443 : difx=dble(i1)/dble(ngfft1)
1444 : igrid = igrid + 1
1445 : rvec(1,igrid)=difx*Cryst%rprimd(1,1)+dify*Cryst%rprimd(1,2)+difz*Cryst%rprimd(1,3)
1446 : rvec(2,igrid)=difx*Cryst%rprimd(2,1)+dify*Cryst%rprimd(2,2)+difz*Cryst%rprimd(2,3)
1447 : rvec(3,igrid)=difx*Cryst%rprimd(3,1)+dify*Cryst%rprimd(3,2)+difz*Cryst%rprimd(3,3)
1448 : end do
1449 : end do
1450 : end do
1451 : if (igrid/=nfft) then
1452 : ABI_ERROR('kxc_ADA: igrid not equal to nfft')
1453 : end if
1454 :
1455 : ! Construct kernel in real space
1456 : do ir=1,nfft
1457 : do irp=ir,nfft
1458 : rmrp(:) = rvec(:,ir)-rvec(:,irp)
1459 : abs_rmrp = sqrt(dot_product(rmrp,rmrp))
1460 : my_fxc_ADA_rrp(ir,irp) = eighth*kappa*kappa*piinv* &
1461 : (kxcr(ir,1)+kxcr(irp,1))* &
1462 : EXP(-kappa*abs_rmrp)/(abs_rmrp+1.e-3_dp)
1463 : ! (a small convergence factor is introduced
1464 : ! to avoid a singularity)
1465 : my_fxc_ADA_rrp(irp,ir) = my_fxc_ADA_rrp(ir,irp)
1466 : end do
1467 : end do
1468 :
1469 : ! Find FFT index for all G
1470 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1471 : ! Use the following indexing (N means ngfft of the adequate direction)
1472 : ! 0 1 2 3 ... N/2 -(N-1)/2 ... -1 <= kg
1473 : ! 1 2 3 4 ....N/2+1 N/2+2 ... N <= index
1474 : do ig=1,npw
1475 : i1=modulo(gvec(1,ig),n1)
1476 : i2=modulo(gvec(2,ig),n2)
1477 : i3=modulo(gvec(3,ig),n3)
1478 : ig_idx_fft(ig)=i1+1+n1*(i2+n2*i3)
1479 : end do
1480 : ! FT kernel to reciprocal space for each q
1481 : do iqbz=1,nqibz
1482 : dummy = CMPLX(my_fxc_ADA_rrp,0.0_dp)
1483 : ! Multiply with q-point phase factors exp(-iq.r)*f_xc(r,r')*exp(iq.r')
1484 : do ir=1,nfft
1485 : do irp=1,nfft
1486 : ! Calculate q (variables defined for other purposes
1487 : ! are being reused as dummy variables)
1488 : q_point(:) = qibz(:,iqbz)
1489 : qpg(:) = two_pi*MATMUL(Cryst%gprimd,q_point(:))
1490 : abs_qpg_sq = dot_product(qpg(:),rvec(:,ir))
1491 : abs_qpgp_sq = dot_product(qpg(:),rvec(:,irp))
1492 : dummy(ir,irp) = EXP(-j_dpc*abs_qpg_sq) * dummy(ir,irp)* EXP(j_dpc*abs_qpgp_sq)
1493 : end do
1494 : end do
1495 :
1496 : call fourdp_6d(2,dummy,-1,MPI_enreg_seq,nfft,ngfft, 0)
1497 :
1498 : do ig=1,npw
1499 : do igp=1,npw
1500 : FT_fxc_ADA_ggpq(ig,igp,iqbz) = dummy(ig_idx_fft(ig),ig_idx_fft(igp))
1501 : end do
1502 : end do
1503 :
1504 : ! Output
1505 : msg=''
1506 : if (iqbz<10) write(msg,'(a,i1,a)') './debug_fxc_ADA_q',iqbz,'.dat'
1507 : if ((iqbz>9).and.(iqbz<100)) write(msg,'(a,i2,a)') './debug_fxc_ADA_q',iqbz,'.dat'
1508 : if ((iqbz>99).and.(iqbz<1000)) write(msg,'(a,i3,a)') './debug_fxc_ADA_q',iqbz,'.dat'
1509 :
1510 : !open(777,file=TRIM(msg),STATUS='REPLACE')
1511 : !do igp=1,npw
1512 : ! do ig=1,npw
1513 : ! write(777,*) ig,igp,REAL(my_fxc_ADA_ggpq(ig,igp,iqbz)),AIMAG(my_fxc_ADA_ggpq(ig,igp,iqbz)), &
1514 : ! REAL(FT_fxc_ADA_ggpq(ig,igp,iqbz)),AIMAG(FT_fxc_ADA_ggpq(ig,igp,iqbz)), &
1515 : ! ABS(ABS(my_fxc_ADA_ggpq(ig,igp,iqbz))-ABS(FT_fxc_ADA_ggpq(ig,igp,iqbz)))
1516 : ! end do
1517 : ! write(777,*) ''
1518 : !end do
1519 : !close(777)
1520 :
1521 : end do ! iqbz
1522 :
1523 : ABI_ERROR('Stopping in kxc_ADA for debugging')
1524 :
1525 : ABI_FREE(rvec)
1526 : ABI_FREE(my_fxc_ADA_rrp)
1527 : ABI_FREE(FT_fxc_ADA_ggpq)
1528 :
1529 : if (xcdata%xclevel==2) then
1530 : ABI_ERROR(" GGA not implemented for kxc_ADA")
1531 : end if !xclevel==2
1532 :
1533 : end if ! Debugging section
1534 :
1535 : ! Revert libxc module to the original settings
1536 0 : if (ixc<0) then
1537 0 : call libxc_functionals_end()
1538 : end if
1539 0 : if (dtset%ixc<0) then
1540 0 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
1541 0 : call libxc_functionals_init(dtset%ixc,dtset%nspden,el_temp=el_temp,xc_tb09_c=dtset%xc_tb09_c)
1542 : end if
1543 :
1544 0 : ABI_FREE(my_kxcg)
1545 0 : ABI_FREE(my_rhor)
1546 0 : ABI_FREE(rhotilder)
1547 0 : ABI_FREE(rhog)
1548 0 : ABI_FREE(kxcr)
1549 :
1550 0 : call destroy_mpi_enreg(MPI_enreg_seq)
1551 :
1552 0 : end subroutine kxc_ADA
1553 : !!***
1554 :
1555 : !----------------------------------------------------------------------
1556 :
1557 : end MODULE m_kxc
|
