« back to all changes in this revision
Viewing changes to Src/xc.f
- browse files at revision 1051
- compare with another revision
- download diff
- download tarball
- view history from revision 1051
- Committer: Alberto Garcia
- Date: 2018-12-18 11:38:26 UTC
- Revision ID: albertog@icmab.es-20181218113826-q1nw0qwlc5smc0xe
Simplify the choice between SiestaXC and BSC's versions of cellxc
The BSC_CELLXC preprocessor blocks have been removed and replaced
by conditional blocks which can be selected at runtime.
The user can control which version of cellxc is used by means of the
fdf logical variable:
XC.Use.BSC.CellXC (T/F)
which is 'false' by default, thus using SiestaXC's cellxc.
The BSC version might be slightly better for GGA operations, and the
SiestaXC version must of course be used when dealing with van der Waals
density functionals.
NOTES:
* BSC's version of cellxc now uses ldaxc, ggaxc, and setxc/getxc from
SiestaXC. This enhances its functionality (more functionals) and
saves on duplicated code. The old 'xc.f' file has been removed.
* In 'meshsubs', 'distriphionmesh' now accepts an extra argument to
flag the need for stencil initilization.
* In 'forhar', a number of 'reord' operations towards 'sequential'
fine-point ordering (which were compiled-in only for the new
interface) have been removed, as all the arrays involved are already
in 'sequential' form.
Tests show that these changes do not seem to affect the results, however.
The 'forhar' interface uses the 'ntm' argument in all cases, to simplify
the code.
* Removed bogus incompatibility check in ldau.F
The BSC_CELLXC preprocessor blocks have been removed and replaced
by conditional blocks which can be selected at runtime.
The user can control which version of cellxc is used by means of the
fdf logical variable:
XC.Use.BSC.CellXC (T/F)
which is 'false' by default, thus using SiestaXC's cellxc.
The BSC version might be slightly better for GGA operations, and the
SiestaXC version must of course be used when dealing with van der Waals
density functionals.
NOTES:
* BSC's version of cellxc now uses ldaxc, ggaxc, and setxc/getxc from
SiestaXC. This enhances its functionality (more functionals) and
saves on duplicated code. The old 'xc.f' file has been removed.
* In 'meshsubs', 'distriphionmesh' now accepts an extra argument to
flag the need for stencil initilization.
* In 'forhar', a number of 'reord' operations towards 'sequential'
fine-point ordering (which were compiled-in only for the new
interface) have been removed, as all the arrays involved are already
in 'sequential' form.
Tests show that these changes do not seem to affect the results, however.
The 'forhar' interface uses the 'ntm' argument in all cases, to simplify
the code.
* Removed bogus incompatibility check in ldau.F
added
removed
Src/xc.f
1
!
2
! Copyright (C) 1996-2016 The SIESTA group
3
! This file is distributed under the terms of the
4
! GNU General Public License: see COPYING in the top directory
5
! or http://www.gnu.org/copyleft/gpl.txt.
6
! See Docs/Contributors.txt for a list of contributors.
7
!
8
! *******************************************************************
9
! This file contains XC subroutines used when siesta is compiled with
10
! option BSC_CELLXC. Otherwise, the SiestaXC library is used.
11
! *******************************************************************
12
13
subroutine atomxc( IREL, NR, MAXR, RMESH, nspin, Dens,
14
. EX, EC, DX, DC, VXC )
15
16
C *******************************************************************
17
C Finds total exchange-correlation energy and potential for a
18
C spherical electron density distribution.
19
C This version implements the Local (spin) Density Approximation and
20
C the Generalized-Gradient-Aproximation with the 'explicit mesh
21
C functional' approach of White & Bird, PRB 50, 4954 (1994).
22
C Gradients are 'defined' by numerical derivatives, using 2*NN+1 mesh
23
C points, where NN is a parameter defined below
24
C Ref: L.C.Balbas et al, PRB 64, 165110 (2001)
25
C Wrtten by J.M.Soler using algorithms developed by
26
C L.C.Balbas, J.L.Martins and J.M.Soler, Dec.1996
27
C ************************* INPUT ***********************************
28
C CHARACTER*(*) FUNCTL : Functional to be used:
29
C 'LDA' or 'LSD' => Local (spin) Density Approximation
30
C 'GGA' => Generalized Gradient Corrections
31
C Uppercase is optional
32
C CHARACTER*(*) AUTHOR : Parametrization desired:
33
C 'CA' or 'PZ' => LSD Perdew & Zunger, PRB 23, 5075 (1981)
34
C 'PW91' => GGA Perdew & Wang, JCP, 100, 1290 (1994)
35
C 'PW92' => LSD Perdew & Wang, PRB, 45, 13244 (1992). This is
36
C the local density limit of the next:
37
C 'PBE' => GGA Perdew, Burke & Ernzerhof, PRL 77, 3865 (1996)
38
C 'RPBE' => GGA Hammer, Hansen & Norskov, PRB 59, 7413 (1999)
39
C 'REVPBE' => GGA Zhang & Yang, PRL 80,890(1998)
40
C 'LYP' => GGA Becke-Lee-Yang-Parr (see subroutine blypxc)
41
C 'WC' => GGA Wu-Cohen (see subroutine wcxc)
42
C 'PBESOL' => GGA Perdew et al, PRL, 100, 136406 (2008)
43
C Uppercase is optional
44
C INTEGER IREL : Relativistic exchange? (0=>no, 1=>yes)
45
C INTEGER NR : Number of radial mesh points
46
C INTEGER MAXR : Physical first dimension of RMESH, Dens and VXC
47
C REAL*8 RMESH(MAXR) : Radial mesh points
48
C INTEGER nspin : nspin=1 => unpolarized; nspin=2 => polarized
49
C REAL*8 Dens(MAXR,nspin) : Total (nspin=1) or spin (nspin=2) electron
50
C density at mesh points
51
C ************************* OUTPUT **********************************
52
C REAL*8 EX : Total exchange energy
53
C REAL*8 EC : Total correlation energy
54
C REAL*8 DX : IntegralOf( rho * (eps_x - v_x) )
55
C REAL*8 DC : IntegralOf( rho * (eps_c - v_c) )
56
C REAL*8 VXC(MAXR,nspin) : (Spin) exch-corr potential
57
C ************************ UNITS ************************************
58
C Distances in atomic units (Bohr).
59
C Densities in atomic units (electrons/Bohr**3)
60
C Energy unit depending of parameter EUNIT below
61
C ********* ROUTINES CALLED *****************************************
62
C GGAXC, LDAXC
63
C *******************************************************************
64
65
use precision, only : dp
66
use bsc_xcmod, only : nXCfunc, XCfunc, XCauth
67
use bsc_xcmod, only : XCweightX, XCweightC
68
use sys, only: die
69
use alloc, only: re_alloc, de_alloc
70
71
C Next line is nonstandard but may be suppressed
72
implicit none
73
74
C Argument types and dimensions
75
integer, intent(in) :: IREL
76
integer, intent(in) :: MAXR
77
integer, intent(in) :: NR
78
integer, intent(in) :: nspin
79
real(dp), intent(in) :: Dens(MAXR,nspin)
80
real(dp), intent(in) :: RMESH(MAXR)
81
real(dp), intent(out) :: VXC(MAXR,nspin)
82
real(dp), intent(out) :: DC
83
real(dp), intent(out) :: DX
84
real(dp), intent(out) :: EC
85
real(dp), intent(out) :: EX
86
87
C Internal parameters
88
C NN : order of the numerical derivatives: the number of radial
89
C points used is 2*NN+1
90
C mspin : must be equal or larger than nspin (4 for non-collinear spin)
91
integer, parameter :: mspin = 4
92
integer, parameter :: NN = 5
93
94
C Fix energy unit: EUNIT=1.0 => Hartrees,
95
C EUNIT=0.5 => Rydbergs,
96
C EUNIT=0.03674903 => eV
97
real(dp), parameter :: EUNIT = 0.5_dp
98
99
C DVMIN is added to differential of volume to avoid division by zero
100
real(dp), parameter :: DVMIN = 1.0d-12
101
102
C Local variables and arrays
103
logical
104
. GGA, GGAfunc
105
integer
106
. IN, IN1, IN2, IR, IS, JN, NF
107
real(dp)
108
. D(mspin), DECDD(mspin), DECDGD(3,mspin),
109
. DEXDD(mspin), DEXDGD(3,mspin),
110
. DGDM(-NN:NN), DGIDFJ(-NN:NN), DRDM, DVol,
111
. DVCDN(mspin,mspin), DVXDN(mspin,mspin),
112
. EPSC, EPSX, F1, F2, GD(3,mspin), PI
113
real(dp), pointer :: Aux(:)
114
external
115
. GGAXC, LDAXC
116
117
C Set GGA switch
118
GGA = .false.
119
do nf = 1,nXCfunc
120
if ( XCfunc(nf).eq.'GGA' .or. XCfunc(nf).eq.'gga') then
121
GGA = .true.
122
else
123
if ( XCfunc(nf).ne.'LDA' .and. XCfunc(nf).ne.'lda' .and.
124
. XCfunc(nf).ne.'LSD' .and. XCfunc(nf).ne.'lsd' ) then
125
call die('ATOMXC: Unknown functional ' // XCfunc(nf))
126
endif
127
endif
128
enddo
129
130
C Initialize output
131
EX = 0.0_dp
132
EC = 0.0_dp
133
DX = 0.0_dp
134
DC = 0.0_dp
135
do IS = 1,nspin
136
do IR = 1,NR
137
VXC(IR,IS) = 0.0_dp
138
enddo
139
enddo
140
141
C Set up workspace array
142
if (GGA) then
143
nullify( Aux )
144
call re_alloc( AUX, 1, NR, 'AUX', 'atomxc' )
145
endif
146
147
C Get number pi
148
PI = 4.0_dp * ATAN(1.0_dp)
149
150
C Loop on mesh points
151
do IR = 1,NR
152
153
C Find interval of neighbour points to calculate derivatives
154
IN1 = MAX( 1, IR-NN ) - IR
155
IN2 = MIN( NR, IR+NN ) - IR
156
157
C Find weights of numerical derivation from Lagrange
158
C interpolation formula
159
do IN = IN1,IN2
160
IF (IN .EQ. 0) THEN
161
DGDM(IN) = 0
162
do JN = IN1,IN2
163
IF (JN.NE.0) DGDM(IN) = DGDM(IN) + 1.D0 / (0 - JN)
164
enddo
165
ELSE
166
F1 = 1
167
F2 = 1
168
do JN = IN1,IN2
169
IF (JN.NE.IN .AND. JN.NE.0) F1 = F1 * (0 - JN)
170
IF (JN.NE.IN) F2 = F2 * (IN - JN)
171
enddo
172
DGDM(IN) = F1 / F2
173
ENDIF
174
enddo
175
176
C Find dr/dmesh
177
DRDM = 0.0_dp
178
do IN = IN1,IN2
179
DRDM = DRDM + RMESH(IR+IN) * DGDM(IN)
180
enddo
181
182
C Find differential of volume. Use trapezoidal integration rule
183
DVol = 4.0_dp * PI * RMESH(IR)**2 * DRDM
184
C DVMIN is a small number added to avoid a division by zero
185
DVol = DVol + DVMIN
186
if (IR.eq.1 .or. IR.eq.NR) DVol = 0.5_dp*DVol
187
if (GGA) Aux(IR) = DVol
188
189
C Find the weights for the derivative d(gradF(i))/d(F(j)), of
190
C the gradient at point i with respect to the value at point j
191
if (GGA) then
192
do IN = IN1,IN2
193
DGIDFJ(IN) = DGDM(IN) / DRDM
194
enddo
195
endif
196
197
C Find density and gradient of density at this point
198
do IS = 1,nspin
199
D(IS) = Dens(IR,IS)
200
enddo
201
if (GGA) then
202
do IS = 1,nspin
203
GD(1,IS) = 0.0_dp
204
GD(2,IS) = 0.0_dp
205
GD(3,IS) = 0.0_dp
206
do IN = IN1,IN2
207
GD(3,IS) = GD(3,IS) + DGIDFJ(IN) * Dens(IR+IN,IS)
208
enddo
209
enddo
210
endif
211
212
C Loop over exchange-correlation functions
213
do nf = 1,nXCfunc
214
215
C Is this a GGA?
216
if (XCfunc(nf).eq.'GGA' .or. XCfunc(nf).eq.'gga') then
217
GGAfunc = .true.
218
else
219
GGAfunc = .false.
220
endif
221
222
C Find exchange and correlation energy densities and their
223
C derivatives with respect to density and density gradient
224
if (GGAfunc) then
225
CALL GGAXC( XCauth(nf), IREL, nspin, D, GD,
226
. EPSX, EPSC, DEXDD, DECDD, DEXDGD, DECDGD )
227
else
228
CALL LDAXC( XCauth(nf), IREL, nspin, D, EPSX, EPSC, DEXDD,
229
. DECDD, DVXDN, DVCDN )
230
endif
231
232
C Scale terms by weights
233
EPSX = EPSX*XCweightX(nf)
234
EPSC = EPSC*XCweightC(nf)
235
DEXDD(1:nspin) = DEXDD(1:nspin)*XCweightX(nf)
236
DECDD(1:nspin) = DECDD(1:nspin)*XCweightC(nf)
237
if (GGAfunc) then
238
DEXDGD(1:3,1:nspin) = DEXDGD(1:3,1:nspin)*XCweightX(nf)
239
DECDGD(1:3,1:nspin) = DECDGD(1:3,1:nspin)*XCweightC(nf)
240
endif
241
242
C Add contributions to exchange-correlation energy and its
243
C derivatives with respect to density at all points
244
do IS = 1,nspin
245
EX = EX + DVol*D(IS)*EPSX
246
EC = EC + DVol*D(IS)*EPSC
247
DX = DX + DVol*D(IS)*(EPSX - DEXDD(IS))
248
DC = DC + DVol*D(IS)*(EPSC - DECDD(IS))
249
if (GGAfunc) then
250
VXC(IR,IS) = VXC(IR,IS) + DVol*(DEXDD(IS) + DECDD(IS))
251
do IN = IN1,IN2
252
DX= DX - DVol*Dens(IR+IN,IS)*DEXDGD(3,IS)*DGIDFJ(IN)
253
DC= DC - DVol*Dens(IR+IN,IS)*DECDGD(3,IS)*DGIDFJ(IN)
254
VXC(IR+IN,IS) = VXC(IR+IN,IS) +
255
. DVol*(DEXDGD(3,IS) + DECDGD(3,IS))*DGIDFJ(IN)
256
enddo
257
else
258
if (GGA) then
259
VXC(IR,IS) = VXC(IR,IS) + DVol*(DEXDD(IS) + DECDD(IS))
260
else
261
VXC(IR,IS) = VXC(IR,IS) + DEXDD(IS) + DECDD(IS)
262
endif
263
endif
264
enddo
265
266
enddo
267
268
enddo
269
270
C Divide by volume element to obtain the potential (per electron)
271
if (GGA) then
272
do IS = 1,NSPIN
273
do IR = 1,NR
274
DVol = AUX(IR)
275
VXC(IR,IS) = VXC(IR,IS) / DVol
276
enddo
277
enddo
278
call de_alloc( AUX, 'AUX', 'atomxc' )
279
endif
280
281
C Divide by energy unit
282
EX = EX / EUNIT
283
EC = EC / EUNIT
284
DX = DX / EUNIT
285
DC = DC / EUNIT
286
do IS = 1,nspin
287
do IR = 1,NR
288
VXC(IR,IS) = VXC(IR,IS) / EUNIT
289
enddo
290
enddo
291
292
end
293
294
295
subroutine exchng( IREL, NSP, DS, EX, VX )
296
297
C *****************************************************************
298
C Finds local exchange energy density and potential
299
C Adapted by J.M.Soler from routine velect of Froyen's
300
C pseudopotential generation program. Madrid, Jan'97. Version 0.5.
301
C **** Input ******************************************************
302
C INTEGER IREL : relativistic-exchange switch (0=no, 1=yes)
303
C INTEGER NSP : spin-polarizations (1=>unpolarized, 2=>polarized)
304
C REAL*8 DS(NSP) : total (nsp=1) or spin (nsp=2) electron density
305
C **** Output *****************************************************
306
C REAL*8 EX : exchange energy density
307
C REAL*8 VX(NSP) : (spin-dependent) exchange potential
308
C **** Units ******************************************************
309
C Densities in electrons/Bohr**3
310
C Energies in Hartrees
311
C *****************************************************************
312
313
use precision, only: dp
314
implicit none
315
316
integer, intent(in) :: nsp, irel
317
real(dp), intent(in) :: DS(NSP)
318
real(dp), intent(out) :: VX(NSP)
319
real(dp), intent(out) :: EX
320
321
real(dp), parameter :: zero = 0.0_dp, one = 1.0_dp
322
real(dp), parameter :: pfive = 0.5_dp, opf = 1.5_dp
323
real(dp), parameter :: C014 = 0.014_dp
324
325
real(dp) :: pi, trd, ftrd, tftm, a0
326
real(dp) :: alp, d1, d2, d, z, fz, fzp, rs, vxp, exp_var
327
real(dp) :: beta, sb, vxf, exf, alb
328
329
PI=4*ATAN(ONE)
330
TRD = ONE/3
331
FTRD = 4*TRD
332
TFTM = 2**FTRD-2
333
A0 = (4/(9*PI))**TRD
334
335
C X-alpha parameter:
336
ALP = 2 * TRD
337
338
IF (NSP .EQ. 2) THEN
339
D1 = MAX(DS(1),0.D0)
340
D2 = MAX(DS(2),0.D0)
341
D = D1 + D2
342
IF (D .LE. ZERO) THEN
343
EX = ZERO
344
VX(1) = ZERO
345
VX(2) = ZERO
346
RETURN
347
ENDIF
348
Z = (D1 - D2) / D
349
FZ = ((1+Z)**FTRD+(1-Z)**FTRD-2)/TFTM
350
FZP = FTRD*((1+Z)**TRD-(1-Z)**TRD)/TFTM
351
ELSE
352
D = DS(1)
353
IF (D .LE. ZERO) THEN
354
EX = ZERO
355
VX(1) = ZERO
356
RETURN
357
ENDIF
358
Z = ZERO
359
FZ = ZERO
360
FZP = ZERO
361
ENDIF
362
RS = (3 / (4*PI*D) )**TRD
363
VXP = -(3*ALP/(2*PI*A0*RS))
364
EXP_VAR = 3*VXP/4
365
IF (IREL .EQ. 1) THEN
366
BETA = C014/RS
367
SB = SQRT(1+BETA*BETA)
368
ALB = LOG(BETA+SB)
369
VXP = VXP * (-PFIVE + OPF * ALB / (BETA*SB))
370
EXP_VAR = EXP_VAR * (ONE-OPF*((BETA*SB-ALB)/BETA**2)**2)
371
ENDIF
372
VXF = 2**TRD*VXP
373
EXF = 2**TRD*EXP_VAR
374
IF (NSP .EQ. 2) THEN
375
VX(1) = VXP + FZ*(VXF-VXP) + (1-Z)*FZP*(EXF-EXP_VAR)
376
VX(2) = VXP + FZ*(VXF-VXP) - (1+Z)*FZP*(EXF-EXP_VAR)
377
EX = EXP_VAR + FZ*(EXF-EXP_VAR)
378
ELSE
379
VX(1) = VXP
380
EX = EXP_VAR
381
ENDIF
382
END
383
384
SUBROUTINE GGAXC( AUTHOR, IREL, nspin, D, GD,
385
. EPSX, EPSC, DEXDD, DECDD, DEXDGD, DECDGD )
386
387
C Finds the exchange and correlation energies at a point, and their
388
C derivatives with respect to density and density gradient, in the
389
C Generalized Gradient Correction approximation.
390
C Lengths in Bohr, energies in Hartrees
391
C Written by L.C.Balbas and J.M.Soler, Dec'96. Version 0.5.
392
C Modified by V.M.Garcia-Suarez to include non-collinear spin. June 2002
393
394
use precision, only : dp
395
use sys, only : die
396
397
implicit none
398
399
CHARACTER*(*) AUTHOR
400
INTEGER IREL, nspin, NS, IS, IX
401
real(dp) THETA, PHI, D(nspin), DECDD(nspin),
402
. DECDGD(3,nspin), DEXDD(nspin), DEXDGD(3,nspin),
403
. EPSC, EPSX, GD(3,nspin),
404
. DD(2), DTOT, DPOL,
405
. GDD(3,2), TINY, DECDN(2), DEXDN(2),
406
. VPOL, DECDGN(3,2), DEXDGN(3,2),
407
. C2, S2, ST, CT, CP, SP, dpolz, dpolxy
408
409
410
PARAMETER ( TINY = 1.D-12 )
411
412
IF (nspin .EQ. 4) THEN
413
C Find eigenvalues of density matrix (up and down densities
414
C along the spin direction)
415
C Note: D(1)=D11, D(2)=D22, D(3)=Real(D12), D(4)=Im(D12)
416
NS = 2
417
DTOT = D(1) + D(2)
418
419
! Explicit calculation of the rotation-matrix elements from
420
! the entries of D
421
422
dpolz= D(1)-D(2)
423
dpolxy= 2.0d0*sqrt(D(3)**2+D(4)**2)
424
DPOL = sqrt( dpolz**2 + dpolxy**2 )
425
if ( DPOL.gt.1.0d-12 ) then
426
THETA = atan2(dpolxy,dpolz)
427
else
428
THETA = 0.0_dp
429
endif
430
C2 = COS(THETA/2)
431
S2 = SIN(THETA/2)
432
ST = SIN(THETA)
433
CT = COS(THETA)
434
PHI = ATAN2(-D(4),D(3))
435
CP = COS(PHI)
436
SP = SIN(PHI)
437
438
DD(1) = 0.5D0 * ( DTOT + DPOL )
439
DD(2) = 0.5D0 * ( DTOT - DPOL )
440
441
C Find diagonal elements of the gradient
442
DO IX = 1,3
443
GDD(IX,1) = GD(IX,1)*C2**2 + GD(IX,2)*S2**2 +
444
. 2.d0*C2*S2*(GD(IX,3)*CP - GD(IX,4)*SP)
445
GDD(IX,2) = GD(IX,1)*S2**2 + GD(IX,2)*C2**2 -
446
. 2.d0*C2*S2*(GD(IX,3)*CP - GD(IX,4)*SP)
447
ENDDO
448
ELSE
449
NS = nspin
450
DO 20 IS = 1,nspin
451
cag Avoid negative densities
452
DD(IS) = max(D(IS),0.0d0)
453
DO 30 IX = 1,3
454
GDD(IX,IS) = GD(IX,IS)
455
30 CONTINUE
456
20 CONTINUE
457
ENDIF
458
459
IF (AUTHOR.EQ.'PBE' .OR. AUTHOR.EQ.'pbe') THEN
460
CALL PBEXC( IREL, NS, DD, GDD,
461
. EPSX, EPSC, DEXDN, DECDN, DEXDGN, DECDGN )
462
cmvfs
463
ELSE IF (AUTHOR.EQ.'RPBE'.OR.AUTHOR.EQ.'rpbe') THEN
464
CALL RPBEXC( IREL, NS, DD, GDD,
465
. EPSX, EPSC, DEXDN, DECDN, DEXDGN, DECDGN )
466
cmvfs
467
ELSE IF (AUTHOR.EQ.'WC'.OR.AUTHOR.EQ.'wc') THEN
468
CALL WCXC( IREL, NS, DD, GDD,
469
. EPSX, EPSC, DEXDN, DECDN, DEXDGN, DECDGN )
470
cea
471
ELSE IF (AUTHOR.EQ.'REVPBE'.OR.AUTHOR.EQ.'revpbe'
472
. .OR.AUTHOR.EQ.'revPBE') THEN
473
CALL REVPBEXC( IREL, NS, DD, GDD,
474
. EPSX, EPSC, DEXDN, DECDN, DEXDGN, DECDGN )
475
cag
476
ELSE IF (AUTHOR.EQ.'LYP'.OR.AUTHOR.EQ.'lyp') THEN
477
CALL BLYPXC( NS, DD, GDD, EPSX, EPSC, dEXdn, dECdn,
478
. DEXDGN, DECDGN)
479
cag
480
ELSEIF (AUTHOR.EQ.'PW91' .OR. AUTHOR.EQ.'pw91') THEN
481
CALL PW91XC( IREL, NS, DD, GDD,
482
. EPSX, EPSC, DEXDN, DECDN, DEXDGN, DECDGN )
483
cjdg
484
ELSEIF (AUTHOR.EQ.'PBESOL'.OR.AUTHOR.EQ.'pbesol'
485
. .OR.AUTHOR.EQ.'PBEsol') THEN
486
CALL PBESOLXC( IREL, NS, DD, GDD,
487
. EPSX, EPSC, DEXDN, DECDN, DEXDGN, DECDGN )
488
ELSE
489
call die('GGAXC: Unknown author ' // trim(AUTHOR))
490
ENDIF
491
492
IF (nspin .EQ. 4) THEN
493
C Find dE/dD(ispin) = dE/dDup * dDup/dD(ispin) +
494
C dE/dDdown * dDown/dD(ispin)
495
C Note convention:
496
C DEDD(1)=dE/dD11, DEDD(2)=dE/dD22,
497
C DEDD(3)=Re(dE/dD12)=Re(dE/dD21),
498
C DEDD(4)=Im(dE/dD12)=-Im(dE/D21)
499
C
500
VPOL = (DEXDN(1)-DEXDN(2)) * CT
501
DEXDD(1) = 0.5D0 * ( DEXDN(1) + DEXDN(2) + VPOL )
502
DEXDD(2) = 0.5D0 * ( DEXDN(1) + DEXDN(2) - VPOL )
503
DEXDD(3) = 0.5d0 * (DEXDN(1)-DEXDN(2)) * ST * CP
504
DEXDD(4) =-0.5d0 * (DEXDN(1)-DEXDN(2)) * ST * SP
505
VPOL = (DECDN(1)-DECDN(2)) * CT
506
DECDD(1) = 0.5D0 * ( DECDN(1) + DECDN(2) + VPOL )
507
DECDD(2) = 0.5D0 * ( DECDN(1) + DECDN(2) - VPOL )
508
DECDD(3) = 0.5d0 * (DECDN(1)-DECDN(2)) * ST * CP
509
DECDD(4) =-0.5d0 * (DECDN(1)-DECDN(2)) * ST * SP
510
C Gradient terms
511
DO 40 IX = 1,3
512
DEXDGD(IX,1) = DEXDGN(IX,1)*C2**2 + DEXDGN(IX,2)*S2**2
513
DEXDGD(IX,2) = DEXDGN(IX,1)*S2**2 + DEXDGN(IX,2)*C2**2
514
DEXDGD(IX,3) = 0.5D0*(DEXDGN(IX,1) - DEXDGN(IX,2))*ST*CP
515
DEXDGD(IX,4) =-0.5D0*(DEXDGN(IX,1) - DEXDGN(IX,2))*ST*SP
516
DECDGD(IX,1) = DECDGN(IX,1)*C2**2 + DECDGN(IX,2)*S2**2
517
DECDGD(IX,2) = DECDGN(IX,1)*S2**2 + DECDGN(IX,2)*C2**2
518
DECDGD(IX,3) = 0.5D0*(DECDGN(IX,1) - DECDGN(IX,2))*ST*CP
519
DECDGD(IX,4) =-0.5D0*(DECDGN(IX,1) - DECDGN(IX,2))*ST*SP
520
40 CONTINUE
521
ELSE
522
DO 60 IS = 1,nspin
523
DEXDD(IS) = DEXDN(IS)
524
DECDD(IS) = DECDN(IS)
525
DO 50 IX = 1,3
526
DEXDGD(IX,IS) = DEXDGN(IX,IS)
527
DECDGD(IX,IS) = DECDGN(IX,IS)
528
50 CONTINUE
529
60 CONTINUE
530
ENDIF
531
532
END
533
534
535
SUBROUTINE LDAXC( AUTHOR, IREL, nspin, D, EPSX, EPSC, VX, VC,
536
. DVXDN, DVCDN )
537
538
C ******************************************************************
539
C Finds the exchange and correlation energies and potentials, in the
540
C Local (spin) Density Approximation.
541
C Written by L.C.Balbas and J.M.Soler, Dec'96.
542
C Non-collinear spin added by J.M.Soler, May'98
543
C *********** INPUT ************************************************
544
C CHARACTER*(*) AUTHOR : Parametrization desired:
545
C 'CA' or 'PZ' => LSD Perdew & Zunger, PRB 23, 5075 (1981)
546
C 'PW92' => LSD Perdew & Wang, PRB, 45, 13244 (1992)
547
C Uppercase is optional
548
C INTEGER IREL : Relativistic exchange? (0=>no, 1=>yes)
549
C INTEGER nspin : nspin=1 => unpolarized; nspin=2 => polarized;
550
C nspin=4 => non-collinear polarization
551
C REAL*8 D(nspin) : Local (spin) density. For non-collinear
552
C polarization, the density matrix is given by:
553
C D(1)=D11, D(2)=D22, D(3)=Real(D12), D(4)=Im(D12)
554
C *********** OUTPUT ***********************************************
555
C REAL*8 EPSX, EPSC : Exchange and correlation energy densities
556
C REAL*8 VX(nspin), VC(nspin) : Exchange and correlation potentials,
557
C defined as dExc/dD(ispin)
558
C REAL*8 DVXDN(nspin,nspin) : Derivative of exchange potential with
559
C respect the charge density, defined
560
C as DVx(spin1)/Dn(spin2)
561
C REAL*8 DVCDN(nspin,nspin) : Derivative of correlation potential
562
C respect the charge density, defined
563
C as DVc(spin1)/Dn(spin2)
564
C *********** UNITS ************************************************
565
C Lengths in Bohr, energies in Hartrees
566
C ******************************************************************
567
568
use precision, only : dp
569
use sys, only : die
570
571
implicit none
572
573
CHARACTER*(*) AUTHOR
574
INTEGER IREL, nspin
575
real(dp) D(nspin), EPSC, EPSX, VX(nspin), VC(nspin),
576
. DVXDN(nspin,nspin), DVCDN(nspin,nspin)
577
578
INTEGER IS, NS, ISPIN1, ISPIN2
579
real(dp) DD(2), DPOL, DTOT, TINY, VCD(2), VPOL, VXD(2)
580
581
PARAMETER ( TINY = 1.D-12 )
582
583
IF (nspin .EQ. 4) THEN
584
C Find eigenvalues of density matrix (up and down densities
585
C along the spin direction)
586
C Note: D(1)=D11, D(2)=D22, D(3)=Real(D12), D(4)=Im(D12)
587
NS = 2
588
DTOT = D(1) + D(2)
589
DPOL = SQRT( (D(1)-D(2))**2 + 4.D0*(D(3)**2+D(4)**2) )
590
DD(1) = 0.5D0 * ( DTOT + DPOL )
591
DD(2) = 0.5D0 * ( DTOT - DPOL )
592
ELSE
593
NS = nspin
594
DO 10 IS = 1,nspin
595
cag Avoid negative densities
596
DD(IS) = max(D(IS),0.0d0)
597
10 CONTINUE
598
ENDIF
599
600
601
DO ISPIN2 = 1, nspin
602
DO ISPIN1 = 1, nspin
603
DVXDN(ISPIN1,ISPIN2) = 0.D0
604
DVCDN(ISPIN1,ISPIN2) = 0.D0
605
ENDDO
606
ENDDO
607
608
IF ( AUTHOR.EQ.'CA' .OR. AUTHOR.EQ.'ca' .OR.
609
. AUTHOR.EQ.'PZ' .OR. AUTHOR.EQ.'pz') THEN
610
CALL PZXC( IREL, NS, DD, EPSX, EPSC, VXD, VCD, DVXDN, DVCDN )
611
ELSEIF ( AUTHOR.EQ.'PW92' .OR. AUTHOR.EQ.'pw92' ) THEN
612
CALL PW92XC( IREL, NS, DD, EPSX, EPSC, VXD, VCD )
613
ELSE
614
call die('LDAXC: Unknown author ' // trim(AUTHOR))
615
ENDIF
616
617
IF (nspin .EQ. 4) THEN
618
C Find dE/dD(ispin) = dE/dDup * dDup/dD(ispin) +
619
C dE/dDdown * dDown/dD(ispin)
620
VPOL = (VXD(1)-VXD(2)) * (D(1)-D(2)) / (DPOL+TINY)
621
VX(1) = 0.5D0 * ( VXD(1) + VXD(2) + VPOL )
622
VX(2) = 0.5D0 * ( VXD(1) + VXD(2) - VPOL )
623
VX(3) = (VXD(1)-VXD(2)) * D(3) / (DPOL+TINY)
624
VX(4) = (VXD(1)-VXD(2)) * D(4) / (DPOL+TINY)
625
VPOL = (VCD(1)-VCD(2)) * (D(1)-D(2)) / (DPOL+TINY)
626
VC(1) = 0.5D0 * ( VCD(1) + VCD(2) + VPOL )
627
VC(2) = 0.5D0 * ( VCD(1) + VCD(2) - VPOL )
628
VC(3) = (VCD(1)-VCD(2)) * D(3) / (DPOL+TINY)
629
VC(4) = (VCD(1)-VCD(2)) * D(4) / (DPOL+TINY)
630
ELSE
631
DO 20 IS = 1,nspin
632
VX(IS) = VXD(IS)
633
VC(IS) = VCD(IS)
634
20 CONTINUE
635
ENDIF
636
END
637
638
639
SUBROUTINE PBEXC( IREL, nspin, Dens, GDens,
640
. EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
641
642
C *********************************************************************
643
C Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation.
644
C Ref: J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996)
645
C Written by L.C.Balbas and J.M.Soler. December 1996. Version 0.5.
646
C ******** INPUT ******************************************************
647
C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
648
C INTEGER nspin : Number of spin polarizations (1 or 2)
649
C REAL*8 Dens(nspin) : Total electron density (if nspin=1) or
650
C spin electron density (if nspin=2)
651
C REAL*8 GDens(3,nspin) : Total or spin density gradient
652
C ******** OUTPUT *****************************************************
653
C REAL*8 EX : Exchange energy density
654
C REAL*8 EC : Correlation energy density
655
C REAL*8 DEXDD(nspin) : Partial derivative
656
C d(DensTot*Ex)/dDens(ispin),
657
C where DensTot = Sum_ispin( Dens(ispin) )
658
C For a constant density, this is the
659
C exchange potential
660
C REAL*8 DECDD(nspin) : Partial derivative
661
C d(DensTot*Ec)/dDens(ispin),
662
C where DensTot = Sum_ispin( Dens(ispin) )
663
C For a constant density, this is the
664
C correlation potential
665
C REAL*8 DEXDGD(3,nspin): Partial derivative
666
C d(DensTot*Ex)/d(GradDens(i,ispin))
667
C REAL*8 DECDGD(3,nspin): Partial derivative
668
C d(DensTot*Ec)/d(GradDens(i,ispin))
669
C ********* UNITS ****************************************************
670
C Lengths in Bohr
671
C Densities in electrons per Bohr**3
672
C Energies in Hartrees
673
C Gradient vectors in cartesian coordinates
674
C ********* ROUTINES CALLED ******************************************
675
C EXCHNG, PW92C
676
C ********************************************************************
677
678
use precision, only : dp
679
680
implicit none
681
INTEGER IREL, nspin
682
real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
683
. DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
684
685
C Internal variables
686
INTEGER
687
. IS, IX
688
689
real(dp)
690
. A, BETA, D(2), DADD, DECUDD, DENMIN,
691
. DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
692
. DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
693
. DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
694
. DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
695
. EC, ECUNIF, EX, EXUNIF,
696
. F, F1, F2, F3, F4, FC, FX, FOUTHD,
697
. GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
698
. H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
699
. T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
700
701
C Lower bounds of density and its gradient to avoid divisions by zero
702
PARAMETER ( DENMIN = 1.D-12 )
703
PARAMETER ( GDMIN = 1.D-12 )
704
705
C Fix some numerical parameters
706
PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
707
. THD=1.D0/3.D0, THRHLF=1.5D0,
708
. TWO=2.D0, TWOTHD=2.D0/3.D0 )
709
710
C Fix some more numerical constants
711
PI = 4 * ATAN(1.D0)
712
BETA = 0.066725D0
713
GAMMA = (1 - LOG(TWO)) / PI**2
714
MU = BETA * PI**2 / 3
715
KAPPA = 0.804D0
716
717
C Translate density and its gradient to new variables
718
IF (nspin .EQ. 1) THEN
719
D(1) = HALF * Dens(1)
720
D(2) = D(1)
721
DT = MAX( DENMIN, Dens(1) )
722
DO 10 IX = 1,3
723
GD(IX,1) = HALF * GDens(IX,1)
724
GD(IX,2) = GD(IX,1)
725
GDT(IX) = GDens(IX,1)
726
10 CONTINUE
727
ELSE
728
D(1) = Dens(1)
729
D(2) = Dens(2)
730
DT = MAX( DENMIN, Dens(1)+Dens(2) )
731
DO 20 IX = 1,3
732
GD(IX,1) = GDens(IX,1)
733
GD(IX,2) = GDens(IX,2)
734
GDT(IX) = GDens(IX,1) + GDens(IX,2)
735
20 CONTINUE
736
ENDIF
737
GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
738
GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
739
GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 )
740
GDMT = MAX( GDMIN, GDMT )
741
742
C Find local correlation energy and potential
743
CALL PW92C( 2, D, ECUNIF, VCUNIF )
744
745
C Find total correlation energy
746
RS = ( 3 / (4*PI*DT) )**THD
747
KF = (3 * PI**2 * DT)**THD
748
KS = SQRT( 4 * KF / PI )
749
ZETA = ( D(1) - D(2) ) / DT
750
ZETA = MAX( -1.D0+DENMIN, ZETA )
751
ZETA = MIN( 1.D0-DENMIN, ZETA )
752
PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
753
T = GDMT / (2 * PHI * KS * DT)
754
F1 = ECUNIF / GAMMA / PHI**3
755
F2 = EXP(-F1)
756
A = BETA / GAMMA / (F2-1)
757
F3 = T**2 + A * T**4
758
F4 = BETA/GAMMA * F3 / (1 + A*F3)
759
H = GAMMA * PHI**3 * LOG( 1 + F4 )
760
FC = ECUNIF + H
761
762
C Find correlation energy derivatives
763
DRSDD = - (THD * RS / DT)
764
DKFDD = THD * KF / DT
765
DKSDD = HALF * KS * DKFDD / KF
766
DZDD(1) = 1 / DT - ZETA / DT
767
DZDD(2) = - (1 / DT) - ZETA / DT
768
DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
769
DO 40 IS = 1,2
770
DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
771
DPDD = DPDZ * DZDD(IS)
772
DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
773
DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
774
DF2DD = (- F2) * DF1DD
775
DADD = (- A) * DF2DD / (F2-1)
776
DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
777
DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
778
DHDD = 3 * H * DPDD / PHI
779
DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
780
DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
781
782
DO 30 IX = 1,3
783
DTDGD = (T / GDMT) * GDT(IX) / GDMT
784
DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
785
DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )
786
DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
787
DFCDGD(IX,IS) = DT * DHDGD
788
30 CONTINUE
789
40 CONTINUE
790
791
C Find exchange energy and potential
792
FX = 0
793
DO 60 IS = 1,2
794
DS(IS) = MAX( DENMIN, 2 * D(IS) )
795
GDMS = MAX( GDMIN, 2 * GDM(IS) )
796
KFS = (3 * PI**2 * DS(IS))**THD
797
S = GDMS / (2 * KFS * DS(IS))
798
F1 = 1 + MU * S**2 / KAPPA
799
F = 1 + KAPPA - KAPPA / F1
800
c
801
c Note nspin=1 in call to exchng...
802
c
803
CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
804
FX = FX + DS(IS) * EXUNIF * F
805
806
DKFDD = THD * KFS / DS(IS)
807
DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
808
DF1DD = 2 * (F1-1) * DSDD / S
809
DFDD = KAPPA * DF1DD / F1**2
810
DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
811
812
DO 50 IX = 1,3
813
GDS = 2 * GD(IX,IS)
814
DSDGD = (S / GDMS) * GDS / GDMS
815
DF1DGD = 2 * MU * S * DSDGD / KAPPA
816
DFDGD = KAPPA * DF1DGD / F1**2
817
DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
818
50 CONTINUE
819
60 CONTINUE
820
FX = HALF * FX / DT
821
822
C Set output arguments
823
EX = FX
824
EC = FC
825
DO 90 IS = 1,nspin
826
DEXDD(IS) = DFXDD(IS)
827
DECDD(IS) = DFCDD(IS)
828
DO 80 IX = 1,3
829
DEXDGD(IX,IS) = DFXDGD(IX,IS)
830
DECDGD(IX,IS) = DFCDGD(IX,IS)
831
80 CONTINUE
832
90 CONTINUE
833
834
END
835
836
SUBROUTINE REVPBEXC( IREL, nspin, Dens, GDens,
837
. EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
838
839
C *********************************************************************
840
C Implements revPBE: revised Perdew-Burke-Ernzerhof GGA.
841
C Ref: Y. Zhang & W. Yang, Phys. Rev. Lett. 80, 890 (1998).
842
C Written by E. Artacho in January 2006 by modifying the PBE routine of
843
C L.C.Balbas and J.M.Soler. December 1996. Version 0.5.
844
C ******** INPUT ******************************************************
845
C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
846
C INTEGER nspin : Number of spin polarizations (1 or 2)
847
C REAL*8 Dens(nspin) : Total electron density (if nspin=1) or
848
C spin electron density (if nspin=2)
849
C REAL*8 GDens(3,nspin) : Total or spin density gradient
850
C ******** OUTPUT *****************************************************
851
C REAL*8 EX : Exchange energy density
852
C REAL*8 EC : Correlation energy density
853
C REAL*8 DEXDD(nspin) : Partial derivative
854
C d(DensTot*Ex)/dDens(ispin),
855
C where DensTot = Sum_ispin( Dens(ispin) )
856
C For a constant density, this is the
857
C exchange potential
858
C REAL*8 DECDD(nspin) : Partial derivative
859
C d(DensTot*Ec)/dDens(ispin),
860
C where DensTot = Sum_ispin( Dens(ispin) )
861
C For a constant density, this is the
862
C correlation potential
863
C REAL*8 DEXDGD(3,nspin): Partial derivative
864
C d(DensTot*Ex)/d(GradDens(i,ispin))
865
C REAL*8 DECDGD(3,nspin): Partial derivative
866
C d(DensTot*Ec)/d(GradDens(i,ispin))
867
C ********* UNITS ****************************************************
868
C Lengths in Bohr
869
C Densities in electrons per Bohr**3
870
C Energies in Hartrees
871
C Gradient vectors in cartesian coordinates
872
C ********* ROUTINES CALLED ******************************************
873
C EXCHNG, PW92C
874
C ********************************************************************
875
876
use precision, only : dp
877
878
implicit none
879
INTEGER IREL, nspin
880
real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
881
. DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
882
883
C Internal variables
884
INTEGER
885
. IS, IX
886
887
real(dp)
888
. A, BETA, D(2), DADD, DECUDD, DENMIN,
889
. DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
890
. DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
891
. DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
892
. DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
893
. EC, ECUNIF, EX, EXUNIF,
894
. F, F1, F2, F3, F4, FC, FX, FOUTHD,
895
. GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
896
. H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
897
. T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
898
899
C Lower bounds of density and its gradient to avoid divisions by zero
900
PARAMETER ( DENMIN = 1.D-12 )
901
PARAMETER ( GDMIN = 1.D-12 )
902
903
C Fix some numerical parameters
904
PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
905
. THD=1.D0/3.D0, THRHLF=1.5D0,
906
. TWO=2.D0, TWOTHD=2.D0/3.D0 )
907
908
C Fix some more numerical constants
909
PI = 4 * ATAN(1.D0)
910
BETA = 0.066725D0
911
GAMMA = (1 - LOG(TWO)) / PI**2
912
MU = BETA * PI**2 / 3
913
cea The only modification w.r.t. PBE in this following line.
914
KAPPA = 1.245D0
915
916
C Translate density and its gradient to new variables
917
IF (nspin .EQ. 1) THEN
918
D(1) = HALF * Dens(1)
919
D(2) = D(1)
920
DT = MAX( DENMIN, Dens(1) )
921
DO 10 IX = 1,3
922
GD(IX,1) = HALF * GDens(IX,1)
923
GD(IX,2) = GD(IX,1)
924
GDT(IX) = GDens(IX,1)
925
10 CONTINUE
926
ELSE
927
D(1) = Dens(1)
928
D(2) = Dens(2)
929
DT = MAX( DENMIN, Dens(1)+Dens(2) )
930
DO 20 IX = 1,3
931
GD(IX,1) = GDens(IX,1)
932
GD(IX,2) = GDens(IX,2)
933
GDT(IX) = GDens(IX,1) + GDens(IX,2)
934
20 CONTINUE
935
ENDIF
936
GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
937
GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
938
GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 )
939
GDMT = MAX( GDMIN, GDMT )
940
941
C Find local correlation energy and potential
942
CALL PW92C( 2, D, ECUNIF, VCUNIF )
943
944
C Find total correlation energy
945
RS = ( 3 / (4*PI*DT) )**THD
946
KF = (3 * PI**2 * DT)**THD
947
KS = SQRT( 4 * KF / PI )
948
ZETA = ( D(1) - D(2) ) / DT
949
ZETA = MAX( -1.D0+DENMIN, ZETA )
950
ZETA = MIN( 1.D0-DENMIN, ZETA )
951
PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
952
T = GDMT / (2 * PHI * KS * DT)
953
F1 = ECUNIF / GAMMA / PHI**3
954
F2 = EXP(-F1)
955
A = BETA / GAMMA / (F2-1)
956
F3 = T**2 + A * T**4
957
F4 = BETA/GAMMA * F3 / (1 + A*F3)
958
H = GAMMA * PHI**3 * LOG( 1 + F4 )
959
FC = ECUNIF + H
960
961
C Find correlation energy derivatives
962
DRSDD = - (THD * RS / DT)
963
DKFDD = THD * KF / DT
964
DKSDD = HALF * KS * DKFDD / KF
965
DZDD(1) = 1 / DT - ZETA / DT
966
DZDD(2) = - (1 / DT) - ZETA / DT
967
DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
968
DO 40 IS = 1,2
969
DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
970
DPDD = DPDZ * DZDD(IS)
971
DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
972
DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
973
DF2DD = (- F2) * DF1DD
974
DADD = (- A) * DF2DD / (F2-1)
975
DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
976
DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
977
DHDD = 3 * H * DPDD / PHI
978
DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
979
DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
980
981
DO 30 IX = 1,3
982
DTDGD = (T / GDMT) * GDT(IX) / GDMT
983
DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
984
DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )
985
DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
986
DFCDGD(IX,IS) = DT * DHDGD
987
30 CONTINUE
988
40 CONTINUE
989
990
C Find exchange energy and potential
991
FX = 0
992
DO 60 IS = 1,2
993
DS(IS) = MAX( DENMIN, 2 * D(IS) )
994
GDMS = MAX( GDMIN, 2 * GDM(IS) )
995
KFS = (3 * PI**2 * DS(IS))**THD
996
S = GDMS / (2 * KFS * DS(IS))
997
F1 = 1 + MU * S**2 / KAPPA
998
F = 1 + KAPPA - KAPPA / F1
999
c
1000
c Note nspin=1 in call to exchng...
1001
c
1002
CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
1003
FX = FX + DS(IS) * EXUNIF * F
1004
1005
DKFDD = THD * KFS / DS(IS)
1006
DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
1007
DF1DD = 2 * (F1-1) * DSDD / S
1008
DFDD = KAPPA * DF1DD / F1**2
1009
DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
1010
1011
DO 50 IX = 1,3
1012
GDS = 2 * GD(IX,IS)
1013
DSDGD = (S / GDMS) * GDS / GDMS
1014
DF1DGD = 2 * MU * S * DSDGD / KAPPA
1015
DFDGD = KAPPA * DF1DGD / F1**2
1016
DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
1017
50 CONTINUE
1018
60 CONTINUE
1019
FX = HALF * FX / DT
1020
1021
C Set output arguments
1022
EX = FX
1023
EC = FC
1024
DO 90 IS = 1,nspin
1025
DEXDD(IS) = DFXDD(IS)
1026
DECDD(IS) = DFCDD(IS)
1027
DO 80 IX = 1,3
1028
DEXDGD(IX,IS) = DFXDGD(IX,IS)
1029
DECDGD(IX,IS) = DFCDGD(IX,IS)
1030
80 CONTINUE
1031
90 CONTINUE
1032
1033
END
1034
1035
1036
SUBROUTINE PW91XC( IREL, nspin, Dens, GDens,
1037
. EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
1038
1039
C *********************************************************************
1040
C Implements Perdew-Wang91 Generalized-Gradient-Approximation.
1041
C Ref: JCP 100, 1290 (1994)
1042
C Written by J.L. Martins August 2000
1043
C ******** INPUT ******************************************************
1044
C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
1045
C INTEGER nspin : Number of spin polarizations (1 or 2)
1046
C REAL*8 Dens(nspin) : Total electron density (if nspin=1) or
1047
C spin electron density (if nspin=2)
1048
C REAL*8 GDens(3,nspin) : Total or spin density gradient
1049
C ******** OUTPUT *****************************************************
1050
C REAL*8 EX : Exchange energy density
1051
C REAL*8 EC : Correlation energy density
1052
C REAL*8 DEXDD(nspin) : Partial derivative
1053
C d(DensTot*Ex)/dDens(ispin),
1054
C where DensTot = Sum_ispin( Dens(ispin) )
1055
C For a constant density, this is the
1056
C exchange potential
1057
C REAL*8 DECDD(nspin) : Partial derivative
1058
C d(DensTot*Ec)/dDens(ispin),
1059
C where DensTot = Sum_ispin( Dens(ispin) )
1060
C For a constant density, this is the
1061
C correlation potential
1062
C REAL*8 DEXDGD(3,nspin): Partial derivative
1063
C d(DensTot*Ex)/d(GradDens(i,ispin))
1064
C REAL*8 DECDGD(3,nspin): Partial derivative
1065
C d(DensTot*Ec)/d(GradDens(i,ispin))
1066
C ********* UNITS ****************************************************
1067
C Lengths in Bohr
1068
C Densities in electrons per Bohr**3
1069
C Energies in Hartrees
1070
C Gradient vectors in cartesian coordinates
1071
C ********* ROUTINES CALLED ******************************************
1072
C EXCHNG, PW92C
1073
C ********************************************************************
1074
1075
use precision, only : dp
1076
1077
implicit none
1078
INTEGER IREL, nspin
1079
real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
1080
. DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
1081
1082
C Internal variables
1083
INTEGER
1084
. IS, IX
1085
real(dp)
1086
. A, BETA, D(2), DADD, DECUDD, DENMIN,
1087
. DF1DD, DF2DD, DF3DD, DF4DD, DF3DGD, DF4DGD,
1088
. DFCDD(2), DFCDGD(3,2), DFDGD, DFXDD(2), DFXDGD(3,2),
1089
. DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
1090
. DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
1091
. EC, ECUNIF, EX, EXUNIF,
1092
. F, F1, F2, F3, F4, FC, FX, FOUTHD,
1093
. GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
1094
. H, HALF, KF, KFS, KS, PHI, PI, RS, S,
1095
. T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
1096
1097
real(dp) F5, F6, F7, F8, ASINHS
1098
real(dp) DF5DD,DF6DD,DF7DD,DF8DD
1099
real(dp) DF1DS, DF2DS, DF3DS, DFDS, DF7DGD
1100
1101
C Lower bounds of density and its gradient to avoid divisions by zero
1102
PARAMETER ( DENMIN = 1.D-12 )
1103
PARAMETER ( GDMIN = 1.D-12 )
1104
1105
C Fix some numerical parameters
1106
PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
1107
. THD=1.D0/3.D0, THRHLF=1.5D0,
1108
. TWO=2.D0, TWOTHD=2.D0/3.D0 )
1109
1110
C Fix some more numerical constants
1111
PI = 4.0_dp * ATAN(1.0_dp)
1112
BETA = 15.75592_dp * 0.004235_dp
1113
GAMMA = BETA**2 / (2.0_dp * 0.09_dp)
1114
1115
C Translate density and its gradient to new variables
1116
IF (nspin .EQ. 1) THEN
1117
D(1) = HALF * Dens(1)
1118
D(2) = D(1)
1119
DT = MAX( DENMIN, Dens(1) )
1120
DO 10 IX = 1,3
1121
GD(IX,1) = HALF * GDens(IX,1)
1122
GD(IX,2) = GD(IX,1)
1123
GDT(IX) = GDens(IX,1)
1124
10 CONTINUE
1125
ELSE
1126
D(1) = Dens(1)
1127
D(2) = Dens(2)
1128
DT = MAX( DENMIN, Dens(1)+Dens(2) )
1129
DO 20 IX = 1,3
1130
GD(IX,1) = GDens(IX,1)
1131
GD(IX,2) = GDens(IX,2)
1132
GDT(IX) = GDens(IX,1) + GDens(IX,2)
1133
20 CONTINUE
1134
ENDIF
1135
GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
1136
GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
1137
GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 )
1138
GDMT = MAX( GDMIN, GDMT )
1139
1140
C Find local correlation energy and potential
1141
CALL PW92C( 2, D, ECUNIF, VCUNIF )
1142
1143
C Find total correlation energy
1144
RS = ( 3 / (4*PI*DT) )**THD
1145
KF = (3 * PI**2 * DT)**THD
1146
KS = SQRT( 4 * KF / PI )
1147
S = GDMT / (2 * KF * DT)
1148
T = GDMT / (2 * KS * DT)
1149
ZETA = ( D(1) - D(2) ) / DT
1150
ZETA = MAX( -1.D0+DENMIN, ZETA )
1151
ZETA = MIN( 1.D0-DENMIN, ZETA )
1152
PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
1153
F1 = ECUNIF / GAMMA / PHI**3
1154
F2 = EXP(-F1)
1155
A = BETA / GAMMA / (F2-1)
1156
F3 = T**2 + A * T**4
1157
F4 = BETA/GAMMA * F3 / (1 + A*F3)
1158
F5 = 0.002568D0 + 0.023266D0*RS + 7.389D-6*RS**2
1159
F6 = 1.0D0 + 8.723D0*RS + 0.472D0*RS**2 + 0.07389D0*RS**3
1160
F7 = EXP(-100.0D0 * S**2 * PHI**4)
1161
F8 = 15.75592D0*(0.001667212D0 + F5/F6 -0.004235D0 +
1162
. 3.0D0*0.001667212D0/7.0D0)
1163
H = GAMMA * PHI**3 * LOG( 1 + F4 ) + F8 * T**2 * F7
1164
FC = ECUNIF + H
1165
1166
C Find correlation energy derivatives
1167
DRSDD = - THD * RS / DT
1168
DKFDD = THD * KF / DT
1169
DKSDD = HALF * KS * DKFDD / KF
1170
DZDD(1) = 1 / DT - ZETA / DT
1171
DZDD(2) = - 1 / DT - ZETA / DT
1172
DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
1173
DO 40 IS = 1,2
1174
DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
1175
DPDD = DPDZ * DZDD(IS)
1176
DTDD = - T * ( DPDD/PHI + DKSDD/KS + 1/DT )
1177
DSDD = - S * ( DPDD/PHI + DKFDD/KF + 1/DT )
1178
DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
1179
DF2DD = - F2 * DF1DD
1180
DADD = - A * DF2DD / (F2-1)
1181
DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
1182
DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
1183
DF5DD = (0.023266D0 + 2.0D0*7.389D-6*RS)*DRSDD
1184
DF6DD = (8.723D0 + 2.0D0*0.472D0*RS
1185
. + 3.0D0*0.07389D0*RS**2)*DRSDD
1186
DF7DD = -200.0D0 * S * PHI**4 * DSDD * F7
1187
. -100.0D0 * S**2 * 4.0D0* PHI**3 * DPDD * F7
1188
DF8DD = 15.75592D0 * DF5DD/F6 - 15.75592D0*F5*DF6DD / F6**2
1189
DHDD = 3 * H * DPDD / PHI
1190
DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
1191
DHDD = DHDD + DF8DD * T**2 * F7
1192
DHDD = DHDD + F8 * 2*T*DTDD *F7
1193
DHDD = DHDD + F8 * T**2 * DF7DD
1194
1195
DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
1196
DO 30 IX = 1,3
1197
DTDGD = (T / GDMT) * GDT(IX) / GDMT
1198
DSDGD = (S / GDMT) * GDT(IX) / GDMT
1199
DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
1200
DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )
1201
DF7DGD = -200.0D0 * S * PHI**4 * DSDGD * F7
1202
DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
1203
DHDGD = DHDGD + F8 * 2*T*DTDGD *F7 + F8 * T**2 *DF7DGD
1204
DFCDGD(IX,IS) = DT * DHDGD
1205
30 CONTINUE
1206
40 CONTINUE
1207
1208
C Find exchange energy and potential
1209
FX = 0
1210
DO 60 IS = 1,2
1211
DS(1) = MAX( DENMIN, 2 * D(IS) )
1212
GDMS = MAX( GDMIN, 2 * GDM(IS) )
1213
KFS = (3 * PI**2 * DS(1))**THD
1214
S = GDMS / (2 * KFS * DS(1))
1215
F4 = SQRT(1.0D0 + (7.7956D0*S)**2)
1216
ASINHS = LOG(7.7956D0*S + F4)
1217
F1 = 1.0D0 + 0.19645D0 * S * ASINHS
1218
F2 = 0.2743D0 - 0.15084D0*EXP(-100.0D0*S*S)
1219
F3 = 1.0D0 / (F1 + 0.004D0 * S*S*S*S)
1220
F = (F1 + F2 * S*S ) * F3
1221
.
1222
CALL EXCHNG( IREL, 1, DS, EXUNIF, VXUNIF )
1223
FX = FX + DS(1) * EXUNIF * F
1224
1225
DKFDD = THD * KFS / DS(1)
1226
DSDD = S * ( -DKFDD/KFS - 1/DS(1) )
1227
DF1DS = 0.19645D0 * ASINHS +
1228
. 0.19645D0 * S * 7.7956D0 / F4
1229
DF2DS = 0.15084D0*200.0D0*S*EXP(-100.0D0*S*S)
1230
DF3DS = - F3*F3 * (DF1DS + 4.0D0*0.004D0 * S*S*S)
1231
DFDS = DF1DS * F3 + DF2DS * S*S * F3 + 2.0D0 * S * F2 * F3
1232
. + (F1 + F2 * S*S ) * DF3DS
1233
DFXDD(IS) = VXUNIF(1) * F + DS(1) * EXUNIF * DFDS * DSDD
1234
1235
DO 50 IX = 1,3
1236
GDS = 2 * GD(IX,IS)
1237
DSDGD = (S / GDMS) * GDS / GDMS
1238
DFDGD = DFDS * DSDGD
1239
DFXDGD(IX,IS) = DS(1) * EXUNIF * DFDGD
1240
50 CONTINUE
1241
60 CONTINUE
1242
FX = HALF * FX / DT
1243
1244
C Set output arguments
1245
EX = FX
1246
EC = FC
1247
DO 90 IS = 1,nspin
1248
DEXDD(IS) = DFXDD(IS)
1249
DECDD(IS) = DFCDD(IS)
1250
DO 80 IX = 1,3
1251
DEXDGD(IX,IS) = DFXDGD(IX,IS)
1252
DECDGD(IX,IS) = DFCDGD(IX,IS)
1253
80 CONTINUE
1254
90 CONTINUE
1255
1256
END
1257
1258
1259
1260
SUBROUTINE PW92C( nspin, Dens, EC, VC )
1261
1262
C ********************************************************************
1263
C Implements the Perdew-Wang'92 local correlation (beyond RPA).
1264
C Ref: J.P.Perdew & Y.Wang, PRB, 45, 13244 (1992)
1265
C Written by L.C.Balbas and J.M.Soler. Dec'96. Version 0.5.
1266
C ********* INPUT ****************************************************
1267
C INTEGER nspin : Number of spin polarizations (1 or 2)
1268
C REAL*8 Dens(nspin) : Local (spin) density
1269
C ********* OUTPUT ***************************************************
1270
C REAL*8 EC : Correlation energy density
1271
C REAL*8 VC(nspin) : Correlation (spin) potential
1272
C ********* UNITS ****************************************************
1273
C Densities in electrons per Bohr**3
1274
C Energies in Hartrees
1275
C ********* ROUTINES CALLED ******************************************
1276
C None
1277
C ********************************************************************
1278
1279
use precision, only : dp
1280
1281
C Next line is nonstandard but may be supressed
1282
implicit none
1283
1284
C Argument types and dimensions
1285
INTEGER nspin
1286
real(dp) Dens(nspin), EC, VC(nspin)
1287
1288
C Internal variable declarations
1289
INTEGER IG
1290
real(dp) A(0:2), ALPHA1(0:2), B, BETA(0:2,4), C,
1291
. DBDRS, DECDD(2), DECDRS, DECDZ, DENMIN, DFDZ,
1292
. DGDRS(0:2), DCDRS, DRSDD, DTOT, DZDD(2),
1293
. F, FPP0, FOUTHD, G(0:2), HALF, ONE,
1294
. P(0:2), PI, RS, THD, THRHLF, ZETA
1295
1296
C Add tiny numbers to avoid numerical errors
1297
PARAMETER ( DENMIN = 1.D-12 )
1298
PARAMETER ( ONE = 1.D0 + 1.D-12 )
1299
1300
C Fix some numerical constants
1301
PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
1302
. THD=1.D0/3.D0, THRHLF=1.5D0 )
1303
1304
C Parameters from Table I of Perdew & Wang, PRB, 45, 13244 (92)
1305
DATA P / 1.00d0, 1.00d0, 1.00d0 /
1306
DATA A / 0.031091d0, 0.015545d0, 0.016887d0 /
1307
DATA ALPHA1 / 0.21370d0, 0.20548d0, 0.11125d0 /
1308
DATA BETA / 7.5957d0, 14.1189d0, 10.357d0,
1309
. 3.5876d0, 6.1977d0, 3.6231d0,
1310
. 1.6382d0, 3.3662d0, 0.88026d0,
1311
. 0.49294d0, 0.62517d0, 0.49671d0 /
1312
1313
C Find rs and zeta
1314
PI = 4 * ATAN(1.D0)
1315
IF (nspin .EQ. 1) THEN
1316
DTOT = MAX( DENMIN, Dens(1) )
1317
ZETA = 0
1318
RS = ( 3 / (4*PI*DTOT) )**THD
1319
C Find derivatives dRs/dDens and dZeta/dDens
1320
DRSDD = (- RS) / DTOT / 3
1321
DZDD(1) = 0
1322
ELSE
1323
DTOT = MAX( DENMIN, Dens(1)+Dens(2) )
1324
ZETA = ( Dens(1) - Dens(2) ) / DTOT
1325
RS = ( 3 / (4*PI*DTOT) )**THD
1326
DRSDD = (- RS) / DTOT / 3
1327
DZDD(1) = (ONE - ZETA) / DTOT
1328
DZDD(2) = - (ONE + ZETA) / DTOT
1329
ENDIF
1330
1331
C Find eps_c(rs,0)=G(0), eps_c(rs,1)=G(1) and -alpha_c(rs)=G(2)
1332
C using eq.(10) of cited reference (Perdew & Wang, PRB, 45, 13244 (92))
1333
DO 20 IG = 0,2
1334
B = BETA(IG,1) * RS**HALF +
1335
. BETA(IG,2) * RS +
1336
. BETA(IG,3) * RS**THRHLF +
1337
. BETA(IG,4) * RS**(P(IG)+1)
1338
DBDRS = BETA(IG,1) * HALF / RS**HALF +
1339
. BETA(IG,2) +
1340
. BETA(IG,3) * THRHLF * RS**HALF +
1341
. BETA(IG,4) * (P(IG)+1) * RS**P(IG)
1342
C = 1 + 1 / (2 * A(IG) * B)
1343
DCDRS = - ( (C-1) * DBDRS / B )
1344
G(IG) = (- 2) * A(IG) * ( 1 + ALPHA1(IG)*RS ) * LOG(C)
1345
DGDRS(IG) = (- 2) *A(IG) * ( ALPHA1(IG) * LOG(C) +
1346
. (1+ALPHA1(IG)*RS) * DCDRS / C )
1347
20 CONTINUE
1348
1349
C Find f''(0) and f(zeta) from eq.(9)
1350
C = 1 / (2**FOUTHD - 2)
1351
FPP0 = 8 * C / 9
1352
F = ( (ONE+ZETA)**FOUTHD + (ONE-ZETA)**FOUTHD - 2 ) * C
1353
DFDZ = FOUTHD * ( (ONE+ZETA)**THD - (ONE-ZETA)**THD ) * C
1354
1355
C Find eps_c(rs,zeta) from eq.(8)
1356
EC = G(0) - G(2) * F / FPP0 * (ONE-ZETA**4) +
1357
. (G(1)-G(0)) * F * ZETA**4
1358
DECDRS = DGDRS(0) - DGDRS(2) * F / FPP0 * (ONE-ZETA**4) +
1359
. (DGDRS(1)-DGDRS(0)) * F * ZETA**4
1360
DECDZ = (- G(2)) / FPP0 * ( DFDZ*(ONE-ZETA**4) - F*4*ZETA**3 ) +
1361
. (G(1)-G(0)) * ( DFDZ*ZETA**4 + F*4*ZETA**3 )
1362
1363
C Find correlation potential
1364
IF (nspin .EQ. 1) THEN
1365
DECDD(1) = DECDRS * DRSDD
1366
VC(1) = EC + DTOT * DECDD(1)
1367
ELSE
1368
DECDD(1) = DECDRS * DRSDD + DECDZ * DZDD(1)
1369
DECDD(2) = DECDRS * DRSDD + DECDZ * DZDD(2)
1370
VC(1) = EC + DTOT * DECDD(1)
1371
VC(2) = EC + DTOT * DECDD(2)
1372
ENDIF
1373
1374
END
1375
1376
1377
1378
SUBROUTINE PW92XC( IREL, nspin, Dens, EPSX, EPSC, VX, VC )
1379
1380
C ********************************************************************
1381
C Implements the Perdew-Wang'92 LDA/LSD exchange correlation
1382
C Ref: J.P.Perdew & Y.Wang, PRB, 45, 13244 (1992)
1383
C Written by L.C.Balbas and J.M.Soler. Dec'96. Version 0.5.
1384
C ********* INPUT ****************************************************
1385
C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
1386
C INTEGER nspin : Number of spin polarizations (1 or 2)
1387
C REAL*8 Dens(nspin) : Local (spin) density
1388
C ********* OUTPUT ***************************************************
1389
C REAL*8 EPSX : Exchange energy density
1390
C REAL*8 EPSC : Correlation energy density
1391
C REAL*8 VX(nspin) : Exchange (spin) potential
1392
C REAL*8 VC(nspin) : Correlation (spin) potential
1393
C ********* UNITS ****************************************************
1394
C Densities in electrons per Bohr**3
1395
C Energies in Hartrees
1396
C ********* ROUTINES CALLED ******************************************
1397
C EXCHNG, PW92C
1398
C ********************************************************************
1399
1400
use precision, only : dp
1401
1402
implicit none
1403
INTEGER IREL, nspin
1404
real(dp) Dens(nspin), EPSX, EPSC, VC(nspin), VX(nspin)
1405
1406
CALL EXCHNG( IREL, nspin, Dens, EPSX, VX )
1407
CALL PW92C( nspin, Dens, EPSC, VC )
1408
END
1409
1410
1411
1412
SUBROUTINE PZXC( IREL, NSP, DS, EX, EC, VX, VC, DVXDN, DVCDN )
1413
1414
C *****************************************************************
1415
C Perdew-Zunger parameterization of Ceperley-Alder exchange and
1416
C correlation. Ref: Perdew & Zunger, Phys. Rev. B 23 5075 (1981).
1417
C Adapted by J.M.Soler from routine velect of Froyen's
1418
C pseudopotential generation program. Madrid, Jan'97.
1419
C **** Input *****************************************************
1420
C INTEGER IREL : relativistic-exchange switch (0=no, 1=yes)
1421
C INTEGER NSP : spin-polarizations (1=>unpolarized, 2=>polarized)
1422
C REAL*8 DS(NSP) : total (nsp=1) or spin (nsp=2) electron density
1423
C **** Output *****************************************************
1424
C REAL*8 EX : exchange energy density
1425
C REAL*8 EC : correlation energy density
1426
C REAL*8 VX(NSP) : (spin-dependent) exchange potential
1427
C REAL*8 VC(NSP) : (spin-dependent) correlation potential
1428
C REAL*8 DVXDN(NSP,NSP): Derivative of the exchange potential
1429
C respect the charge density,
1430
C Dvx(spin1)/Dn(spin2)
1431
C REAL*8 DVCDN(NSP,NSP): Derivative of the correlation potential
1432
C respect the charge density,
1433
C Dvc(spin1)/Dn(spin2)
1434
C **** Units *******************************************************
1435
C Densities in electrons/Bohr**3
1436
C Energies in Hartrees
1437
C *****************************************************************
1438
1439
use precision, only: dp
1440
1441
implicit none
1442
1443
integer :: nsp, irel, isp1, isp2, isp
1444
real(dp) :: DS(NSP), VX(NSP), VC(NSP),
1445
. DVXDN(NSP,NSP), DVCDN(NSP,NSP)
1446
real(dp), parameter ::
1447
$ ZERO=0.D0,ONE=1.D0,PFIVE=.5D0,OPF=1.5D0,PNN=.99D0,
1448
$ PTHREE=0.3D0,PSEVF=0.75D0,C0504=0.0504D0,
1449
$ C0254=0.0254D0,C014=0.014D0,C0406=0.0406D0,
1450
$ C15P9=15.9D0,C0666=0.0666D0,C11P4=11.4D0,
1451
$ C045=0.045D0,C7P8=7.8D0,C88=0.88D0,C20P59=20.592D0,
1452
$ C3P52=3.52D0,C0311=0.0311D0,C0014=0.0014D0,
1453
$ C0538=0.0538D0,C0096=0.0096D0,C096=0.096D0,
1454
$ C0622=0.0622D0,C004=0.004D0,C0232=0.0232D0,
1455
$ C1686=0.1686D0,C1P398=1.3981D0,C2611=0.2611D0,
1456
$ C2846=0.2846D0,C1P053=1.0529D0,C3334=0.3334D0
1457
1458
C Ceperly-Alder 'ca' constants. Internal energies in Rydbergs.
1459
real(dp), parameter ::
1460
$ CON1=1.D0/6, CON2=0.008D0/3, CON3=0.3502D0/3,
1461
$ CON4=0.0504D0/3, CON5=0.0028D0/3, CON6=0.1925D0/3,
1462
$ CON7=0.0206D0/3, CON8=9.7867D0/6, CON9=1.0444D0/3,
1463
$ CON10=7.3703D0/6, CON11=1.3336D0/3
1464
1465
C X-alpha parameter:
1466
real(dp), PARAMETER :: ALP = 2.D0 / 3.D0
1467
1468
C Other variables converted into parameters by J.M.Soler
1469
real(dp), parameter ::
1470
$ TINY = 1.D-6 ,
1471
$ PI = 3.14159265358979312_dp,
1472
$ TWO = 2.0D0,
1473
$ HALF = 0.5D0,
1474
$ TRD = 1.D0 / 3.D0,
1475
$ FTRD = 4.D0 / 3.D0,
1476
$ TFTM = 0.51984209978974638D0,
1477
$ A0 = 0.52106176119784808D0,
1478
$ CRS = 0.620350490899400087D0,
1479
$ CXP = (- 3.D0) * ALP / (PI*A0),
1480
$ CXF = 1.25992104989487319D0
1481
1482
real(dp) :: d1, d2, d, z, fz, fzp
1483
real(dp) :: ex, ec, dfzpdn, rs, vxp, exp_var
1484
real(dp) :: beta, sb, alb, vxf, exf, dvxpdn
1485
real(dp) :: dvxfdn, sqrs, te, be, ecp, vcp
1486
real(dp) :: dtedn, be2, dbedn, dvcpdn, decpdn
1487
real(dp) :: ecf, vcf, dvcfdn, decfdn, rslog
1488
1489
1490
C Find density and polarization
1491
IF (NSP .EQ. 2) THEN
1492
D1 = MAX(DS(1),ZERO)
1493
D2 = MAX(DS(2),ZERO)
1494
D = D1 + D2
1495
IF (D .LE. ZERO) THEN
1496
EX = ZERO
1497
EC = ZERO
1498
VX(1) = ZERO
1499
VX(2) = ZERO
1500
VC(1) = ZERO
1501
VC(2) = ZERO
1502
RETURN
1503
ENDIF
1504
c
1505
c Robustness enhancement by Jose Soler (August 2002)
1506
c
1507
Z = (D1 - D2) / D
1508
IF (Z .LE. -ONE) THEN
1509
FZ = (TWO**FTRD-TWO)/TFTM
1510
FZP = -FTRD*TWO**TRD/TFTM
1511
DFZPDN = FTRD*TRD*TWO**(-ALP)/TFTM
1512
ELSEIF (Z .GE. ONE) THEN
1513
FZ = (TWO**FTRD-TWO)/TFTM
1514
FZP = FTRD*TWO**TRD/TFTM
1515
DFZPDN = FTRD*TRD*TWO**(-ALP)/TFTM
1516
ELSE
1517
FZ = ((ONE+Z)**FTRD+(ONE-Z)**FTRD-TWO)/TFTM
1518
FZP = FTRD*((ONE+Z)**TRD-(ONE-Z)**TRD)/TFTM
1519
DFZPDN = FTRD*TRD*((ONE+Z)**(-ALP) + (ONE-Z)**(-ALP))/TFTM
1520
ENDIF
1521
ELSE
1522
D = DS(1)
1523
IF (D .LE. ZERO) THEN
1524
EX = ZERO
1525
EC = ZERO
1526
VX(1) = ZERO
1527
VC(1) = ZERO
1528
RETURN
1529
ENDIF
1530
Z = ZERO
1531
FZ = ZERO
1532
FZP = ZERO
1533
ENDIF
1534
RS = CRS / D**TRD
1535
1536
C Exchange
1537
VXP = CXP / RS
1538
EXP_VAR = 0.75D0 * VXP
1539
IF (IREL .EQ. 1) THEN
1540
BETA = C014/RS
1541
IF (BETA .LT. TINY) THEN
1542
SB = ONE + HALF*BETA**2
1543
ALB = BETA
1544
ELSE
1545
SB = SQRT(1+BETA*BETA)
1546
ALB = LOG(BETA+SB)
1547
ENDIF
1548
VXP = VXP * (-PFIVE + OPF * ALB / (BETA*SB))
1549
EXP_VAR = EXP_VAR *(ONE-OPF*((BETA*SB-ALB)/BETA**2)**2)
1550
ENDIF
1551
VXF = CXF * VXP
1552
EXF = CXF * EXP_VAR
1553
DVXPDN = TRD * VXP / D
1554
DVXFDN = TRD * VXF / D
1555
1556
C Correlation
1557
IF (RS .GT. ONE) THEN
1558
SQRS=SQRT(RS)
1559
TE = ONE+CON10*SQRS+CON11*RS
1560
BE = ONE+C1P053*SQRS+C3334*RS
1561
ECP = -(C2846/BE)
1562
VCP = ECP*TE/BE
1563
DTEDN = ((CON10 * SQRS *HALF) + CON11 * RS)*(-TRD/D)
1564
BE2 = BE * BE
1565
DBEDN = ((C1P053 * SQRS *HALF) + C3334 * RS)*(-TRD/D)
1566
DVCPDN = -(C2846/BE2)*(DTEDN - 2.0D0 * TE * DBEDN/BE)
1567
DECPDN = (C2846/BE2)*DBEDN
1568
TE = ONE+CON8*SQRS+CON9*RS
1569
BE = ONE+C1P398*SQRS+C2611*RS
1570
ECF = -(C1686/BE)
1571
VCF = ECF*TE/BE
1572
DTEDN = ((CON8 * SQRS * HALF) + CON9 * RS)*(-TRD/D)
1573
BE2 = BE * BE
1574
DBEDN = ((C1P398 * SQRS * HALF) + C2611 * RS)*(-TRD/D)
1575
DVCFDN = -(C1686/BE2)*(DTEDN - 2.0D0 * TE * DBEDN/BE)
1576
DECFDN = (C1686/BE2)*DBEDN
1577
ELSE
1578
RSLOG=LOG(RS)
1579
ECP=(C0622+C004*RS)*RSLOG-C096-C0232*RS
1580
VCP=(C0622+CON2*RS)*RSLOG-CON3-CON4*RS
1581
DVCPDN = (CON2*RS*RSLOG + (CON2-CON4)*RS + C0622)*(-TRD/D)
1582
DECPDN = (C004*RS*RSLOG + (C004-C0232)*RS + C0622)*(-TRD/D)
1583
ECF=(C0311+C0014*RS)*RSLOG-C0538-C0096*RS
1584
VCF=(C0311+CON5*RS)*RSLOG-CON6-CON7*RS
1585
DVCFDN = (CON5*RS*RSLOG + (CON5-CON7)*RS + C0311)*(-TRD/D)
1586
DECFDN = (C0014*RS*RSLOG + (C0014-C0096)*RS + C0311)*(-TRD/D)
1587
ENDIF
1588
1589
ISP1 = 1
1590
ISP2 = 2
1591
1592
C Find up and down potentials
1593
IF (NSP .EQ. 2) THEN
1594
EX = EXP_VAR + FZ*(EXF-EXP_VAR)
1595
EC = ECP + FZ*(ECF-ECP)
1596
VX(1) = VXP + FZ*(VXF-VXP) + (ONE-Z)*FZP*(EXF-EXP_VAR)
1597
VX(2) = VXP + FZ*(VXF-VXP) - (ONE+Z)*FZP*(EXF-EXP_VAR)
1598
VC(1) = VCP + FZ*(VCF-VCP) + (ONE-Z)*FZP*(ECF-ECP)
1599
VC(2) = VCP + FZ*(VCF-VCP) - (ONE+Z)*FZP*(ECF-ECP)
1600
1601
C Derivatives of exchange potential respect the density
1602
1603
DVXDN(ISP1,ISP1) =
1604
. DVXPDN
1605
. + FZP*(VXF-VXP-EXF+EXP_VAR)*( 2.D0*D2/(D*D) )
1606
. + FZ*(DVXFDN-DVXPDN)+(1-Z)*FZP*(VXF-VXP)/(4.D0*D)
1607
. + (1-Z)*DFZPDN*(EXF-EXP_VAR)*( 2.D0*D2/(D*D) )
1608
DVXDN(ISP1,ISP2) =
1609
. DVXPDN
1610
. + FZP*(VXF-VXP-EXF+EXP_VAR)*(-2.D0*D1/(D*D) )
1611
. + FZ*(DVXFDN-DVXPDN)+(1-Z)*FZP*(VXF-VXP)/(4.D0*D)
1612
. + (1-Z)*DFZPDN*(EXF-EXP_VAR)*( -2.D0*D1/(D*D) )
1613
DVXDN(ISP2,ISP1) =
1614
. DVXPDN
1615
. + FZP*(VXF-VXP-EXF+EXP_VAR)*( 2.D0*D2/(D*D) )
1616
. + FZ*(DVXFDN-DVXPDN)-(1+Z)*FZP*(VXF-VXP)/(4.D0*D)
1617
. - (1+Z)*DFZPDN*(EXF-EXP_VAR)*( 2.D0*D2/(D*D) )
1618
DVXDN(ISP2,ISP2) =
1619
. DVXPDN
1620
. + FZP*(VXF-VXP-EXF+EXP_VAR)*(-2.D0*D1/(D*D) )
1621
. + FZ*(DVXFDN-DVXPDN)-(1+Z)*FZP*(VXF-VXP)/(4.D0*D)
1622
. - (1+Z)*DFZPDN*(EXF-EXP_VAR)*( -2.D0*D1/(D*D) )
1623
1624
C Derivatives of correlation potential respect the density
1625
1626
DVCDN(ISP1,ISP1) =
1627
. DVCPDN
1628
. + FZP*(VCF-VCP-ECF+ECP)*( 2.D0*D2/(D*D) )
1629
. + FZ*(DVCFDN-DVCPDN)+ (1-Z)*FZP*(DECFDN-DECPDN)
1630
. + (1-Z)*DFZPDN*(ECF-ECP)*( 2.D0*D2/(D*D) )
1631
DVCDN(ISP1,ISP2) =
1632
. DVCPDN
1633
. + FZP*(VCF-VCP-ECF+ECP)*(-2.D0*D1/(D*D) )
1634
. + FZ*(DVCFDN-DVCPDN)+ (1-Z)*FZP*(DECFDN-DECPDN)
1635
. + (1-Z)*DFZPDN*(ECF-ECP)*( -2.D0*D1/(D*D) )
1636
DVCDN(ISP2,ISP1) =
1637
. DVCPDN
1638
. + FZP*(VCF-VCP-ECF+ECP)*( 2.D0*D2/(D*D) )
1639
. + FZ*(DVCFDN-DVCPDN)- (1+Z)*FZP*(DECFDN-DECPDN)
1640
. - (1+Z)*DFZPDN*(ECF-ECP)*( 2.D0*D2/(D*D) )
1641
DVCDN(ISP2,ISP2) =
1642
. DVCPDN
1643
. + FZP*(VCF-VCP-ECF+ECP)*(-2.D0*D1/(D*D) )
1644
. + FZ*(DVCFDN-DVCPDN)- (1+Z)*FZP*(DECFDN-DECPDN)
1645
. - (1+Z)*DFZPDN*(ECF-ECP)*( -2.D0*D1/(D*D) )
1646
1647
ELSE
1648
EX = EXP_VAR
1649
EC = ECP
1650
VX(1) = VXP
1651
VC(1) = VCP
1652
DVXDN(1,1) = DVXPDN
1653
DVCDN(1,1) = DVCPDN
1654
ENDIF
1655
1656
C Change from Rydbergs to Hartrees
1657
EX = HALF * EX
1658
EC = HALF * EC
1659
DO 10 ISP = 1,NSP
1660
VX(ISP) = HALF * VX(ISP)
1661
VC(ISP) = HALF * VC(ISP)
1662
DO 5 ISP2 = 1,NSP
1663
DVXDN(ISP,ISP2) = HALF * DVXDN(ISP,ISP2)
1664
DVCDN(ISP,ISP2) = HALF * DVCDN(ISP,ISP2)
1665
5 CONTINUE
1666
10 CONTINUE
1667
END
1668
1669
subroutine blypxc(nspin,dens,gdens,EX,EC,
1670
. dEXdd,dECdd,dEXdgd,dECdgd)
1671
c ***************************************************************
1672
c Implements Becke gradient exchange functional (A.D.
1673
c Becke, Phys. Rev. A 38, 3098 (1988)) and Lee, Yang, Parr
1674
c correlation functional (C. Lee, W. Yang, R.G. Parr, Phys. Rev. B
1675
c 37, 785 (1988)), as modificated by Miehlich,Savin,Stoll and Preuss,
1676
c Chem. Phys. Lett. 157,200 (1989). See also Johnson, Gill and Pople,
1677
c J. Chem. Phys. 98, 5612 (1993). Some errors were detected in this
1678
c last paper, so not all of the expressions correspond exactly to those
1679
c implemented here.
1680
c Written by Maider Machado. July 1998.
1681
c **************** INPUT ********************************************
1682
c integer nspin : Number of spin polarizations (1 or 2)
1683
c real*8 dens(nspin) : Total electron density (if nspin=1) or
1684
c spin electron density (if nspin=2)
1685
c real*8 gdens(3,nspin) : Total or spin density gradient
1686
c ******** OUTPUT *****************************************************
1687
c real*8 ex : Exchange energy density
1688
c real*8 ec : Correlation energy density
1689
c real*8 dexdd(nspin) : Partial derivative
1690
c d(DensTot*Ex)/dDens(ispin),
1691
c where DensTot = Sum_ispin( Dens(ispin) )
1692
c For a constant density, this is the
1693
c exchange potential
1694
c real*8 decdd(nspin) : Partial derivative
1695
c d(DensTot*Ec)/dDens(ispin),
1696
c where DensTot = Sum_ispin( Dens(ispin) )
1697
c For a constant density, this is the
1698
c correlation potential
1699
c real*8 dexdgd(3,nspin): Partial derivative
1700
c d(DensTot*Ex)/d(GradDens(i,ispin))
1701
c real*8 decdgd(3,nspin): Partial derivative
1702
c d(DensTot*Ec)/d(GradDens(i,ispin))
1703
c ********* UNITS ****************************************************
1704
c Lengths in Bohr
1705
c Densities in electrons per Bohr**3
1706
c Energies in Hartrees
1707
c Gradient vectors in cartesian coordinates
1708
c ********************************************************************
1709
1710
use precision, only : dp
1711
1712
implicit none
1713
1714
integer nspin
1715
real(dp) dens(nspin), gdens(3,nspin), EX, EC,
1716
. dEXdd(nspin), dECdd(nspin), dEXdgd(3,nspin),
1717
. dECdgd(3,nspin)
1718
1719
c Internal variables
1720
integer is,ix
1721
real(dp) pi, beta, thd, tthd, thrhlf, half, fothd,
1722
. d(2),gd(3,2),dmin, ash,gdm(2),denmin,dt,
1723
. g(2),x(2),a,b,c,dd,onzthd,gdmin,
1724
. ga, gb, gc,becke,dbecgd(3,2),
1725
. dgdx(2), dgdxa, dgdxb, dgdxc,dgdxd,dbecdd(2),
1726
. den,omega, domega, delta, ddelta,cf,
1727
. gam11, gam12, gam22, LYPa, LYPb1,
1728
. LYPb2,dLYP11,dLYP12,dLYP22,LYP,
1729
. dd1g11,dd1g12,dd1g22,dd2g12,dd2g11,dd2g22,
1730
. dLYPdd(2),dg11dd(3,2),dg22dd(3,2),
1731
. dLYPgd(3,2)
1732
1733
c Lower bounds of density and its gradient to avoid divisions by zero
1734
parameter ( denmin=1.d-8 )
1735
parameter (gdmin=1.d-8)
1736
parameter (dmin=1.d-5)
1737
1738
c Fix some numerical parameters
1739
parameter ( thd = 1.d0/3.d0, tthd=2.d0/3.d0 )
1740
parameter ( thrhlf=1.5d0, half=0.5d0,
1741
. fothd=4.d0/3.d0, onzthd=11.d0/3.d0)
1742
1743
c Empirical parameter for Becke exchange functional (a.u.)
1744
parameter(beta= 0.0042d0)
1745
1746
c Constants for LYP functional (a.u.)
1747
parameter(a=0.04918d0, b=0.132d0, c=0.2533d0, dd=0.349d0)
1748
1749
pi= 4*atan(1.d0)
1750
1751
1752
c Translate density and its gradient to new variables
1753
if (nspin .eq. 1) then
1754
d(1) = half * dens(1)
1755
d(1) = max(denmin,d(1))
1756
d(2) = d(1)
1757
dt = max( denmin, dens(1) )
1758
do ix = 1,3
1759
gd(ix,1) = half * gdens(ix,1)
1760
gd(ix,2) = gd(ix,1)
1761
enddo
1762
else
1763
d(1) = dens(1)
1764
d(2) = dens(2)
1765
do is=1,2
1766
d(is) = max (denmin,d(is))
1767
enddo
1768
dt = max( denmin, dens(1)+dens(2) )
1769
do ix = 1,3
1770
gd(ix,1) = gdens(ix,1)
1771
gd(ix,2) = gdens(ix,2)
1772
enddo
1773
endif
1774
1775
gdm(1) = sqrt( gd(1,1)**2 + gd(2,1)**2 + gd(3,1)**2 )
1776
gdm(2) = sqrt( gd(1,2)**2 + gd(2,2)**2 + gd(3,2)**2 )
1777
1778
do is=1,2
1779
gdm(is)= max(gdm(is),gdmin)
1780
enddo
1781
1782
c Find Becke exchange energy
1783
ga = -thrhlf*(3.d0/4.d0/pi)**thd
1784
do is=1,2
1785
if(d(is).lt.dmin) then
1786
g(is)=ga
1787
else
1788
x(is) = gdm(is)/d(is)**fothd
1789
gb = beta*x(is)**2
1790
ash=log(x(is)+sqrt(x(is)**2+1))
1791
gc = 1+6*beta*x(is)*ash
1792
g(is) = ga-gb/gc
1793
endif
1794
enddo
1795
1796
c Density of energy
1797
becke=(g(1)*d(1)**fothd+g(2)*d(2)**fothd)/dt
1798
1799
1800
c Exchange energy derivatives
1801
do is=1,2
1802
if(d(is).lt.dmin)then
1803
dbecdd(is)=0.
1804
do ix=1,3
1805
dbecgd(ix,is)=0.
1806
enddo
1807
else
1808
dgdxa=6*beta**2*x(is)**2
1809
ash=log(x(is)+sqrt(x(is)**2+1))
1810
dgdxb=x(is)/sqrt(x(is)**2+1)-ash
1811
dgdxc=-2*beta*x(is)
1812
dgdxd=(1+6*beta*x(is)*ash)**2
1813
dgdx(is)=(dgdxa*dgdxb+dgdxc)/dgdxd
1814
dbecdd(is)=fothd*d(is)**thd*(g(is)-x(is)*dgdx(is))
1815
do ix=1,3
1816
dbecgd(ix,is)=d(is)**(-fothd)*dgdx(is)*gd(ix,is)/x(is)
1817
enddo
1818
endif
1819
enddo
1820
1821
c Lee-Yang-Parr correlation energy
1822
den=1+dd*dt**(-thd)
1823
omega=dt**(-onzthd)*exp(-c*dt**(-thd))/den
1824
delta=c*dt**(-thd)+dd*dt**(-thd)/den
1825
cf=3.*(3*pi**2)**tthd/10.
1826
gam11=gdm(1)**2
1827
gam12=gd(1,1)*gd(1,2)+gd(2,1)*gd(2,2)+gd(3,1)*gd(3,2)
1828
gam22=gdm(2)**2
1829
LYPa=-4*a*d(1)*d(2)/(den*dt)
1830
LYPb1=2**onzthd*cf*a*b*omega*d(1)*d(2)
1831
LYPb2=d(1)**(8./3.)+d(2)**(8./3.)
1832
dLYP11=-a*b*omega*(d(1)*d(2)/9.*(1.-3.*delta-(delta-11.)
1833
.*d(1)/dt)-d(2)**2)
1834
dLYP12=-a*b*omega*(d(1)*d(2)/9.*(47.-7.*delta)
1835
.-fothd*dt**2)
1836
dLYP22=-a*b*omega*(d(1)*d(2)/9.*(1.-3.*delta-(delta-11.)*
1837
.d(2)/dt)-d(1)**2)
1838
1839
c Density of energy
1840
LYP=(LYPa-LYPb1*LYPb2+dLYP11*gam11+dLYP12*gam12
1841
.+dLYP22*gam22)/dt
1842
1843
c Correlation energy derivatives
1844
domega=-thd*dt**(-fothd)*omega*(11.*dt**thd-c-dd/den)
1845
ddelta=thd*(dd**2*dt**(-5./3.)/den**2-delta/dt)
1846
1847
c Second derivatives with respect to the density
1848
dd1g11=domega/omega*dLYP11-a*b*omega*(d(2)/9.*
1849
. (1.-3.*delta-2*(delta-11.)*d(1)/dt)-d(1)*d(2)/9.*
1850
. ((3.+d(1)/dt)*ddelta-(delta-11.)*d(1)/dt**2))
1851
1852
dd1g12=domega/omega*dLYP12-a*b*omega*(d(2)/9.*
1853
. (47.-7.*delta)-7./9.*d(1)*d(2)*ddelta-8./3.*dt)
1854
1855
dd1g22=domega/omega*dLYP22-a*b*omega*(1./9.*d(2)
1856
. *(1.-3.*delta-(delta-11.)*d(2)/dt)-d(1)*d(2)/9.*
1857
. ((3.+d(2)/dt)*ddelta-(delta-11.)*d(2)/dt**2)-2*d(1))
1858
1859
1860
dd2g22=domega/omega*dLYP22-a*b*omega*(d(1)/9.*
1861
. (1.-3.*delta-2*(delta-11.)*d(2)/dt)-d(1)*d(2)/9.*
1862
. ((3+d(2)/dt)*ddelta-(delta-11.)*d(2)/dt**2))
1863
1864
1865
dd2g12=domega/omega*dLYP12-a*b*omega*(d(1)/9.*
1866
. (47.-7.*delta)-7./9.*d(1)*d(2)*ddelta-8./3.*dt)
1867
1868
dd2g11=domega/omega*dLYP11-a*b*omega*(1./9.*d(1)
1869
. *(1.-3.*delta-(delta-11.)*d(1)/dt)-d(1)*d(2)/9.*
1870
. ((3.+d(1)/dt)*ddelta-(delta-11.)*d(1)/dt**2)-2*d(2))
1871
1872
1873
dLYPdd(1)=-4*a/den*d(1)*d(2)/dt*
1874
. (thd*dd*dt**(-fothd)/den
1875
. +1./d(1)-1./dt)-2**onzthd*cf*a*b*(domega*d(1)*d(2)*
1876
. (d(1)**(8./3.)+d(2)**(8./3.))+omega*d(2)*(onzthd*
1877
. d(1)**(8./3.)+d(2)**(8./3.)))+dd1g11*gam11+
1878
. dd1g12*gam12+dd1g22*gam22
1879
1880
1881
dLYPdd(2)=-4*a/den*d(1)*d(2)/dt*(thd*dd*dt**(-fothd)/den
1882
. +1./d(2)-1./dt)-2**onzthd*cf*a*b*(domega*d(1)*d(2)*
1883
. (d(1)**(8./3.)+d(2)**(8./3.))+omega*d(1)*(onzthd*
1884
. d(2)**(8./3.)+d(1)**(8./3.)))+dd2g22*gam22+
1885
. dd2g12*gam12+dd2g11*gam11
1886
1887
1888
c second derivatives with respect to the density gradient
1889
1890
do is=1,2
1891
do ix=1,3
1892
dg11dd(ix,is)=2*gd(ix,is)
1893
dg22dd(ix,is)=2*gd(ix,is)
1894
enddo
1895
enddo
1896
do ix=1,3
1897
dLYPgd(ix,1)=dLYP11*dg11dd(ix,1)+dLYP12*gd(ix,2)
1898
dLYPgd(ix,2)=dLYP22*dg22dd(ix,2)+dLYP12*gd(ix,1)
1899
enddo
1900
1901
1902
EX=becke
1903
EC=LYP
1904
do is=1,nspin
1905
dEXdd(is)=dbecdd(is)
1906
dECdd(is)=dLYPdd(is)
1907
do ix=1,3
1908
dEXdgd(ix,is)=dbecgd(ix,is)
1909
dECdgd(ix,is)=dLYPgd(ix,is)
1910
enddo
1911
enddo
1912
end
1913
1914
SUBROUTINE RPBEXC( IREL, nspin, Dens, GDens,
1915
. EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
1916
1917
C *********************************************************************
1918
C Implements Hammer's RPBE Generalized-Gradient-Approximation (GGA).
1919
C A revision of PBE (Perdew-Burke-Ernzerhof)
1920
C Ref: Hammer, Hansen & Norskov, PRB 59, 7413 (1999) and
1921
C J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996)
1922
C
1923
C Written by M.V. Fernandez-Serra. March 2004. On the PBE routine of
1924
C L.C.Balbas and J.M.Soler. December 1996. Version 0.5.
1925
C ******** INPUT ******************************************************
1926
C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
1927
C INTEGER nspin : Number of spin polarizations (1 or 2)
1928
C REAL*8 Dens(nspin) : Total electron density (if nspin=1) or
1929
C spin electron density (if nspin=2)
1930
C REAL*8 GDens(3,nspin) : Total or spin density gradient
1931
C ******** OUTPUT *****************************************************
1932
C REAL*8 EX : Exchange energy density
1933
C REAL*8 EC : Correlation energy density
1934
C REAL*8 DEXDD(nspin) : Partial derivative
1935
C d(DensTot*Ex)/dDens(ispin),
1936
C where DensTot = Sum_ispin( Dens(ispin) )
1937
C For a constant density, this is the
1938
C exchange potential
1939
C REAL*8 DECDD(nspin) : Partial derivative
1940
C d(DensTot*Ec)/dDens(ispin),
1941
C where DensTot = Sum_ispin( Dens(ispin) )
1942
C For a constant density, this is the
1943
C correlation potential
1944
C REAL*8 DEXDGD(3,nspin): Partial derivative
1945
C d(DensTot*Ex)/d(GradDens(i,ispin))
1946
C REAL*8 DECDGD(3,nspin): Partial derivative
1947
C d(DensTot*Ec)/d(GradDens(i,ispin))
1948
C ********* UNITS ****************************************************
1949
C Lengths in Bohr
1950
C Densities in electrons per Bohr**3
1951
C Energies in Hartrees
1952
C Gradient vectors in cartesian coordinates
1953
C ********* ROUTINES CALLED ******************************************
1954
C EXCHNG, PW92C
1955
C ********************************************************************
1956
1957
use precision, only : dp
1958
1959
implicit none
1960
INTEGER IREL, nspin
1961
real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
1962
. DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
1963
1964
C Internal variables
1965
INTEGER
1966
. IS, IX
1967
1968
real(dp)
1969
. A, BETA, D(2), DADD, DECUDD, DENMIN,
1970
. DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
1971
. DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
1972
. DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
1973
. DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
1974
. EC, ECUNIF, EX, EXUNIF,
1975
. F, F1, F2, F3, F4, FC, FX, FOUTHD,
1976
. GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
1977
. H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
1978
. T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
1979
1980
C Lower bounds of density and its gradient to avoid divisions by zero
1981
PARAMETER ( DENMIN = 1.D-12 )
1982
PARAMETER ( GDMIN = 1.D-12 )
1983
1984
C Fix some numerical parameters
1985
PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
1986
. THD=1.D0/3.D0, THRHLF=1.5D0,
1987
. TWO=2.D0, TWOTHD=2.D0/3.D0 )
1988
1989
C Fix some more numerical constants
1990
PI = 4 * ATAN(1.D0)
1991
BETA = 0.066725D0
1992
GAMMA = (1 - LOG(TWO)) / PI**2
1993
MU = BETA * PI**2 / 3
1994
KAPPA = 0.804D0
1995
1996
C Translate density and its gradient to new variables
1997
IF (nspin .EQ. 1) THEN
1998
D(1) = HALF * Dens(1)
1999
D(2) = D(1)
2000
DT = MAX( DENMIN, Dens(1) )
2001
DO 10 IX = 1,3
2002
GD(IX,1) = HALF * GDens(IX,1)
2003
GD(IX,2) = GD(IX,1)
2004
GDT(IX) = GDens(IX,1)
2005
10 CONTINUE
2006
ELSE
2007
D(1) = Dens(1)
2008
D(2) = Dens(2)
2009
DT = MAX( DENMIN, Dens(1)+Dens(2) )
2010
DO 20 IX = 1,3
2011
GD(IX,1) = GDens(IX,1)
2012
GD(IX,2) = GDens(IX,2)
2013
GDT(IX) = GDens(IX,1) + GDens(IX,2)
2014
20 CONTINUE
2015
ENDIF
2016
GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
2017
GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
2018
GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 )
2019
GDMT = MAX( GDMIN, GDMT )
2020
2021
C Find local correlation energy and potential
2022
CALL PW92C( 2, D, ECUNIF, VCUNIF )
2023
2024
C Find total correlation energy
2025
RS = ( 3 / (4*PI*DT) )**THD
2026
KF = (3 * PI**2 * DT)**THD
2027
KS = SQRT( 4 * KF / PI )
2028
ZETA = ( D(1) - D(2) ) / DT
2029
ZETA = MAX( -1.D0+DENMIN, ZETA )
2030
ZETA = MIN( 1.D0-DENMIN, ZETA )
2031
PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
2032
T = GDMT / (2 * PHI * KS * DT)
2033
F1 = ECUNIF / GAMMA / PHI**3
2034
F2 = EXP(-F1)
2035
A = BETA / GAMMA / (F2-1)
2036
F3 = T**2 + A * T**4
2037
F4 = BETA/GAMMA * F3 / (1 + A*F3)
2038
H = GAMMA * PHI**3 * LOG( 1 + F4 )
2039
FC = ECUNIF + H
2040
2041
C Find correlation energy derivatives
2042
DRSDD = - (THD * RS / DT)
2043
DKFDD = THD * KF / DT
2044
DKSDD = HALF * KS * DKFDD / KF
2045
DZDD(1) = 1 / DT - ZETA / DT
2046
DZDD(2) = - (1 / DT) - ZETA / DT
2047
DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
2048
DO 40 IS = 1,2
2049
DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
2050
DPDD = DPDZ * DZDD(IS)
2051
DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
2052
DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
2053
DF2DD = (- F2) * DF1DD
2054
DADD = (- A) * DF2DD / (F2-1)
2055
DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
2056
DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
2057
DHDD = 3 * H * DPDD / PHI
2058
DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
2059
DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
2060
2061
DO 30 IX = 1,3
2062
DTDGD = (T / GDMT) * GDT(IX) / GDMT
2063
DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
2064
DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )
2065
DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
2066
DFCDGD(IX,IS) = DT * DHDGD
2067
30 CONTINUE
2068
40 CONTINUE
2069
2070
C Find exchange energy and potential
2071
FX = 0
2072
DO 60 IS = 1,2
2073
DS(IS) = MAX( DENMIN, 2 * D(IS) )
2074
GDMS = MAX( GDMIN, 2 * GDM(IS) )
2075
KFS = (3 * PI**2 * DS(IS))**THD
2076
S = GDMS / (2 * KFS * DS(IS))
2077
cea Hammer's RPBE (Hammer, Hansen & Norskov PRB 59 7413 (99)
2078
cea F1 = DEXP( - MU * S**2 / KAPPA)
2079
cea F = 1 + KAPPA * (1 - F1)
2080
cea Following is standard PBE
2081
cea F1 = 1 + MU * S**2 / KAPPA
2082
cea F = 1 + KAPPA - KAPPA / F1
2083
cea (If revPBE Zhang & Yang, PRL 80,890(1998),change PBE's KAPPA to 1.245)
2084
F1 = DEXP( - MU * S**2 / KAPPA)
2085
F = 1 + KAPPA * (1 - F1)
2086
2087
c Note nspin=1 in call to exchng...
2088
2089
CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
2090
FX = FX + DS(IS) * EXUNIF * F
2091
2092
cMVFS The derivatives of F also need to be changed for Hammer's RPBE.
2093
cMVFS DF1DD = 2 * F1 * DSDD * ( - MU * S / KAPPA)
2094
cMVFS DF1DGD= 2 * F1 * DSDGD * ( - MU * S / KAPPA)
2095
cMVFS DFDD = -1 * KAPPA * DF1DD
2096
cMVFS DFDGD = -1 * KAPPA * DFDGD
2097
2098
DKFDD = THD * KFS / DS(IS)
2099
DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
2100
c DF1DD = 2 * (F1-1) * DSDD / S
2101
c DFDD = KAPPA * DF1DD / F1**2
2102
DF1DD = 2* F1 * DSDD * ( - MU * S / KAPPA)
2103
DFDD = -1 * KAPPA * DF1DD
2104
DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
2105
2106
DO 50 IX = 1,3
2107
GDS = 2 * GD(IX,IS)
2108
DSDGD = (S / GDMS) * GDS / GDMS
2109
c DF1DGD = 2 * MU * S * DSDGD / KAPPA
2110
c DFDGD = KAPPA * DF1DGD / F1**2
2111
DF1DGD =2*F1 * DSDGD * ( - MU * S / KAPPA)
2112
DFDGD = -1 * KAPPA * DF1DGD
2113
DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
2114
50 CONTINUE
2115
60 CONTINUE
2116
FX = HALF * FX / DT
2117
2118
C Set output arguments
2119
EX = FX
2120
EC = FC
2121
DO 90 IS = 1,nspin
2122
DEXDD(IS) = DFXDD(IS)
2123
DECDD(IS) = DFCDD(IS)
2124
DO 80 IX = 1,3
2125
DEXDGD(IX,IS) = DFXDGD(IX,IS)
2126
DECDGD(IX,IS) = DFCDGD(IX,IS)
2127
80 CONTINUE
2128
90 CONTINUE
2129
2130
END
2131
SUBROUTINE WCXC( IREL, nspin, Dens, GDens,
2132
. EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
2133
2134
C *********************************************************************
2135
C Implements Wu-Cohen Generalized-Gradient-Approximation.
2136
C Ref: Z. Wu and R. E. Cohen PRB 73, 235116 (2006)
2137
C Written by Marivi Fernandez-Serra, with contributions by
2138
C Julian Gale and Alberto Garcia,
2139
C over the PBEXC subroutine of L.C.Balbas and J.M.Soler.
2140
C September, 2006.
2141
C ******** INPUT ******************************************************
2142
C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
2143
C INTEGER nspin : Number of spin polarizations (1 or 2)
2144
C REAL*8 Dens(nspin) : Total electron density (if nspin=1) or
2145
C spin electron density (if nspin=2)
2146
C REAL*8 GDens(3,nspin) : Total or spin density gradient
2147
C ******** OUTPUT *****************************************************
2148
C REAL*8 EX : Exchange energy density
2149
C REAL*8 EC : Correlation energy density
2150
C REAL*8 DEXDD(nspin) : Partial derivative
2151
C d(DensTot*Ex)/dDens(ispin),
2152
C where DensTot = Sum_ispin( Dens(ispin) )
2153
C For a constant density, this is the
2154
C exchange potential
2155
C REAL*8 DECDD(nspin) : Partial derivative
2156
C d(DensTot*Ec)/dDens(ispin),
2157
C where DensTot = Sum_ispin( Dens(ispin) )
2158
C For a constant density, this is the
2159
C correlation potential
2160
C REAL*8 DEXDGD(3,nspin): Partial derivative
2161
C d(DensTot*Ex)/d(GradDens(i,ispin))
2162
C REAL*8 DECDGD(3,nspin): Partial derivative
2163
C d(DensTot*Ec)/d(GradDens(i,ispin))
2164
C ********* UNITS ****************************************************
2165
C Lengths in Bohr
2166
C Densities in electrons per Bohr**3
2167
C Energies in Hartrees
2168
C Gradient vectors in cartesian coordinates
2169
C ********* ROUTINES CALLED ******************************************
2170
C EXCHNG, PW92C
2171
C ********************************************************************
2172
2173
use precision, only : dp
2174
2175
implicit none
2176
INTEGER IREL, nspin
2177
real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
2178
. DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
2179
2180
C Internal variables
2181
INTEGER
2182
. IS, IX
2183
2184
real(dp)
2185
. A, BETA, D(2), DADD, DECUDD, DENMIN,
2186
. DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
2187
. DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
2188
. DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
2189
. DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
2190
. XWC, DXWCDS, CWC,
2191
. EC, ECUNIF, EX, EXUNIF,
2192
. F, F1, F2, F3, F4, FC, FX, FOUTHD,
2193
. GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
2194
. H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
2195
. TEN81,
2196
. T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
2197
2198
C Lower bounds of density and its gradient to avoid divisions by zero
2199
PARAMETER ( DENMIN = 1.D-12 )
2200
PARAMETER ( GDMIN = 1.D-12 )
2201
2202
C Fix some numerical parameters
2203
PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
2204
. THD=1.D0/3.D0, THRHLF=1.5D0,
2205
. TWO=2.D0, TWOTHD=2.D0/3.D0 )
2206
PARAMETER ( TEN81 = 10.0d0/81.0d0 )
2207
2208
C Fix some more numerical constants
2209
PI = 4 * ATAN(1.D0)
2210
BETA = 0.066725D0
2211
GAMMA = (1 - LOG(TWO)) / PI**2
2212
MU = BETA * PI**2 / 3
2213
KAPPA = 0.804D0
2214
CWC = 0.0079325D0
2215
2216
C Translate density and its gradient to new variables
2217
IF (nspin .EQ. 1) THEN
2218
D(1) = HALF * Dens(1)
2219
D(2) = D(1)
2220
DT = MAX( DENMIN, Dens(1) )
2221
DO 10 IX = 1,3
2222
GD(IX,1) = HALF * GDens(IX,1)
2223
GD(IX,2) = GD(IX,1)
2224
GDT(IX) = GDens(IX,1)
2225
10 CONTINUE
2226
ELSE
2227
D(1) = Dens(1)
2228
D(2) = Dens(2)
2229
DT = MAX( DENMIN, Dens(1)+Dens(2) )
2230
DO 20 IX = 1,3
2231
GD(IX,1) = GDens(IX,1)
2232
GD(IX,2) = GDens(IX,2)
2233
GDT(IX) = GDens(IX,1) + GDens(IX,2)
2234
20 CONTINUE
2235
ENDIF
2236
GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
2237
GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
2238
GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 )
2239
GDMT = MAX( GDMIN, GDMT )
2240
2241
C Find local correlation energy and potential
2242
CALL PW92C( 2, D, ECUNIF, VCUNIF )
2243
2244
C Find total correlation energy
2245
RS = ( 3 / (4*PI*DT) )**THD
2246
KF = (3 * PI**2 * DT)**THD
2247
KS = SQRT( 4 * KF / PI )
2248
ZETA = ( D(1) - D(2) ) / DT
2249
ZETA = MAX( -1.D0+DENMIN, ZETA )
2250
ZETA = MIN( 1.D0-DENMIN, ZETA )
2251
PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
2252
T = GDMT / (2 * PHI * KS * DT)
2253
F1 = ECUNIF / GAMMA / PHI**3
2254
F2 = EXP(-F1)
2255
A = BETA / GAMMA / (F2-1)
2256
F3 = T**2 + A * T**4
2257
F4 = BETA/GAMMA * F3 / (1 + A*F3)
2258
H = GAMMA * PHI**3 * LOG( 1 + F4 )
2259
FC = ECUNIF + H
2260
2261
C Find correlation energy derivatives
2262
DRSDD = - (THD * RS / DT)
2263
DKFDD = THD * KF / DT
2264
DKSDD = HALF * KS * DKFDD / KF
2265
DZDD(1) = 1 / DT - ZETA / DT
2266
DZDD(2) = - (1 / DT) - ZETA / DT
2267
DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
2268
DO 40 IS = 1,2
2269
DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
2270
DPDD = DPDZ * DZDD(IS)
2271
DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
2272
DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
2273
DF2DD = (- F2) * DF1DD
2274
DADD = (- A) * DF2DD / (F2-1)
2275
DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
2276
DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
2277
DHDD = 3 * H * DPDD / PHI
2278
DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
2279
DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
2280
2281
DO 30 IX = 1,3
2282
DTDGD = (T / GDMT) * GDT(IX) / GDMT
2283
DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
2284
DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )
2285
DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
2286
DFCDGD(IX,IS) = DT * DHDGD
2287
30 CONTINUE
2288
40 CONTINUE
2289
2290
C Find exchange energy and potential
2291
FX = 0
2292
DO 60 IS = 1,2
2293
DS(IS) = MAX( DENMIN, 2 * D(IS) )
2294
GDMS = MAX( GDMIN, 2 * GDM(IS) )
2295
KFS = (3 * PI**2 * DS(IS))**THD
2296
S = GDMS / (2 * KFS * DS(IS))
2297
c
2298
c For PBE:
2299
c
2300
c x = MU * S**2
2301
c dxds = 2*MU*S
2302
c
2303
c Wu-Cohen form:
2304
c
2305
XWC= TEN81 * s**2 + (MU- TEN81) *
2306
. S**2 * exp(-S**2) + log(1+ CWC * S**4)
2307
DXWCDS = 2 * TEN81 * S + (MU - TEN81) * exp(-S**2) *
2308
. 2*S * (1 - S*S) + 4 * CWC * S**3 / (1 + CWC * S**4)
2309
c-------------------
2310
2311
F1 = 1 + XWC / KAPPA
2312
F = 1 + KAPPA - KAPPA / F1
2313
c
2314
c Note nspin=1 in call to exchng...
2315
c
2316
CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
2317
FX = FX + DS(IS) * EXUNIF * F
2318
2319
DKFDD = THD * KFS / DS(IS)
2320
DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
2321
DF1DD = DXWCDS * DSDD / KAPPA
2322
DFDD = KAPPA * DF1DD / F1**2
2323
DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
2324
2325
DO 50 IX = 1,3
2326
GDS = 2 * GD(IX,IS)
2327
DSDGD = (S / GDMS) * GDS / GDMS
2328
DF1DGD = DXWCDS * DSDGD / KAPPA
2329
DFDGD = KAPPA * DF1DGD / F1**2
2330
DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
2331
50 CONTINUE
2332
60 CONTINUE
2333
FX = HALF * FX / DT
2334
2335
C Set output arguments
2336
EX = FX
2337
EC = FC
2338
DO 90 IS = 1,nspin
2339
DEXDD(IS) = DFXDD(IS)
2340
DECDD(IS) = DFCDD(IS)
2341
DO 80 IX = 1,3
2342
DEXDGD(IX,IS) = DFXDGD(IX,IS)
2343
DECDGD(IX,IS) = DFCDGD(IX,IS)
2344
80 CONTINUE
2345
90 CONTINUE
2346
2347
END
2348
2349
SUBROUTINE PBESOLXC( IREL, nspin, Dens, GDens,
2350
. EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
2351
2352
C *********************************************************************
2353
C Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation.
2354
C with the revised parameters for solids (PBEsol).
2355
C Ref: J.P.Perdew et al, PRL 100, 136406 (2008)
2356
C Written by L.C.Balbas and J.M.Soler for PBE. December 1996.
2357
C Modified by J.D. Gale for PBEsol. May 2009.
2358
C ******** INPUT ******************************************************
2359
C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
2360
C INTEGER nspin : Number of spin polarizations (1 or 2)
2361
C REAL*8 Dens(nspin) : Total electron density (if nspin=1) or
2362
C spin electron density (if nspin=2)
2363
C REAL*8 GDens(3,nspin) : Total or spin density gradient
2364
C ******** OUTPUT *****************************************************
2365
C REAL*8 EX : Exchange energy density
2366
C REAL*8 EC : Correlation energy density
2367
C REAL*8 DEXDD(nspin) : Partial derivative
2368
C d(DensTot*Ex)/dDens(ispin),
2369
C where DensTot = Sum_ispin( Dens(ispin) )
2370
C For a constant density, this is the
2371
C exchange potential
2372
C REAL*8 DECDD(nspin) : Partial derivative
2373
C d(DensTot*Ec)/dDens(ispin),
2374
C where DensTot = Sum_ispin( Dens(ispin) )
2375
C For a constant density, this is the
2376
C correlation potential
2377
C REAL*8 DEXDGD(3,nspin): Partial derivative
2378
C d(DensTot*Ex)/d(GradDens(i,ispin))
2379
C REAL*8 DECDGD(3,nspin): Partial derivative
2380
C d(DensTot*Ec)/d(GradDens(i,ispin))
2381
C ********* UNITS ****************************************************
2382
C Lengths in Bohr
2383
C Densities in electrons per Bohr**3
2384
C Energies in Hartrees
2385
C Gradient vectors in cartesian coordinates
2386
C ********* ROUTINES CALLED ******************************************
2387
C EXCHNG, PW92C
2388
C ********************************************************************
2389
2390
use precision, only : dp
2391
2392
implicit none
2393
INTEGER IREL, nspin
2394
real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
2395
. DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
2396
2397
C Internal variables
2398
INTEGER
2399
. IS, IX
2400
2401
real(dp)
2402
. A, BETA, D(2), DADD, DECUDD, DENMIN,
2403
. DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
2404
. DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
2405
. DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
2406
. DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
2407
. EC, ECUNIF, EX, EXUNIF,
2408
. F, F1, F2, F3, F4, FC, FX, FOUTHD,
2409
. GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
2410
. H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
2411
. T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
2412
2413
C Lower bounds of density and its gradient to avoid divisions by zero
2414
PARAMETER ( DENMIN = 1.D-12 )
2415
PARAMETER ( GDMIN = 1.D-12 )
2416
2417
C Fix some numerical parameters
2418
PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
2419
. THD=1.D0/3.D0, THRHLF=1.5D0,
2420
. TWO=2.D0, TWOTHD=2.D0/3.D0 )
2421
2422
C Fix some more numerical constants
2423
PI = 4 * ATAN(1.D0)
2424
BETA = 0.046d0
2425
GAMMA = (1 - LOG(TWO)) / PI**2
2426
MU = 10.0d0/81.0d0
2427
KAPPA = 0.804D0
2428
2429
C Translate density and its gradient to new variables
2430
IF (nspin .EQ. 1) THEN
2431
D(1) = HALF * Dens(1)
2432
D(2) = D(1)
2433
DT = MAX( DENMIN, Dens(1) )
2434
DO 10 IX = 1,3
2435
GD(IX,1) = HALF * GDens(IX,1)
2436
GD(IX,2) = GD(IX,1)
2437
GDT(IX) = GDens(IX,1)
2438
10 CONTINUE
2439
ELSE
2440
D(1) = Dens(1)
2441
D(2) = Dens(2)
2442
DT = MAX( DENMIN, Dens(1)+Dens(2) )
2443
DO 20 IX = 1,3
2444
GD(IX,1) = GDens(IX,1)
2445
GD(IX,2) = GDens(IX,2)
2446
GDT(IX) = GDens(IX,1) + GDens(IX,2)
2447
20 CONTINUE
2448
ENDIF
2449
GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
2450
GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
2451
GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 )
2452
GDMT = MAX( GDMIN, GDMT )
2453
2454
C Find local correlation energy and potential
2455
CALL PW92C( 2, D, ECUNIF, VCUNIF )
2456
2457
C Find total correlation energy
2458
RS = ( 3 / (4*PI*DT) )**THD
2459
KF = (3 * PI**2 * DT)**THD
2460
KS = SQRT( 4 * KF / PI )
2461
ZETA = ( D(1) - D(2) ) / DT
2462
ZETA = MAX( -1.D0+DENMIN, ZETA )
2463
ZETA = MIN( 1.D0-DENMIN, ZETA )
2464
PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
2465
T = GDMT / (2 * PHI * KS * DT)
2466
F1 = ECUNIF / GAMMA / PHI**3
2467
F2 = EXP(-F1)
2468
A = BETA / GAMMA / (F2-1)
2469
F3 = T**2 + A * T**4
2470
F4 = BETA/GAMMA * F3 / (1 + A*F3)
2471
H = GAMMA * PHI**3 * LOG( 1 + F4 )
2472
FC = ECUNIF + H
2473
2474
C Find correlation energy derivatives
2475
DRSDD = - (THD * RS / DT)
2476
DKFDD = THD * KF / DT
2477
DKSDD = HALF * KS * DKFDD / KF
2478
DZDD(1) = 1 / DT - ZETA / DT
2479
DZDD(2) = - (1 / DT) - ZETA / DT
2480
DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
2481
DO 40 IS = 1,2
2482
DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
2483
DPDD = DPDZ * DZDD(IS)
2484
DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
2485
DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
2486
DF2DD = (- F2) * DF1DD
2487
DADD = (- A) * DF2DD / (F2-1)
2488
DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
2489
DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
2490
DHDD = 3 * H * DPDD / PHI
2491
DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
2492
DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
2493
2494
DO 30 IX = 1,3
2495
DTDGD = (T / GDMT) * GDT(IX) / GDMT
2496
DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
2497
DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )
2498
DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
2499
DFCDGD(IX,IS) = DT * DHDGD
2500
30 CONTINUE
2501
40 CONTINUE
2502
2503
C Find exchange energy and potential
2504
FX = 0
2505
DO 60 IS = 1,2
2506
DS(IS) = MAX( DENMIN, 2 * D(IS) )
2507
GDMS = MAX( GDMIN, 2 * GDM(IS) )
2508
KFS = (3 * PI**2 * DS(IS))**THD
2509
S = GDMS / (2 * KFS * DS(IS))
2510
F1 = 1 + MU * S**2 / KAPPA
2511
F = 1 + KAPPA - KAPPA / F1
2512
c
2513
c Note nspin=1 in call to exchng...
2514
c
2515
CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
2516
FX = FX + DS(IS) * EXUNIF * F
2517
2518
DKFDD = THD * KFS / DS(IS)
2519
DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
2520
DF1DD = 2 * (F1-1) * DSDD / S
2521
DFDD = KAPPA * DF1DD / F1**2
2522
DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
2523
2524
DO 50 IX = 1,3
2525
GDS = 2 * GD(IX,IS)
2526
DSDGD = (S / GDMS) * GDS / GDMS
2527
DF1DGD = 2 * MU * S * DSDGD / KAPPA
2528
DFDGD = KAPPA * DF1DGD / F1**2
2529
DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
2530
50 CONTINUE
2531
60 CONTINUE
2532
FX = HALF * FX / DT
2533
2534
C Set output arguments
2535
EX = FX
2536
EC = FC
2537
DO 90 IS = 1,nspin
2538
DEXDD(IS) = DFXDD(IS)
2539
DECDD(IS) = DFCDD(IS)
2540
DO 80 IX = 1,3
2541
DEXDGD(IX,IS) = DFXDGD(IX,IS)
2542
DECDGD(IX,IS) = DFCDGD(IX,IS)
2543
80 CONTINUE
2544
90 CONTINUE
2545
2546
END
2547
Loggerhead is a web-based interface for Breezy
Version: 2.0.1