Line data Source code
1 : #define POTENTIAL
2 : MODULE param
3 : USE m_constants
4 : IMPLICIT NONE
5 : !
6 :
7 : REAL,PARAMETER :: au2A =0.529177249
8 : REAL,PARAMETER :: Ha2eV=hartree_to_ev_const
9 : REAL,PARAMETER :: pi=3.14159265358979323846
10 : !
11 : INTEGER, PARAMETER :: jnlout=6
12 : INTEGER, PARAMETER :: ftable=20,fdebug=21,fxcrysd=22,fwarn=23
13 : !
14 : ! parameters used to calculate the vdW-DF functional
15 : !
16 : REAL,PARAMETER :: Zab_v1=-0.8491 ! vdW-DF1
17 : REAL,PARAMETER :: Zab_v2=-1.887 ! vdW-DF2
18 : !
19 : ! parameters used to calculate the LDA correlation energy
20 : !
21 : REAL,PARAMETER :: A = 0.031091 ! some fitted parameters for the
22 : REAL,PARAMETER :: alpha_1 = 0.2137 ! evaluation of LDA correlation
23 : REAL,PARAMETER :: beta_1 = 7.5957 ! energy according to:
24 : REAL,PARAMETER :: beta_2 = 3.5876 !
25 : REAL,PARAMETER :: beta_3 = 1.6382 ! J. P. Perdew and Yue Wang,
26 : REAL,PARAMETER :: beta_4 = 0.49294 ! Phys. Rev. B 45, 13244 (1992)
27 : REAL,PARAMETER :: beta_H = 0.066725 ! for H function from PRL 77, 3865 (1996)
28 : REAL,PARAMETER :: gamma_H = 0.031091 ! for H function
29 : !
30 : ! parameters used to calculate the PBE exchange energy
31 : !
32 : REAL,PARAMETER :: kappa_revPBE = 1.245 ! revPBE
33 : REAL,PARAMETER :: kappa_PBE = 0.804 ! PBE from PRL 77, 3865 (1996)
34 : REAL,PARAMETER :: mu = 0.21951 ! PBE from PRL 77, 3865 (1996)
35 : !
36 : END MODULE param
37 : !
38 : ! module containing the main non-local variables and subroutines
39 : !
40 : MODULE nonlocal_data
41 :
42 : IMPLICIT NONE
43 : !
44 : REAL :: Zab ! This is the Zab really used. Zab_v1
45 : ! is the default. Can be switched to
46 : ! Zab_v2 by invoking the program with
47 : ! command line option 'vdw2'
48 : !
49 : INTEGER :: nx, ny, nz ! number of grid points in x, y, z direction
50 : INTEGER :: n_grid ! total number of grid points
51 : INTEGER :: n_k ! number of k-points for which the kernel was
52 : ! tabulated
53 : !
54 : REAL :: r_max ! maximum r for which the kernel has been generated
55 : !
56 : REAL,ALLOCATABLE :: q_alpha(:) ! q mesh for the interpolation
57 : REAL,ALLOCATABLE :: phi(:,:,:) !
58 : REAL,ALLOCATABLE :: d2phi_dk2(:,:,:) !
59 : !
60 : INTEGER :: n_gvectors ! number of G vectors
61 : REAL :: G_cut ! radius of cutoff sphere
62 : REAL, ALLOCATABLE :: G(:,:) ! the G vectors
63 : INTEGER,ALLOCATABLE :: G_ind(:) ! the index of the G vectors
64 : !
65 : REAL :: a1(3),a2(3),a3(3) ! lattice vectors
66 : REAL :: b1(3),b2(3),b3(3) ! reciprocal lattice vectors
67 : !
68 : INTEGER :: n_alpha ! number of q points
69 : REAL :: q_cut ! maximum q value
70 : INTEGER :: m_c ! maximum m for the saturation of q_0
71 : !
72 : REAL :: omega ! volume of the unit cell
73 : REAL :: tpibya ! (2*pi/omega)
74 : !
75 : REAL :: lambda ! parameter for the logarithmic q mesh
76 : REAL :: dk ! spacing of the uniform radial k grid
77 : !
78 : INTEGER :: time1(8),time2(8) ! arrays for measuring the execution time
79 : !
80 : END MODULE nonlocal_data
81 : !
82 : MODULE driver_fft
83 : !
84 : CONTAINS
85 : !
86 : !==========================================================================================
87 : !
88 : ! this subroutine is an interface for the fft. It will fourier transform in place the
89 : ! complex array fftin. idir = -1 means forward and idir = +1 means backward fourier
90 : ! transform. At the moment it is designed to use fftw3.
91 :
92 0 : SUBROUTINE inplfft( fftin, idir )
93 : USE m_juDFT
94 : USE param, ONLY: jnlout
95 : USE nonlocal_data,ONLY: nx,ny,nz,n_grid
96 :
97 : implicit none
98 :
99 :
100 : complex, intent(inout) :: fftin(n_grid)
101 :
102 : integer :: idir
103 :
104 : integer*8 :: plan
105 : #ifdef CPP_FFTW
106 : include 'fftw3.f' ! some definitions for fftw3
107 :
108 0 : if ( idir == -1 ) then
109 :
110 0 : CALL dfftw_plan_dft_3d(plan,nz,ny,nx,fftin,fftin,FFTW_FORWARD,FFTW_ESTIMATE)
111 0 : CALL dfftw_execute_dft(plan,fftin,fftin)
112 0 : CALL dfftw_destroy_plan(plan)
113 :
114 0 : fftin = fftin / real(n_grid) ! rescale as fftw3 puts a factor n_grid on forward FFT
115 :
116 0 : elseif( idir == 1 ) then
117 :
118 0 : CALL dfftw_plan_dft_3d(plan,nz,ny,nx,fftin,fftin,FFTW_BACKWARD,FFTW_ESTIMATE)
119 0 : CALL dfftw_execute_dft(plan,fftin,fftin)
120 0 : CALL dfftw_destroy_plan(plan)
121 :
122 : else
123 :
124 : write(jnlout,*) 'ERROR during FFT: neither FORWARD &
125 0 : nor BACKWARD FFT was chosen.'
126 0 : STOP 'Error in FFT'
127 :
128 : endif
129 : #else
130 : CALL judft_error("VdW functionals only available if compiled with CPP_FFTW")
131 : #endif
132 0 : END SUBROUTINE inplfft
133 : !
134 : END MODULE driver_fft
135 : !
136 : MODULE functionals
137 : !
138 : CONTAINS
139 : !
140 : !================================================================================
141 : !
142 0 : SUBROUTINE calc_PBE_correlation(n,grad_n,Ec_PBE,Ec_LDA,e_cLDA,e_cSL)
143 : !
144 : USE param, ONLY: Ha2eV,pi, &
145 : jnlout, &
146 : beta_H,gamma_H, & ! Parameters for LDA_c
147 : A,alpha_1,beta_1,beta_2,beta_2,beta_3,beta_4 ! Parameters for LDA_c
148 : USE nonlocal_data,ONLY: n_grid,omega
149 : !
150 : IMPLICIT NONE
151 : !
152 : REAL,INTENT(in) :: n(n_grid) ! charge density
153 : REAL,INTENT(in) :: grad_n(n_grid,3) ! gradient of charge density
154 : REAL,INTENT(out) :: Ec_PBE ! PBE-correlation energy
155 : REAL,INTENT(out) :: Ec_LDA ! LDA-correlation energy
156 : REAL,INTENT(out) :: e_cLDA(n_grid) ! LDA correlation energy density
157 : REAL,INTENT(out) :: e_cSL(n_grid) ! SL correlation energy density in Ha
158 : !
159 : REAL :: r_s ! intermediate values needed for the
160 : REAL :: k_F ! formulas given below.
161 : REAL :: zeta !
162 : REAL :: phi_zeta !
163 : REAL :: grad_n_squ !
164 : REAL :: t ! -- "" --
165 : REAL :: k_s !
166 : REAL :: AA !
167 : REAL :: H !
168 : !
169 : REAL :: sqrt_r_s
170 : REAL :: eLDA_c
171 : REAL :: LDA_dummy_1
172 : REAL :: LDA_dummy_2
173 : !
174 : integer :: i1
175 : !
176 0 : e_cLDA(:)=0.0
177 0 : e_cSL(:) =0.0
178 : !
179 0 : Ec_PBE=0.0
180 0 : Ec_LDA=0.0
181 : !
182 : ! Formula for PBE correlation ( Perdew, Burke and Ernzerhof PRL 77, 18 (1996) eqn. [3]):
183 : ! Ec_PBE = int d^3 r n(r) (ec_LSDA (r_s, zeta) + H (r_s, zeta, t))
184 : !
185 : ! r_s = 3 / (4* pi*n)**(1/3) Seitz radius.
186 : ! zeta is relative spin polarization will be set to 0 as vdW-DF only works without spin
187 : ! t = | grad_n | / (2*phi(zeta)*k_s*n) is a dimensionless gradient
188 : ! phi(zeta) = [ (1 + zeta)**(2/3) + (1 - zeta)**(2/3) ] / 2
189 : ! k_s = sqrt( 4*k_F / (pi*a_0) )
190 : ! a_0 = ( hbar / (m*e^2) ) = 1 in a. u.
191 : !
192 : ! H = (e**2/a_0) g * phi**3
193 : ! * ln(1 + beta_H/gamma_H * t**2 * [1 + AA*t**2 / (1 + AA*t**2 + AA**2 * t**4 )] )
194 : ! AA = beta_H/gamma_H * [ exp( -ec_LSDA / ( gamma_H * phi**3 * e**2 / a_0)) - 1 ]**(-1)
195 : !
196 : !
197 0 : zeta = 0.0 ! non-spin polarized case
198 : !
199 : phi_zeta = ( (1.0 + zeta)**(2.0/3.0) + &
200 0 : (1.0 - zeta)**(2.0/3.0))/2.0
201 : !
202 0 : do i1=1,n_grid
203 : !
204 0 : if(n(i1) < 1e-12) cycle
205 :
206 0 : r_s = (3.0/(4.0*pi*n(i1)))**(1.0/3.0)
207 :
208 : ! we also need LDA correlation per particle so I repeat here the formulas from calc_q0
209 :
210 0 : sqrt_r_s = sqrt(r_s)
211 :
212 0 : LDA_dummy_1 = (8.0*pi/3.0)*A*(1.0+ alpha_1*r_s)
213 : LDA_dummy_2 = 2.0*A*(beta_1*sqrt_r_s + beta_2*r_s + &
214 0 : beta_3*sqrt_r_s*r_s + beta_4*r_s*r_s)
215 :
216 0 : eLDA_c = (-3.0/(4.0*pi))*LDA_dummy_1*log(1.0+1.0/LDA_dummy_2)
217 :
218 : ! now start calculation of PBE correlation. first check wether density is very small (approx
219 : ! zero) or negative which is also not correct
220 : !
221 0 : k_F = (3.0*pi*pi*n(i1))**(1.0/3.0)
222 :
223 0 : k_s = sqrt( 4.0 * k_F / pi)
224 :
225 0 : grad_n_squ = grad_n(i1,1)**2 + grad_n(i1,2)**2 + grad_n(i1,3)**2
226 :
227 0 : t = sqrt(grad_n_squ)/(2.0*phi_zeta*k_s*n(i1))
228 :
229 : AA = (beta_H/gamma_H)* &
230 0 : (1.0/(exp(-1.0*eLDA_c/(gamma_H*phi_zeta**3))- 1.0 ))
231 :
232 : H = gamma_H*phi_zeta**3 * &
233 : log(1.0+ &
234 : (beta_H/gamma_H)*(t**2)* &
235 0 : (1.0+AA*t**2)/(1.0+AA*t**2+AA**2 * t**4))
236 : !
237 : ! LDA and SL correlation energy density
238 : !
239 0 : e_cLDA(i1) = eLDA_c*n(i1)
240 0 : e_cSL(i1) = H*n(i1)
241 : !
242 : #if 0
243 : Ec_PBE = Ec_PBE + n(i1)*eLDA_c + n(i1)*H
244 : Ec_LDA = Ec_LDA + n(i1)*eLDA_c
245 : #else
246 0 : Ec_PBE = Ec_PBE + e_cLDA(i1) + e_cSL(i1)
247 0 : Ec_LDA = Ec_LDA + e_cLDA(i1)
248 : #endif
249 : enddo
250 : !
251 0 : Ec_PBE = Ec_PBE*omega/real(n_grid)
252 0 : Ec_LDA = Ec_LDA*omega/real(n_grid)
253 : !
254 0 : END SUBROUTINE calc_PBE_correlation
255 : !
256 : !================================================================================
257 : !
258 0 : SUBROUTINE calc_GGA_exchange(n,grad_n,Ex_PBE,Ex_revPBE,Ex_PW86,Ex_LDA)
259 : !
260 : USE param, ONLY: pi, &
261 : mu,kappa_PBE,kappa_revPBE ! Param. for PBE_ex
262 : USE nonlocal_data,ONLY: n_grid,omega
263 : !
264 : IMPLICIT NONE
265 : !
266 : REAL,INTENT(IN) :: n(n_grid) ! charge density
267 : REAL,INTENT(IN) :: grad_n(n_grid,3) ! gradient of charge density
268 : REAL,INTENT(OUT) :: Ex_PBE ! PBE-exchange energy in Ha
269 : REAL,INTENT(OUT) :: Ex_revPBE ! revPBE-exchange energy in Ha
270 : REAL,INTENT(OUT) :: Ex_PW86 ! refit PW86-exchange energy
271 : ! in Ha
272 : REAL,INTENT(OUT) :: Ex_LDA ! LDA part of the PBE exchange
273 : ! energy Ex_PBE
274 : !
275 : INTEGER :: i1
276 : REAL :: k_F ! formulas given below.
277 : REAL :: ss ! reduced density gradient
278 : REAL :: sqrt_grad_n
279 : !
280 : ! ss is actually the reduced density gradient
281 : ! ss = |\grad n|/2(3\pi^2)^1/3n^4/3
282 : !
283 0 : Ex_PBE =0.0
284 0 : Ex_revPBE=0.0
285 0 : Ex_PW86 =0.0
286 0 : Ex_LDA =0.0
287 : !
288 : ! here the charge density has already the correct indices
289 : !
290 0 : do i1=1,n_grid
291 : !
292 0 : if(n(i1) < 1e-12) cycle
293 : !
294 0 : k_F = (3.0*pi*pi*n(i1))**(1.0/3.0)
295 : !
296 0 : sqrt_grad_n = sqrt(grad_n(i1,1)**2 + grad_n(i1,2)**2 + grad_n(i1,3)**2)
297 : !
298 0 : ss=sqrt_grad_n/(2.0*k_F*n(i1))
299 : !
300 : ! LDA part of the PBE exchange energy density
301 : !
302 : Ex_LDA = Ex_LDA + &
303 : (-3.0/4.0)*(3.0/pi)**(1.0/3.0)* &
304 0 : (n(i1)**(4.0/3.0))
305 : !
306 : ! PBE exchange energy
307 : !
308 : Ex_PBE = Ex_PBE + &
309 : (-3.0/4.0)*(3.0/pi)**(1.0/3.0)* &
310 : (n(i1)**(4.0/3.0))* &
311 : (1.0+kappa_PBE- &
312 0 : kappa_PBE/(1.0+mu*ss**2.0/kappa_PBE))
313 : !
314 : ! revPBE exchange energy
315 : !
316 : Ex_revPBE = Ex_revPBE + &
317 : (-3.0/4.0)*(3.0/pi)**(1.0/3.0)* &
318 : (n(i1)**(4.0/3.0))* &
319 : (1.0+kappa_revPBE- &
320 0 : kappa_revPBE/(1.0+mu*ss**2.0/kappa_revPBE))
321 : !
322 : ! use refit PW86
323 : !
324 : Ex_PW86 = Ex_PW86 + &
325 : (-3.0/4.0)*(3.0/pi)**(1.0/3.0)* &
326 : (n(i1)**(4.0/3.0))* &
327 : (1.0+15.0*0.1234*ss**2.0 + &
328 0 : 17.33*ss**4.0+0.163*ss**6.0)**(1.0/15.0)
329 : enddo
330 : !
331 0 : Ex_PBE = Ex_PBE *omega/real(n_grid)
332 0 : Ex_revPBE = Ex_revPBE*omega/real(n_grid)
333 0 : Ex_PW86 = Ex_PW86 *omega/real(n_grid)
334 0 : Ex_LDA = Ex_LDA *omega/real(n_grid)
335 : !
336 0 : END SUBROUTINE calc_GGA_exchange
337 : !
338 : !==========================================================================================
339 : !
340 0 : SUBROUTINE calc_ehartree(n)
341 : USE param, ONLY: Ha2ev,pi, &
342 : jnlout
343 : USE nonlocal_data,ONLY: n_grid,n_gvectors, &
344 : omega, &
345 : G,G_ind
346 : USE driver_fft
347 : !
348 : IMPLICIT NONE
349 : !
350 : REAL,INTENT(in) :: n(n_grid)
351 : !
352 : INTEGER :: i1,ind
353 : REAL :: E_Hartree
354 : REAL :: fac,size
355 0 : COMPLEX,ALLOCATABLE :: n_cmplx(:)
356 : !
357 : ! E_Hartree= 2*pi*omega \sum_{G/=0} n(G)^2/G^2
358 : !
359 0 : allocate(n_cmplx(n_grid))
360 : !
361 0 : n_cmplx = (0.0,0.0)
362 0 : n_cmplx(1:n_grid) = cmplx( n(1:n_grid), 0.0 )
363 : !
364 0 : call inplfft( n_cmplx, -1 )
365 : !
366 0 : E_Hartree=0.0
367 : !
368 0 : do i1 = 1, n_gvectors
369 : !
370 0 : ind = G_ind(i1)
371 0 : size=sum(G(i1,1:3)**2)
372 0 : if (size > 1e-12) then
373 0 : fac=1.0/size
374 : else
375 0 : print'(A,F10.3)','Total nr. of electrons: ',omega*real(n_cmplx(ind))
376 0 : fac=0.0
377 : endif
378 : !
379 : E_Hartree=E_Hartree+ &
380 0 : (real(n_cmplx(ind))**2+aimag(n_cmplx(ind))**2)*fac
381 : !
382 : enddo
383 : !
384 0 : deallocate(n_cmplx)
385 : !
386 0 : E_Hartree = E_Hartree*omega*2.0*pi
387 : !
388 0 : write(jnlout,'(/,A)') '-------------- Hartree energy --------------'
389 0 : write(jnlout,'(A)') 'E_Hartree (Ha/Ry/eV):'
390 0 : write(jnlout,'(3F24.15)') E_Hartree,E_Hartree*2.0,E_Hartree*Ha2eV
391 : !
392 0 : END SUBROUTINE calc_ehartree
393 : !
394 : END MODULE functionals
395 : !
396 : MODULE plot_functions
397 : !
398 : CONTAINS
399 : !
400 : !==========================================================================================
401 : !
402 0 : SUBROUTINE write_xcrysden_LDA(file_name,ene_dens)
403 : !
404 : USE param, ONLY: au2A,fxcrysd
405 : USE nonlocal_data,ONLY: nx,ny,nz,n_grid, &
406 : a1,a2,a3
407 : !
408 : IMPLICIT NONE
409 : !
410 : CHARACTER(LEN=*),INTENT(IN) :: file_name
411 : REAL, INTENT(IN) :: ene_dens(n_grid)
412 : INTEGER :: i1,ind,ix,iy,iz
413 0 : REAL,ALLOCATABLE :: ene_temp(:,:,:)
414 : !
415 0 : allocate(ene_temp(nx,ny,nz))
416 : !
417 0 : open(fxcrysd,FILE=trim(file_name))
418 : !
419 0 : write(fxcrysd,'(A)') 'CRYSTAL'
420 0 : write(fxcrysd,'(A)') 'PRIMVEC'
421 0 : write(fxcrysd,'(3F16.10)') a1(:)*au2A
422 0 : write(fxcrysd,'(3F16.10)') a2(:)*au2A
423 0 : write(fxcrysd,'(3F16.10)') a3(:)*au2A
424 : !
425 0 : write(fxcrysd,'(A)') 'PRIMCOORD'
426 0 : write(fxcrysd,'(A,2x,I1)') 'XXX',1
427 0 : write(fxcrysd,'(A)') 'Atomic_positions'
428 : !
429 0 : write(fxcrysd,'(A)') 'BEGIN_BLOCK_DATAGRID_3D'
430 0 : write(fxcrysd,'(A)') '3D_PWSCF'
431 0 : write(fxcrysd,'(A)') 'DATAGRID_3D_UNKNOWN'
432 0 : write(fxcrysd,'(3(3x,I5))') nx,ny,nz
433 0 : write(fxcrysd,'(3F16.10)') 0.0,0.0,0.0 ! origin of the system
434 0 : write(fxcrysd,'(3F16.10)') a1(:)*au2A
435 0 : write(fxcrysd,'(3F16.10)') a2(:)*au2A
436 0 : write(fxcrysd,'(3F16.10)') a3(:)*au2A
437 : !
438 0 : do ix=1,nx
439 0 : do iy=1,ny
440 0 : do iz=1,nz
441 0 : ind = ix + nx*(iy - 1) + nx*ny*(iz - 1)
442 0 : ene_temp(ix,iy,iz)=ene_dens(ind)
443 : ! in the case of the semi-local corr. ene. dens.
444 : ! ene_temp(ix,iy,iz) can be positive
445 : !if(ene_temp(ix,iy,iz) > 0.0) then
446 : ! print'(A)','WARNING: ene_temp(ix,iy,iz) > 0.0!'
447 : !endif
448 : enddo
449 : enddo
450 : enddo
451 : !
452 0 : do i1=1,nz
453 0 : write(fxcrysd,'(5(1x,E15.8))') ene_temp(1:nx,1:ny,i1)
454 : enddo
455 : !
456 : ! write the end of the XSF file
457 : !
458 0 : write(fxcrysd,'(A)') 'END_DATAGRID_3D'
459 0 : write(fxcrysd,'(A)') 'END_BLOCK_DATAGRID_3D'
460 : !
461 0 : close(fxcrysd)
462 : !
463 0 : deallocate(ene_temp)
464 : !
465 0 : END SUBROUTINE write_xcrysden_LDA
466 : !
467 : !==========================================================================================
468 : !
469 0 : SUBROUTINE write_xcrysden_NL(file_name,ene_dens)
470 : !
471 : USE param, ONLY: au2A,fxcrysd
472 : USE nonlocal_data,ONLY: nx,ny,nz, &
473 : n_grid,n_alpha, &
474 : a1,a2,a3
475 : !
476 : IMPLICIT NONE
477 : !
478 : CHARACTER(LEN=*),INTENT(IN) :: file_name
479 : REAL, INTENT(IN) :: ene_dens(n_grid,n_alpha)
480 : INTEGER :: i1,ind,ix,iy,iz
481 0 : REAL,ALLOCATABLE :: ene_temp(:,:,:)
482 : !
483 0 : allocate(ene_temp(nx,ny,nz))
484 : !
485 0 : open(fxcrysd,FILE=trim(file_name))
486 : !
487 0 : write(fxcrysd,'(A)') 'CRYSTAL'
488 0 : write(fxcrysd,'(A)') 'PRIMVEC'
489 0 : write(fxcrysd,'(3F16.10)') a1(:)*au2A
490 0 : write(fxcrysd,'(3F16.10)') a2(:)*au2A
491 0 : write(fxcrysd,'(3F16.10)') a3(:)*au2A
492 : !
493 0 : write(fxcrysd,'(A)') 'PRIMCOORD'
494 0 : write(fxcrysd,'(A,2x,I1)') 'XXX',1
495 0 : write(fxcrysd,'(A)') 'Atomic_positions'
496 : !
497 0 : write(fxcrysd,'(A)') 'BEGIN_BLOCK_DATAGRID_3D'
498 0 : write(fxcrysd,'(A)') '3D_PWSCF'
499 0 : write(fxcrysd,'(A)') 'DATAGRID_3D_UNKNOWN'
500 0 : write(fxcrysd,'(3(3x,I5))') nx,ny,nz
501 0 : write(fxcrysd,'(3F16.10)') 0.0,0.0,0.0 ! origin of the system
502 0 : write(fxcrysd,'(3F16.10)') a1(:)*au2A
503 0 : write(fxcrysd,'(3F16.10)') a2(:)*au2A
504 0 : write(fxcrysd,'(3F16.10)') a3(:)*au2A
505 : !
506 0 : do ix=1,nx
507 0 : do iy=1,ny
508 0 : do iz=1,nz
509 0 : ind = ix + nx*(iy - 1) + nx*ny*(iz - 1)
510 0 : ene_temp(ix,iy,iz)=sum(ene_dens(ind,1:n_alpha))
511 : enddo
512 : enddo
513 : enddo
514 : !
515 0 : do i1=1,nz
516 0 : write(fxcrysd,'(5(1x,E15.8))') ene_temp(1:nx,1:ny,i1)
517 : enddo
518 : !
519 : ! write the end of the XSF file
520 : !
521 0 : write(fxcrysd,'(A)') 'END_DATAGRID_3D'
522 0 : write(fxcrysd,'(A)') 'END_BLOCK_DATAGRID_3D'
523 : !
524 0 : close(fxcrysd)
525 : !
526 0 : deallocate(ene_temp)
527 : !
528 0 : END SUBROUTINE write_xcrysden_NL
529 : !
530 : END MODULE plot_functions
531 : !
532 : MODULE nonlocal_funct
533 : PRIVATE
534 : PUBLIC :: soler
535 : !
536 : CONTAINS
537 : !
538 : !========================================================================================
539 : !
540 0 : SUBROUTINE soler(n, Ecnl, v_nl)
541 : USE param, ONLY: Ha2eV,au2A,jnlout
542 : USE nonlocal_data,ONLY: n_grid,n_alpha, &
543 : omega, &
544 : q_alpha,phi, &
545 : G,G_ind, &
546 : time1,time2
547 : USE functionals
548 : USE driver_fft
549 : USE plot_functions
550 : !
551 : IMPLICIT NONE
552 : !
553 : REAL, INTENT(INOUT) :: n(n_grid) ! charge density
554 : REAL, INTENT(INOUT) :: v_nl(n_grid)
555 : REAL :: Ecnl ! the non-local correlation
556 :
557 : CHARACTER(LEN=132) :: option
558 : !
559 0 : REAL,ALLOCATABLE :: grad_n(:,:) ! gradient of charge density
560 0 : REAL,ALLOCATABLE :: e_cLDA(:) ! LDA correlation energy density
561 0 : REAL,ALLOCATABLE :: e_cSL(:) ! semi local correlation energy density
562 0 : REAL,ALLOCATABLE :: e_cNL(:,:) ! non local correlation energy density
563 0 : REAL,ALLOCATABLE :: q_0(:) ! array holding the q0 on the grid
564 :
565 : #ifdef POTENTIAL
566 :
567 0 : REAL,ALLOCATABLE :: dq0_dn(:)
568 0 : REAL,ALLOCATABLE :: dq0_dgrad_n(:)
569 :
570 : #endif
571 :
572 0 : COMPLEX,ALLOCATABLE :: theta_alpha(:,:)! array holding theta_i
573 0 : COMPLEX,ALLOCATABLE :: u_a(:,:) ! quantity needed for the non local
574 : ! correlation energy density
575 : REAL :: Ecnl_rsp ! Ecnl calculated in real space
576 : REAL :: Ec_LDA_rsp ! Ec_LDA calculated in real space
577 : REAL :: Ec_LDA ! LDA correlation
578 : REAL :: Ec_PBE ! PBE correlation energy. will be replaced by
579 : ! LDA correlation + non local correlation
580 : REAL :: Ex_PBE ! PBE exchange energy
581 : REAL :: Ex_revPBE ! revPBE echange energy
582 : REAL :: Ex_PW86 ! refit PW86 exchange energy
583 : REAL :: Ex_LDA ! LDA part of the PBE exchange
584 : ! energy Ex_PBE
585 : INTEGER :: i, j, k
586 : INTEGER :: ind ,mem_size
587 : REAL :: sum_test
588 : !
589 : ! the actual program starts here. allocation for n and grad_n have to be reconsidered when
590 : ! putting this as a subroutine into a different program.
591 : !
592 : ! read the pregenerated kernel from file 'vdW_kernel_table'. This will set variables n_alpha
593 : ! n_k and dk. The arrays phi(n_alpha,n_alpha,n_k) and d2phi_dk2(n_alpha,n_alpha,n_k) needed for the
594 : ! later interpolation of the kernel will be allocated and read.
595 : !
596 0 : call read_kernel()
597 : !
598 0 : allocate(theta_alpha(n_grid,n_alpha))
599 0 : allocate(u_a(n_grid,n_alpha))
600 0 : allocate(grad_n(n_grid,3),q_0(n_grid))
601 0 : allocate(e_cLDA(n_grid),e_cSL(n_grid),e_cNL(n_grid,n_alpha))
602 : !
603 : mem_size=2*n_grid*n_alpha + 2*n_grid*n_alpha + n_grid*3 + n_grid + &
604 0 : n_grid + n_grid + n_grid*n_alpha
605 : write(jnlout,'(A,F12.3,A,/)') &
606 0 : 'Memory required: ',real(mem_size)*8.0/(1024.0*1024.0),' MB'
607 : !
608 : ! setup the G vectors which we need for the calculation of the gradient, the non local correlation
609 : ! energy and the non local contribution to the potential.
610 : !
611 0 : call DATE_AND_TIME(values=time1)
612 : !
613 0 : call setup_g_vectors()
614 : !
615 0 : call DATE_AND_TIME(values=time2)
616 : !
617 0 : write(jnlout,'(A)') 'time to setup G vectors:'
618 0 : call timing( time2 - time1 )
619 : !
620 0 : write(jnlout,*)
621 : !
622 : ! calculate the gradient of the density numerically by FFT.
623 : !
624 0 : call DATE_AND_TIME(values=time1)
625 : !
626 0 : call calc_gradient( n, grad_n )
627 : !
628 0 : call DATE_AND_TIME(values=time2)
629 : !
630 0 : write(jnlout,'(A)') 'time to calculate gradient of n numerically:'
631 0 : call timing( time2 - time1 )
632 : !
633 : ! calculate q_0 for every grid point according to equations (11),(12) from Dion and (7)
634 : ! from Soler. The latter cares for the saturation.
635 : !
636 0 : call DATE_AND_TIME(values=time1)
637 : !
638 0 : call calc_q0(n, grad_n, q_0)
639 : !
640 0 : call DATE_AND_TIME(values=time2)
641 : !
642 0 : write(jnlout,'(A)') 'time to calculate q0:'
643 0 : call timing( time2 - time1 )
644 : !
645 : ! calculate theta_i defined as theta_i = n(r_i) * p_alpha(q_i) where p_alpha are the
646 : ! polynomials obtained by interpolating dirac delta. These polynomials will also be
647 : ! derived in this subroutine
648 : !
649 0 : call DATE_AND_TIME(values=time1)
650 : !
651 0 : call calc_theta_i(n, q_alpha, q_0, theta_alpha)
652 : !
653 0 : call DATE_AND_TIME(values=time2)
654 : !
655 0 : write(jnlout,'(A)') 'time to setup Thetas:'
656 0 : call timing( time2 - time1 )
657 : !
658 : ! to calculate the non local energy density we need theta_alpha_i in real space. As they are
659 : ! overwritten during the fourier transform we save them here.
660 : !
661 0 : e_cNL = theta_alpha
662 : !
663 : ! Fourier transform the theta_alpha_i in order to get theta_alpha_k.
664 : !
665 0 : call DATE_AND_TIME(values=time1)
666 : !
667 0 : do i=1,n_alpha
668 :
669 0 : call inplfft( theta_alpha(:,i), -1 )
670 :
671 : enddo
672 : !
673 0 : call DATE_AND_TIME(values=time2)
674 : !
675 0 : write(jnlout,'(A)') 'time spent to Fourier transform Thetas:'
676 0 : call timing( time2 - time1 )
677 : !
678 : ! now we can evaluate the integral in eqn. (8) of Soler.
679 : !
680 0 : call calc_ecnl(Ecnl, theta_alpha, u_a, e_cNL)
681 : !
682 : #ifdef POTENTIAL
683 :
684 0 : allocate(dq0_dn(n_grid))
685 0 : allocate(dq0_dgrad_n(n_grid))
686 :
687 0 : call calc_dq0(n, grad_n, dq0_dn, dq0_dgrad_n)
688 :
689 0 : call calc_potential(v_nl, u_a, q_0, dq0_dn, dq0_dgrad_n, n, grad_n)
690 :
691 0 : deallocate( dq0_dn, dq0_dgrad_n)
692 :
693 : #endif
694 :
695 : ! In vdW-DF GGA-correlation is been replaced by LDA-correlation + non local correlation. Thus
696 : ! we have to calculate Ec_LDA and Ec_PBE now to output th energy term which has to be added
697 : ! to the total energy. We already have the LDA correlation energy density from calc_q0 so we
698 : ! just have to integrate that array and afterwards call a subroutine calculating the PBE.
699 : ! omega/n_grid is the according volume element.
700 : !
701 0 : call calc_PBE_correlation(n,grad_n,Ec_PBE,Ec_LDA,e_cLDA,e_cSL)
702 : !
703 0 : write(jnlout,'(/,A)') 'PBE_c (Ha/Ry/eV):'
704 0 : write(jnlout,'(3F24.15)') Ec_PBE,Ec_PBE*2.0,Ec_PBE*Ha2eV
705 : !
706 0 : write(jnlout,'(/,A)') 'LDA_c (Ha/Ry/eV):'
707 0 : write(jnlout,'(3F24.15)') Ec_LDA,Ec_LDA*2.0,Ec_LDA*Ha2eV
708 : !
709 : ! for testing write the sum over LDA energy density.
710 : !
711 0 : Ec_LDA_rsp = sum(e_cLDA)*omega/real(n_grid)
712 : write(jnlout,'(/,A)') &
713 0 : '(check) LDA_c evaluated in real space (Ha/Ry):'
714 0 : write(jnlout,'(2F24.15)') Ec_LDA_rsp,Ec_LDA_rsp*2.0
715 : !
716 0 : write(jnlout,'(/,A)') 'E_c,nl (Ha/Ry/eV):'
717 0 : write(jnlout,'(3F24.15)') Ecnl,Ecnl*2.0,Ecnl*Ha2eV
718 : !
719 : ! for testing write the sum over non local energy density.
720 : ! it should be equal to Ec_nl
721 : !
722 0 : Ecnl_rsp=sum(e_cNL)*omega/real(n_grid)
723 : write(jnlout,'(/,A)') &
724 0 : '(check) E_c,nl evaluated in real space (Ha/Ry):'
725 0 : write(jnlout,'(2F24.15)') 0.5*Ecnl_rsp,Ecnl_rsp
726 : !
727 0 : call write_xcrysden_LDA('Ec_LDA.xsf',e_cLDA)
728 0 : call write_xcrysden_LDA('Ec_SL.xsf' ,e_cSL)
729 0 : call write_xcrysden_NL ('Ec_NL.xsf' ,e_cNL)
730 : !
731 0 : call calc_GGA_exchange(n,grad_n,Ex_PBE,Ex_revPBE,Ex_PW86,Ex_LDA)
732 : !
733 0 : write(jnlout,'(/,A)') 'Ex_LDA (Ha/Ry/eV):'
734 0 : write(jnlout,'(3F24.15)') Ex_LDA,Ex_LDA*2.0,Ex_LDA*Ha2eV
735 : !
736 0 : write(jnlout,'(/,A)') '-------------- PBE exchange --------------'
737 0 : write(jnlout,'(A)') 'Ex_PBE (Ha/Ry/eV):'
738 0 : write(jnlout,'(3F24.15)') Ex_PBE,Ex_PBE*2.0,Ex_PBE*Ha2eV
739 : !
740 0 : write(jnlout,'(A)') 'Ex_PBE+Ec_LDA+Ecnl (Ha/Ry/eV):'
741 0 : sum_test=Ex_PBE+Ec_LDA+Ecnl
742 0 : write(jnlout,'(3F24.15)') sum_test,sum_test*2.0,sum_test*Ha2eV
743 : !
744 0 : write(jnlout,'(/,A)') '-------------- revPBE exchange --------------'
745 0 : write(jnlout,'(A)') 'Ex_revPBE (Ha/Ry/eV):'
746 0 : write(jnlout,'(3F24.15)') Ex_revPBE,Ex_revPBE*2.0,Ex_revPBE*Ha2eV
747 : !
748 0 : write(jnlout,'(A)') 'Ex_revPBE+Ec_LDA+Ecnl (Ha/Ry/eV):'
749 0 : sum_test=Ex_revPBE+Ec_LDA+Ecnl
750 0 : write(jnlout,'(3F24.15)') sum_test,sum_test*2.0,sum_test*Ha2eV
751 : !
752 0 : if (trim(option) == 'vdw2' .or. trim(option) == 'vdW2') then
753 : !
754 : ! in the case of vdW-DF2 use refit PW86 exchange
755 : !
756 0 : write(jnlout,'(/,A)') '-------------- for vdW-DF2 --------------'
757 0 : write(jnlout,'(A)') 'refit PW86 exchange:'
758 0 : write(jnlout,'(A)') 'Ex_PW86 (Ha/Ry/eV):'
759 0 : write(jnlout,'(3F24.15)') Ex_PW86,Ex_PW86*2.0,Ex_PW86*Ha2eV
760 : !
761 0 : write(jnlout,'(A)') 'Ex_PW86+Ec_LDA+Ecnl (Ha/Ry/eV):'
762 0 : sum_test=Ex_PW86+Ec_LDA+Ecnl
763 0 : write(jnlout,'(3F24.15)') sum_test,sum_test*2.0,sum_test*Ha2eV
764 : !
765 : endif
766 : !
767 0 : call calc_ehartree(n)
768 : !
769 0 : deallocate(grad_n,G,G_ind)
770 0 : deallocate(e_cNL,e_cSL,e_cLDA)
771 0 : deallocate(q_0,theta_alpha,u_a)
772 : !
773 0 : END SUBROUTINE soler
774 : !
775 : !==========================================================================================
776 : !
777 : ! We take the kernel from the Thonhauser implementation (Espresso). This subroutine to
778 : ! read the kernel thus is also taken from there.
779 :
780 0 : SUBROUTINE read_kernel()
781 : USE param, ONLY: pi,ftable,jnlout
782 : USE nonlocal_data,ONLY: n_alpha,n_k, &
783 : r_max,q_cut,dk, &
784 : q_alpha,phi,d2phi_dk2
785 : implicit none
786 :
787 : character(len=30) :: double_format = '(1p4e23.14)'
788 : integer :: q1, q2
789 0 : if (allocated(q_alpha)) return
790 0 : write(jnlout,*) 'reading kernel table'
791 :
792 0 : open(ftable, file='vdW_kernel_table')
793 :
794 0 : read(ftable, '(2i5)' ) n_alpha, n_k
795 0 : read(ftable, double_format) r_max
796 :
797 0 : allocate(q_alpha(n_alpha), phi(0:n_k,n_alpha,n_alpha), d2phi_dk2(0:n_k,n_alpha,n_alpha))
798 :
799 0 : read(ftable, double_format) q_alpha
800 :
801 0 : q_cut = q_alpha(n_alpha)
802 :
803 0 : write(jnlout,*)'n_a:', n_alpha
804 0 : write(jnlout,*)'q_c:', q_cut
805 0 : write(jnlout,*)
806 :
807 0 : do q1 = 1, n_alpha
808 0 : do q2 = 1, q1
809 :
810 0 : read(ftable, double_format ) phi(0:n_k, q1, q2)
811 0 : phi(0:n_k,q2, q1) = phi(0:n_k, q1, q2)
812 :
813 : end do
814 : end do
815 :
816 0 : do q1 = 1, n_alpha
817 0 : do q2 = 1, q1
818 :
819 0 : read(ftable, double_format ) d2phi_dk2(0:n_k,q1, q2)
820 0 : d2phi_dk2(0:n_k,q2, q1) = d2phi_dk2(0:n_k,q1, q2)
821 :
822 : end do
823 : end do
824 :
825 0 : dk = 2.0*pi/r_max
826 : !
827 0 : close(ftable)
828 : !
829 : END SUBROUTINE read_kernel
830 : !
831 : !==========================================================================================
832 : !
833 0 : SUBROUTINE setup_g_vectors()
834 : USE param, ONLY: jnlout
835 : USE nonlocal_data,ONLY: nx,ny,nz,n_gvectors,n_grid, &
836 : G_cut,b1,b2,b3, &
837 : G,G_ind
838 : implicit none
839 : !
840 : integer :: igx,igy,igz
841 : integer :: tmpx,tmpy,tmpz
842 : integer :: ind
843 : !
844 : real :: G_squ
845 : real :: gx, gy, gz
846 :
847 : ! We need the G vectors in a number of subroutines of this code. So the idea is to set them up
848 : ! once and for all to minimize possible sources of errors. The array G(n_gvectors,3) will hold
849 : ! the components of the g_vectors and G_ind their index on the fft mesh. We have to do the loop
850 : ! over all G vectors twice. First time to count number of G_vectors inside cut off sphere and
851 : ! second time to set them up.
852 :
853 0 : n_gvectors = 0
854 :
855 0 : do igz = -(nz-1)/2,(nz-1)/2
856 0 : do igy = -(ny-1)/2,(ny-1)/2
857 0 : do igx = -(nx-1)/2,(nx-1)/2
858 :
859 0 : gx = igx * b1(1) + igy * b2(1) + igz * b3(1)
860 0 : gy = igx * b1(2) + igy * b2(2) + igz * b3(2)
861 0 : gz = igx * b1(3) + igy * b2(3) + igz * b3(3)
862 :
863 0 : G_squ = gx**2 + gy**2 + gz**2
864 :
865 0 : if ( G_squ .le. G_cut ) n_gvectors = n_gvectors + 1
866 :
867 : enddo
868 : enddo
869 : enddo
870 0 : if (allocated(g)) return
871 0 : allocate(G(n_gvectors, 3), G_ind(n_gvectors))
872 :
873 0 : write(jnlout,*) "#g-vectors:",n_gvectors," outof ",nz*ny*nx
874 :
875 0 : G = 0.0
876 0 : G_ind = 0
877 :
878 0 : n_gvectors = 0
879 :
880 0 : do igz = -(nz-1)/2,(nz-1)/2
881 0 : do igy = -(ny-1)/2,(ny-1)/2
882 0 : do igx = -(nx-1)/2,(nx-1)/2
883 :
884 0 : gx = igx * b1(1) + igy * b2(1) + igz * b3(1)
885 0 : gy = igx * b1(2) + igy * b2(2) + igz * b3(2)
886 0 : gz = igx * b1(3) + igy * b2(3) + igz * b3(3)
887 :
888 0 : G_squ = gx**2 + gy**2 + gz**2
889 :
890 0 : tmpx = 0
891 0 : tmpy = 0
892 0 : tmpz = 0
893 :
894 0 : IF (igx .LT. 0) tmpx = nx
895 0 : IF (igy .LT. 0) tmpy = ny
896 0 : IF (igz .LT. 0) tmpz = nz
897 :
898 : ind = 1 + ( igx + tmpx ) + &
899 : nx*( igy + tmpy ) + &
900 0 : nx*ny*( igz + tmpz )
901 :
902 0 : if ( G_squ .le. G_cut ) then
903 :
904 0 : n_gvectors = n_gvectors + 1
905 :
906 0 : G(n_gvectors, 1) = gx
907 0 : G(n_gvectors, 2) = gy
908 0 : G(n_gvectors, 3) = gz
909 :
910 0 : G_ind(n_gvectors) = ind
911 :
912 : endif
913 :
914 : enddo
915 : enddo
916 : enddo
917 :
918 : end subroutine setup_g_vectors
919 : !
920 : !==========================================================================================
921 : !
922 0 : SUBROUTINE calc_gradient(n,grad_n)
923 : USE param, ONLY: jnlout
924 : USE nonlocal_data,ONLY: n_grid,n_gvectors,nx,ny,nz, &
925 : a1,a2,a3,G,G_ind
926 : USE driver_fft
927 : !
928 : implicit none
929 : !
930 : real, intent(in) :: n(n_grid)
931 : real, intent(inout) :: grad_n(n_grid,3)
932 :
933 0 : complex, allocatable :: n_cmplx(:)
934 0 : complex, allocatable :: grad_n_cmplx(:,:)
935 :
936 : integer :: i
937 :
938 : integer :: ind
939 :
940 : integer*8 :: plan
941 :
942 0 : grad_n = 0.0
943 :
944 : ! A derivative d/dx f(x) in real space corresponds to i*G*f(G) in reciprocal space. So we have
945 : ! to fourier transform the density, find the according G vectors, calculate the product and back
946 : ! fourier transform the gradient
947 :
948 0 : allocate(n_cmplx(n_grid), grad_n_cmplx(n_grid,3))
949 :
950 0 : n_cmplx = (0.0, 0.0)
951 0 : grad_n_cmplx = (0.0, 0.0)
952 :
953 0 : do i = 1, n_grid
954 :
955 0 : n_cmplx(i) = cmplx( n(i), 0.0)
956 :
957 : enddo
958 :
959 0 : call inplfft( n_cmplx, -1 )
960 :
961 0 : do i = 1, n_gvectors
962 :
963 0 : ind = G_ind(i)
964 :
965 0 : grad_n_cmplx(ind,:) = (0.0,1.0) * G(i,:) * n_cmplx(ind)
966 :
967 : enddo
968 :
969 0 : do i=1,3
970 :
971 0 : call inplfft( grad_n_cmplx(:,i), 1)
972 :
973 : enddo
974 : !
975 0 : grad_n(:,:) = real(grad_n_cmplx(:,:))
976 :
977 0 : deallocate(n_cmplx, grad_n_cmplx)
978 : !
979 0 : END SUBROUTINE calc_gradient
980 : !
981 : !==========================================================================================
982 : !
983 : ! This subroutine calculates q_0 following eqns. 11 and 12 of Dion and saturates it
984 : ! following eqn. 7 of Soler
985 : !
986 0 : SUBROUTINE calc_q0(n, grad_n, q_0)
987 : USE param, ONLY: pi, &
988 : A,alpha_1,beta_1,beta_2,beta_2,beta_3,beta_4 ! Parameters for LDA_c
989 : USE nonlocal_data,ONLY: n_grid,q_cut,m_c,Zab
990 : !
991 : implicit none
992 : !
993 : real, intent(in) :: n(n_grid) ! charge density
994 : real, intent(in) :: grad_n(n_grid,3) ! gradient of charge density
995 : real, intent(out) :: q_0(n_grid) ! array holding g_0 on the grid
996 :
997 : real :: q ! q before saturation
998 :
999 : real :: k_F ! fermi wave vector
1000 : real :: r_s !
1001 : real :: sqrt_r_s !
1002 : real :: eLDA_c ! LDA correlation
1003 : real :: eLDA_x ! LDA exchange
1004 : real :: grad_n_squ ! squared gradient of n
1005 :
1006 : real :: LDA_dummy_1
1007 : real :: LDA_dummy_2
1008 :
1009 : real :: sum_m
1010 :
1011 : integer :: i, m ! counters
1012 : !
1013 : ! loop over all grid points
1014 : !
1015 0 : q_0(:) = q_cut
1016 : !
1017 0 : do i=1,n_grid
1018 :
1019 : ! check if charge density is negative. If so treat it like zero charge density. Zero
1020 : ! density corresponds to high q_0. Thats why we set it to q_cut the highest possible
1021 : ! q_0 and continue with the next grid point.
1022 :
1023 0 : if ( n(i) < 1e-12 ) then
1024 0 : q_0(i) = q_cut
1025 0 : cycle
1026 : end if
1027 : !
1028 : ! calculate q according to eqns (11) and (12) from Dion.
1029 : !
1030 0 : k_F = (3.0*(pi**2.0)*n(i))**(1.0/3.0)
1031 0 : r_s = (3.0/(4.0*pi*n(i)))**(1.0/3.0)
1032 0 : sqrt_r_s = r_s**(1.0/2.0)
1033 :
1034 0 : grad_n_squ = grad_n(i,1)**2.0 + grad_n(i,2)**2.0 + grad_n(i,3)**2.0
1035 :
1036 0 : LDA_dummy_1 = (8.0*pi/3.0)*A*(1.0+ alpha_1*r_s)
1037 : LDA_dummy_2 = 2.0*A*(beta_1*sqrt_r_s + beta_2*r_s + &
1038 0 : beta_3*sqrt_r_s*r_s + beta_4*r_s*r_s)
1039 :
1040 0 : eLDA_x = (-3.0/(4.0*pi))*k_F
1041 0 : eLDA_c = (-3.0/(4.0*pi))*LDA_dummy_1*log(1.0+1.0/LDA_dummy_2)
1042 : !
1043 : ! LDA correlation energy density is epsilon_c[n(r)]*n(r)
1044 : !
1045 : ! e_cLDA(i) = eLDA_c*n(i)
1046 : !
1047 : q = (1.0 + ( eLDA_c / eLDA_x) - (Zab/9.0)*grad_n_squ &
1048 0 : /(4.0*(n(i)**2.0)*(k_F**2.0)))*k_F
1049 :
1050 : ! now saturate q according to eq. (7) from Soler.
1051 :
1052 0 : sum_m = 0.0
1053 :
1054 0 : do m=1,m_c
1055 :
1056 0 : sum_m = sum_m + (q/q_cut)**m / real(m)
1057 :
1058 : enddo
1059 :
1060 0 : q_0(i) = q_cut*(1.0 - exp(-sum_m))
1061 :
1062 : enddo
1063 :
1064 0 : END SUBROUTINE calc_q0
1065 : !
1066 : !==========================================================================================
1067 : !
1068 : !
1069 0 : SUBROUTINE calc_theta_i(n, q_alpha, q_0, theta_alpha)
1070 :
1071 : USE nonlocal_data,ONLY: n_grid, n_alpha
1072 : !
1073 : implicit none
1074 : !
1075 : real, intent(in) :: n(n_grid) ! charge density
1076 : real, intent(in) :: q_alpha(n_alpha) ! q-mesh for interpolation
1077 : real, intent(in) :: q_0(n_grid) ! q values on the grid
1078 :
1079 : complex, intent(inout) :: theta_alpha(n_grid,n_alpha) ! thetas on real space grid
1080 :
1081 : integer :: i ! counter
1082 :
1083 : ! this call will interpolate dirac delta which gives p_alpha(q_i) according to a method
1084 : ! from numerical recipes. Whithin this subroutine also the initial setup of the second
1085 : ! derivatives needed for spline interpolation will be done once.
1086 :
1087 0 : call splint( q_alpha, q_0, theta_alpha )
1088 :
1089 0 : do i=1,n_grid
1090 :
1091 : ! theta_alpha will hold at this stage the p_alpha(q_i)
1092 :
1093 0 : theta_alpha(i,:) = n(i)*theta_alpha(i,:)
1094 :
1095 : enddo
1096 :
1097 0 : END SUBROUTINE calc_theta_i
1098 : !
1099 : !==========================================================================================
1100 : !
1101 0 : SUBROUTINE calc_ecnl(Ecnl, theta_alpha, u_a, e_cNL)
1102 : USE param ,ONLY: pi,jnlout
1103 : USE nonlocal_data,ONLY: nx,ny,nz,n_grid,n_alpha,n_gvectors, &
1104 : omega,time1,time2, &
1105 : G, G_ind
1106 : USE driver_fft
1107 : !
1108 : implicit none
1109 : !
1110 : real, intent(out) :: Ecnl
1111 : complex,intent(in) :: theta_alpha(n_grid,n_alpha)
1112 : complex,intent(out) :: u_a(n_grid,n_alpha) ! function u_a(r) needed for the construction
1113 : ! of the potential and also for non local
1114 : real, intent(inout) :: e_cNL(n_grid,n_alpha) ! non local correlation energy density ! correlation energy density.
1115 : !
1116 : integer :: i, ind
1117 : integer :: alpha,beta
1118 : real :: k
1119 0 : real :: phi_k(n_alpha,n_alpha)
1120 : complex :: integral
1121 : !
1122 0 : u_a = (0.0, 0.0)
1123 0 : Ecnl = 0.0
1124 : !
1125 : ! integration of theta_a*theta_b*phi_ab
1126 : !
1127 0 : write (jnlout,*) 'calculating E_c,nl:'
1128 : !
1129 0 : call DATE_AND_TIME(values=time1)
1130 : !
1131 : ! The difference to older versions of this code is the way how the G-vectors are set up. The
1132 : ! index is not any more calculated on the fly but stored in G_ind.
1133 : !
1134 0 : do i = 1, n_gvectors
1135 :
1136 0 : ind = G_ind(i)
1137 :
1138 0 : k = sqrt(dble( G(i,1)**2 + G(i,2)**2 + G(i,3)**2))
1139 :
1140 0 : call interpolate_kernel(k, phi_k)
1141 :
1142 0 : integral = (0.0, 0.0)
1143 :
1144 0 : do alpha = 1,n_alpha
1145 0 : do beta = 1,n_alpha
1146 :
1147 : integral = integral + conjg(theta_alpha(ind,alpha)) * &
1148 0 : theta_alpha(ind,beta)*cmplx(phi_k(alpha,beta),0.0) !&
1149 : !
1150 : ! the array u_a(k,a) = sum_b theta_b(k) * phi_ab(k) equals the fourier transform
1151 : ! of the function u_a(r) = sum_b int d^3 r' theta_b(r') phi_ab(r - r') which we need for
1152 : ! the calculation of the potential and the non local correlation energy density.
1153 : !
1154 : u_a(ind,alpha) = u_a(ind,alpha) + &
1155 0 : theta_alpha(ind,beta)*cmplx(phi_k(alpha,beta),0.0)
1156 :
1157 : enddo
1158 : enddo
1159 :
1160 0 : Ecnl = Ecnl + real(integral)
1161 :
1162 : enddo
1163 : !
1164 : ! backward fourier transform of u_a(k,a) in order to get u_a(r) the non local correlation
1165 : ! energy density.
1166 : !
1167 0 : do i=1,n_alpha
1168 :
1169 0 : call inplfft( u_a(:,i), 1 )
1170 :
1171 : enddo
1172 : !
1173 : ! at this moment e_cNL(i,alpha) holds the theta_alpha_i which we have to multiply by the convolution
1174 : ! u_a(r) in order to get the non local energy density as a function of alpha
1175 : !
1176 0 : e_cNL = e_cNL*real(u_a)
1177 : !
1178 0 : call DATE_AND_TIME(values=time2)
1179 : !
1180 0 : write(jnlout,'(A)') 'time for integration:'
1181 0 : call timing( time2 - time1 )
1182 : !
1183 : ! taking care of the units. (2*pi)^3/omega is the volume element of the reciprocal unit
1184 : ! cell. omega^2 we get when replacing the integrals by sums over a real space grid. The
1185 : ! (2*pi)^3 cancels with an according factor from the radial fourier transform on the kernel
1186 : !
1187 0 : Ecnl = Ecnl*0.5*omega
1188 : !
1189 0 : END SUBROUTINE calc_ecnl
1190 : !
1191 : !==========================================================================================
1192 : !
1193 : SUBROUTINE make_q_alpha(q_alpha)
1194 :
1195 :
1196 : USE nonlocal_data,ONLY: n_alpha,q_cut,lambda ! n_alpha: number of q points, q_cut: maximum
1197 : ! q value, lambda: parameter for logarithmic
1198 : ! mesh.
1199 : implicit none
1200 : !
1201 : real, intent(out) :: q_alpha(n_alpha) ! array holding the qs
1202 :
1203 : real :: q1 ! auxiliary vairable holding the first
1204 : ! q which can be calculated explicitely
1205 :
1206 : integer :: i ! counter
1207 :
1208 : ! first we have to setup the qs starting from the maximum value q_cut on a logarithmic mesh
1209 : ! which satisfies ( q_a+1 - q_a ) = lambda*( q_a - q_a-1 ). Lambda > 1 is a parameter.
1210 : ! In GPAW lambda = 1.2 is used.
1211 :
1212 : q1 = q_cut * ( lambda - 1.0) / (lambda**( n_alpha - 1.0 ) - 1.0 )
1213 :
1214 : do i = 1,n_alpha
1215 :
1216 : q_alpha(i) = q1 * (lambda**( i - 1.0 ) - 1.0 ) / ( lambda - 1.0 )
1217 :
1218 : enddo
1219 :
1220 : END SUBROUTINE make_q_alpha
1221 : !
1222 : !==========================================================================================
1223 : !
1224 : ! From Numerical Recipes
1225 :
1226 0 : SUBROUTINE splint( x_i, x, p_iofx )
1227 :
1228 :
1229 : USE nonlocal_data,ONLY: n_grid, n_alpha
1230 : !
1231 : implicit none
1232 : !
1233 : real, intent(in) :: x_i(n_alpha) ! the x_i values where the function is known
1234 : real, intent(in) :: x(n_grid) ! the x values for which the function
1235 : ! shall be interpolated
1236 :
1237 : complex, intent(inout) :: p_iofx(n_grid,n_alpha) ! the function values for each x
1238 :
1239 0 : real, allocatable :: d2y_dx2(:,:) ! second derivatives needed for the splines
1240 :
1241 0 : real, allocatable :: y(:)
1242 :
1243 : real :: a, b, c, d, h ! some intermediate variables for the interpolation
1244 :
1245 : integer :: i, j ! counters
1246 : integer :: u_bound, l_bound ! variables for finding the section of x_i
1247 : ! in which x is located.
1248 : integer :: alpha ! index of the found section
1249 :
1250 0 : allocate(y(n_alpha))
1251 :
1252 0 : allocate(d2y_dx2(n_alpha,n_alpha))
1253 :
1254 0 : call setup_spline(x_i,d2y_dx2)
1255 :
1256 0 : do i=1,n_grid
1257 :
1258 : ! first find the correct section of x_i in which x is located by bisectioning. According to
1259 : ! numerical recipes this is efficient if the x values are random. In our case there might
1260 : ! be some correlation as the density and its gradient are smooth.
1261 :
1262 : l_bound = 1
1263 : u_bound = n_alpha
1264 :
1265 0 : do while ((u_bound - l_bound) > 1)
1266 :
1267 0 : alpha = (u_bound + l_bound) / 2
1268 :
1269 0 : if ( x(i) > x_i(alpha)) then
1270 : l_bound = alpha
1271 : else
1272 0 : u_bound = alpha
1273 : endif
1274 : enddo
1275 :
1276 0 : h = x_i(u_bound) - x_i(l_bound)
1277 :
1278 0 : a = (x_i(u_bound) - x(i))/h
1279 0 : b = (x(i) - x_i(l_bound))/h
1280 0 : c = ((a**3.0 - a)*h**2.0)/6.0
1281 0 : d = ((b**3.0 - b)*h**2.0)/6.0
1282 :
1283 0 : do j=1,n_alpha
1284 :
1285 0 : y = 0.0
1286 0 : y(j) = 1.0
1287 :
1288 : p_iofx(i,j) = a*y(l_bound) + b*y(u_bound) + &
1289 0 : c*d2y_dx2(j,l_bound) + d*d2y_dx2(j,u_bound)
1290 :
1291 : enddo
1292 : enddo
1293 :
1294 0 : deallocate(y, d2y_dx2)
1295 :
1296 0 : END SUBROUTINE splint
1297 : !
1298 : !==========================================================================================
1299 : !
1300 : ! From Numerical Recipes
1301 :
1302 0 : SUBROUTINE setup_spline(x_i,d2y_dx2)
1303 :
1304 :
1305 : USE nonlocal_data,ONLY: n_alpha
1306 : !
1307 : implicit none
1308 : !
1309 : real, intent(in) :: x_i(n_alpha)
1310 : real, intent(inout) :: d2y_dx2(n_alpha,n_alpha)
1311 :
1312 0 : real, allocatable :: y(:) ! this array holds the function values at x_i which
1313 : ! are going to be interpolated.
1314 0 : real, allocatable :: dy_dx(:) ! temporary array holding the first derivative
1315 :
1316 : real :: sig, p ! temporary values for the interpolation
1317 :
1318 : integer :: i, j ! counter
1319 :
1320 0 : allocate(y(n_alpha), dy_dx(n_alpha))
1321 :
1322 0 : do i=1,n_alpha
1323 :
1324 : ! as we are interpolating dirac delta set the according function values:
1325 :
1326 0 : y=0.0
1327 0 : y(i)=1.0
1328 :
1329 : ! now according to numerical recipes we set the boundary conditions which will give
1330 : ! so called "natural" splines.
1331 :
1332 0 : d2y_dx2(i,1) = 0.0
1333 0 : dy_dx(1) = 0.0
1334 :
1335 0 : do j=2,n_alpha-1
1336 :
1337 0 : sig = (x_i(j) - x_i(j-1)) / (x_i(j+1) - x_i(j-1))
1338 0 : p = sig*d2y_dx2(i,j-1) + 2.0
1339 0 : d2y_dx2(i,j) = (sig - 1.0)/p
1340 0 : dy_dx(j) = (y(j+1) - y(j))/(x_i(j+1) - x_i(j)) - (y(j) - y(j-1))/(x_i(j) - x_i(j-1))
1341 0 : dy_dx(j) = (6.0*dy_dx(j)/(x_i(j+1) - x_i(j-1)) - sig*dy_dx(j-1))/p
1342 :
1343 : enddo
1344 :
1345 0 : d2y_dx2(i,n_alpha) = 0.0
1346 :
1347 0 : do j=n_alpha - 1, 1, -1
1348 :
1349 0 : d2y_dx2(i,j) = d2y_dx2(i,j)*d2y_dx2(i,j+1) + dy_dx(j)
1350 :
1351 : enddo
1352 : enddo
1353 :
1354 0 : deallocate(y, dy_dx)
1355 :
1356 0 : END SUBROUTINE setup_spline
1357 : !
1358 : !===========================================================================================
1359 : !
1360 : ! similar to the Thonhauser implementation
1361 :
1362 0 : SUBROUTINE interpolate_kernel(k, phi_k)
1363 :
1364 :
1365 : USE nonlocal_data,ONLY: n_alpha, dk, n_k, phi, d2phi_dk2
1366 : !
1367 : implicit none
1368 : !
1369 : real, intent(in) :: k
1370 : real, intent(inout) :: phi_k(n_alpha,n_alpha)
1371 :
1372 : real :: a, b, c, d
1373 :
1374 : integer :: q1, q2, k_i
1375 :
1376 : ! the algorithm for interpolation will be more or less the same as in splint().
1377 0 : phi_k = 0.0
1378 :
1379 : ! find the interval in which our k lies
1380 :
1381 0 : k_i = int(k/dk)
1382 :
1383 : ! simple case when k equals one of the values we have tabulated the kernel for
1384 :
1385 0 : if (mod(k,dk) == 0.0) then
1386 :
1387 0 : do q1=1,n_alpha
1388 0 : do q2=1,q1
1389 :
1390 0 : phi_k(q1, q2) = phi(k_i, q1, q2)
1391 0 : phi_k(q2, q1) = phi(k_i, q2, q1)
1392 :
1393 : enddo
1394 : enddo
1395 :
1396 : return
1397 :
1398 : endif
1399 :
1400 0 : a = (dk*(k_i+1.0) - k)/dk
1401 0 : b = (k - dk*k_i)/dk
1402 0 : c = (a**3.0-a)*dk**2.0/6.0
1403 0 : d = (b**3.0-b)*dk**2.0/6.0
1404 :
1405 0 : do q1 = 1, n_alpha
1406 0 : do q2 = 1, q1
1407 :
1408 : phi_k(q1, q2) = a*phi(k_i, q1, q2) + b*phi(k_i+1, q1, q2) &
1409 0 : +(c*d2phi_dk2(k_i, q1, q2) + d*d2phi_dk2(k_i+1,q1, q2))
1410 :
1411 0 : phi_k(q2, q1) = phi_k(q1, q2)
1412 :
1413 : end do
1414 : end do
1415 :
1416 : END SUBROUTINE interpolate_kernel
1417 : !
1418 : !=================================================================================================
1419 : !
1420 : #ifdef POTENTIAL
1421 :
1422 0 : subroutine calc_dq0(n, grad_n, dq0_dn, dq0_dgrad_n)
1423 :
1424 : USE param, ONLY: pi, &
1425 : A,alpha_1,beta_1,beta_2,beta_2,beta_3,beta_4 ! Parameters for LDA_c
1426 : USE nonlocal_data,ONLY: n_grid,q_cut,m_c,Zab
1427 :
1428 :
1429 : implicit none
1430 :
1431 : real, intent(in) :: n(n_grid) ! charge density
1432 : real, intent(in) :: grad_n(n_grid,3) ! gradient of charge density
1433 :
1434 : real, intent(out) :: dq0_dn(n_grid) ! derivative of q0 w.r.t. the density
1435 : real, intent(out) :: dq0_dgrad_n(n_grid) ! derivative of q0 w.r.t. the gradient
1436 :
1437 : real :: q ! q before saturation
1438 :
1439 : real :: dq0_dq ! derivative needed for dq0_dn and dq0_dgrad_n
1440 :
1441 : real :: k_F ! fermi wave vector
1442 : real :: r_s !
1443 : real :: sqrt_r_s !
1444 : real :: eLDA_c ! LDA correlation
1445 : real :: eLDA_x ! LDA exchange
1446 : real :: grad_n_squ ! squared gradient of n
1447 :
1448 : real :: LDA_dummy_1
1449 : real :: LDA_dummy_2
1450 :
1451 : real :: sum_m
1452 :
1453 : integer :: i, m ! counters
1454 :
1455 : ! loop over all grid points
1456 :
1457 0 : dq0_dn(:) = 0.0
1458 0 : dq0_dgrad_n(:) = 0.0
1459 :
1460 0 : do i=1,n_grid
1461 :
1462 0 : if ( n(i) < 1e-12 ) then
1463 : cycle ! NOTE: The derivatives will be zero at these points
1464 : end if
1465 : !
1466 : ! calculate q according to eqns (11) and (12) from Dion.
1467 : !
1468 0 : k_F = (3.0*(pi**2.0)*n(i))**(1.0/3.0)
1469 0 : r_s = (3.0/(4.0*pi*n(i)))**(1.0/3.0)
1470 0 : sqrt_r_s = r_s**(1.0/2.0)
1471 :
1472 0 : grad_n_squ = grad_n(i,1)**2.0 + grad_n(i,2)**2.0 + grad_n(i,3)**2.0
1473 :
1474 0 : LDA_dummy_1 = (8.0*pi/3.0)*A*(1.0+ alpha_1*r_s)
1475 : LDA_dummy_2 = 2.0*A*(beta_1*sqrt_r_s + beta_2*r_s + &
1476 0 : beta_3*sqrt_r_s*r_s + beta_4*r_s*r_s)
1477 :
1478 0 : eLDA_x = (-3.0/(4.0*pi))*k_F
1479 0 : eLDA_c = (-3.0/(4.0*pi))*LDA_dummy_1*log(1.0+1.0/LDA_dummy_2)
1480 :
1481 : q = (1.0 + ( eLDA_c / eLDA_x) - (Zab/9.0)*grad_n_squ &
1482 0 : /(4.0*(n(i)**2.0)*(k_F**2.0)))*k_F
1483 :
1484 : ! now saturate q according to eq. (7) from Soler.
1485 :
1486 0 : sum_m = 0.0
1487 0 : dq0_dq = 0.0
1488 :
1489 0 : do m=1,m_c
1490 :
1491 0 : sum_m = sum_m + (q/q_cut)**m / real(m)
1492 :
1493 0 : dq0_dq = dq0_dq + ((q/q_cut)**(m-1))
1494 :
1495 : enddo
1496 :
1497 0 : dq0_dq = dq0_dq * exp(-sum_m)
1498 :
1499 : ! here we calculate the derivatives needed for the potential later. i.e. we calculate
1500 : ! n*dq0/dn and n*dq0/dgrad_n so that we do not have to divide by n which might be very
1501 : ! small at some points and thus produce numerical errors.
1502 :
1503 : dq0_dn(i) = dq0_dq * ( k_F/3.0 &
1504 : + k_F*7.0/3.0*(Zab/9.0)*grad_n_squ &
1505 : /(4.0*(n(i)**2.0)*(k_F**2.0)) &
1506 : - (8.0*pi/9.0)*A*alpha_1*r_s*log(1.0+1.0/LDA_dummy_2) &
1507 : + LDA_dummy_1/(LDA_dummy_2*(1.0 + LDA_dummy_2)) &
1508 : * (2.0*A*(beta_1/6.0*sqrt_r_s + beta_2/3.0*r_s &
1509 0 : + beta_3/2.0*r_s*sqrt_r_s + 2.0*beta_4/3.0*r_s**2)))
1510 :
1511 : dq0_dgrad_n(i) = dq0_dq * 2.0 * (-1.0*Zab/9.0) &
1512 0 : * sqrt(grad_n_squ) / (4.0*n(i)*k_F)
1513 :
1514 : enddo
1515 :
1516 0 : end subroutine calc_dq0
1517 : !
1518 : !=======================================================================================
1519 : !
1520 : ! similar to the Thonhauser implementation
1521 :
1522 0 : subroutine calc_potential(v_nl, u_a, q_0, dq0_dn, dq0_dgrad_n, n, grad_n)
1523 :
1524 :
1525 : use nonlocal_data, only : n_grid, n_alpha, q_alpha, q_cut, nx, ny, nz, G, G_ind, n_gvectors
1526 :
1527 : USE driver_fft
1528 :
1529 : implicit none
1530 :
1531 : real, intent(out) :: v_nl(n_grid) ! non local part of the potential
1532 :
1533 : real, intent(in) :: q_0(n_grid)
1534 :
1535 : real, intent(in) :: dq0_dn(n_grid) ! n*dq0/dn
1536 : real, intent(in) :: dq0_dgrad_n(n_grid) ! n*dq0/dgrad_n
1537 :
1538 : real, intent(in) :: n(n_grid)
1539 : real, intent(in) :: grad_n(n_grid,3)
1540 :
1541 : complex, intent(in) :: u_a(n_grid, n_alpha)
1542 :
1543 0 : complex, allocatable :: h_j(:,:)
1544 :
1545 0 : real, allocatable :: d2y_dx2(:,:) ! second derivatives needed for the splines
1546 :
1547 0 : real, allocatable :: y(:)
1548 :
1549 : real :: a, b, c, d, d1, d2, h ! some intermediate variables
1550 : ! for the interpolation of p_a
1551 :
1552 : real :: p_a
1553 : real :: dpa_dq0
1554 :
1555 : real :: grad_n_abs ! | grad_n |
1556 :
1557 : integer :: i, alpha
1558 :
1559 : integer :: l_bound, u_bound
1560 :
1561 : integer :: ind
1562 :
1563 0 : allocate(h_j(n_grid,3))
1564 0 : allocate(y(n_alpha))
1565 :
1566 0 : h_j = (0.0, 0.0)
1567 :
1568 0 : allocate( d2y_dx2(n_alpha, n_alpha) )
1569 :
1570 0 : call setup_spline( q_alpha, d2y_dx2 )
1571 :
1572 : ! The potential will be calculated using FFT following White and Bird. PRB 50, 4954
1573 : ! (1994) eqn. (10). v^xc(r) = d f_xc / d n(r) - (1/N) * sum_G,r' i * G * ( grad_n(r')
1574 : ! / |grad_n(r')| ) * ( d f_xc / d |grad_n(r')| ) * e^(i*G*(r - r')) with E_c,NL =
1575 : ! int d^3 r f_xc .
1576 :
1577 0 : do i=1,n_grid
1578 :
1579 : ! first we have to interpolate the polynomials p_a and their derivatives with respect
1580 : ! to q0 again as we need them for the calculation of the potential. This will be the
1581 : ! same procedure as in splint.
1582 :
1583 : l_bound = 1
1584 : u_bound = n_alpha
1585 :
1586 0 : do while ((u_bound - l_bound) > 1)
1587 :
1588 0 : ind = (u_bound + l_bound) / 2
1589 :
1590 0 : if ( q_0(i) > q_alpha(ind)) then
1591 : l_bound = ind
1592 : else
1593 0 : u_bound = ind
1594 : endif
1595 : enddo
1596 :
1597 0 : h = q_alpha(u_bound) - q_alpha(l_bound)
1598 :
1599 0 : a = (q_alpha(u_bound) - q_0(i))/h
1600 0 : b = (q_0(i) - q_alpha(l_bound))/h
1601 0 : c = ((a**3.0 - a)*h**2.0)/6.0
1602 0 : d = ((b**3.0 - b)*h**2.0)/6.0
1603 0 : d1 = (3.0*a**2.0 - 1.0)*h/6.0
1604 0 : d2 = (3.0*b**2.0 - 1.0)*h/6.0
1605 :
1606 0 : do alpha = 1, n_alpha
1607 :
1608 0 : y = 0.0
1609 0 : y(alpha) = 1.0
1610 :
1611 : p_a = a*y(l_bound) + b*y(u_bound) + &
1612 0 : c*d2y_dx2(alpha,l_bound) + d*d2y_dx2(alpha,u_bound)
1613 :
1614 : dpa_dq0 = (y(u_bound) - y(l_bound))/h &
1615 0 : - d1*d2y_dx2(alpha,l_bound) + d2*d2y_dx2(alpha,u_bound)
1616 :
1617 : ! first term sum_a u_ai * ( p_ai + n_i * dpai/dqi * dqi/dni ). the factor n(i) is already
1618 : ! contained in dq0_dn
1619 :
1620 0 : v_nl(i) = v_nl(i) + u_a(i,alpha) * (p_a + dpa_dq0 * dq0_dn(i) )
1621 :
1622 : ! the following sum we need for h_j which will be fourier transformed later. The IF
1623 : ! condition excludes the contributions of high q values.
1624 :
1625 0 : grad_n_abs = sqrt(grad_n(i,1)**2.0 + grad_n(i,2)**2.0 + grad_n(i,3)**2.0)
1626 :
1627 0 : if ( q_0(i) .ne. q_cut ) then
1628 :
1629 : h_j(i,:) = h_j(i,:) + u_a(i,alpha)*cmplx((grad_n(i,:) / grad_n_abs) &
1630 0 : * dpa_dq0 * dq0_dgrad_n(i), 0.0)
1631 :
1632 : endif
1633 : enddo
1634 : enddo
1635 :
1636 : ! now fourier transform h_j and carry out the sum over G.
1637 :
1638 0 : do i = 1,3
1639 :
1640 0 : call inplfft( h_j(:,i), -1 )
1641 :
1642 : enddo
1643 :
1644 0 : do i = 1, n_gvectors
1645 :
1646 0 : ind = G_ind(i)
1647 :
1648 0 : h_j(ind,:) = (0.0,1.0) * G(i,:) * h_j(ind,:)
1649 :
1650 : enddo
1651 :
1652 : ! back fourier transform h_j and add it to the potential
1653 :
1654 0 : do i = 1,3
1655 :
1656 0 : call inplfft( h_j(:,i), 1 )
1657 :
1658 : enddo
1659 :
1660 0 : v_nl(:) = v_nl(:) - real( h_j(:,1) + h_j(:,2) + h_j(:,3) )
1661 :
1662 0 : deallocate(h_j, y,d2y_dx2)
1663 :
1664 0 : end subroutine
1665 : !
1666 : !==========================================================================================
1667 : !
1668 : #endif
1669 :
1670 0 : SUBROUTINE timing ( time )
1671 : USE param,ONLY: jnlout
1672 : !
1673 : implicit none
1674 : !
1675 : integer, intent(in) :: time(8)
1676 : integer :: tmp(8), time_out(8)
1677 :
1678 : integer :: i
1679 :
1680 : ! In the older versions of this program the time written to output file could be negative. This
1681 : ! subroutine will fix this problem.
1682 :
1683 0 : tmp = 0
1684 :
1685 0 : tmp(8) = 1000
1686 0 : tmp(7) = 60
1687 0 : tmp(6) = 60
1688 0 : tmp(5) = 24
1689 :
1690 0 : do i=1,8
1691 :
1692 0 : if (time(i) .lt. 0) then
1693 0 : time_out(i) = time(i) + tmp(i)
1694 0 : time_out(i-1) = time(i-1) - 1
1695 : else
1696 0 : time_out(i) = time(i)
1697 : endif
1698 :
1699 : enddo
1700 : !
1701 0 : write(jnlout, '(I4,A,I4,A,I4,A,I4,A)') time_out(5)," h ", &
1702 0 : time_out(6)," min", &
1703 0 : time_out(7)," s ", &
1704 0 : time_out(8)," ms "
1705 : !
1706 0 : END SUBROUTINE timing
1707 : !
1708 : ! As the hole thing is intended to be interfaced with different codes I put here just an
1709 : ! example program for testing, which is capable of reading an input file, reading a
1710 : ! charge density, which otherwise has to be given by the calling program, and then call
1711 : ! the subroutine soler(n, ... ).
1712 : !
1713 : END MODULE nonlocal_funct
|
