Branch data 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
|
