1
Subroutine xc_hcth(tol_rho, xfac, lxfac, nlxfac,
2
, cfac, lcfac, nlcfac,rho, delrho,
3
& Amat, Cmat, nq, ipol, Ex, Ec, qwght,
6
c$Id: xc_hcth.F 20247 2011-04-28 18:58:49Z d3y133 $
12
logical ldew ! [input]
13
logical lcfac, nlcfac, lxfac, nlxfac ! [input]
14
double precision func(*) ![input/output]
15
double precision cfac, xfac ![input]
16
character*4 funcname ! functional name [input]
18
integer ipol ! no. of spin states [input]
19
integer nq ! no. of quadrature pts [input]
20
double precision tol_rho! [input]!threshold on density
21
double precision Ec ! Correlation energy [input/output]
22
double precision Ex ! Exchange energy [input/output]
23
double precision rho(nq,ipol*(ipol+1)/2)! Charge Density [input]
24
double precision delrho(nq,3,ipol) ! Charge Density Gradient[input]
25
double precision qwght(nq) ! Quadrature Weights [input]
26
double precision Amat(nq,ipol) !Sampling Matrices for the XC [output]
27
double precision Cmat(nq,*)!Potential & Energy [output]
30
c F.A.Hamprecht, A.J.Cohen, D.J.Tozer and N.C.Handy,
31
c J. Chem. Phys. 109, 6264-6271 (1998)
32
c subroutine supplied by Fred Hamprecht fah@igc.phys.chem.ethz.ch
35
double precision gammaval
41
double precision hE_x, hE_c
42
double precision dfdrax, dfdrac
43
double precision dfdrbx, dfdrbc
44
double precision dfdzac,dfdzax,dfdza
45
double precision dfdzbc,dfdzbx,dfdzb
49
if(rho(n,1).gt.tol_rho) then
55
gammaval=0.25d0*(delrho(n,1,1)*delrho(n,1,1) +
56
& delrho(n,2,1)*delrho(n,2,1) +
57
& delrho(n,3,1)*delrho(n,3,1))
58
if(gammaval.gt.tol_rho) za=sqrt(gammaval)
60
call hcth(ipol,funcname,
61
* dfdrax,dfdrac, dfdzax,dfdzac,
62
* dfdrbx,dfdrbc, dfdzbx,dfdzbc,
64
& za, zb, hE_x, hE_c, tol_rho)
65
if(ldew) func(n)=func(n)+hE_c*cfac+hE_x*xfac
66
Ec=Ec+hE_c*qwght(n)*cfac
67
Ex=Ex+hE_x*qwght(n)*xfac
68
Amat(n,1) = Amat(n,1)+dfdrac*cfac+dfdrax*xfac
69
dfdza=(dfdzac*cfac+dfdzax*xfac)*0.5d0
70
Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dfdza
75
if(rho(n,1).gt.tol_rho) then
79
gammaval=delrho(n,1,1)*delrho(n,1,1) +
80
& delrho(n,2,1)*delrho(n,2,1) +
81
& delrho(n,3,1)*delrho(n,3,1)
82
if(gammaval.gt.tol_rho) za=sqrt(gammaval)
84
gammaval=delrho(n,1,2)*delrho(n,1,2) +
85
& delrho(n,2,2)*delrho(n,2,2) +
86
& delrho(n,3,2)*delrho(n,3,2)
87
if(gammaval.gt.tol_rho) zb=sqrt(gammaval)
88
call hcth(ipol,funcname,
89
* dfdrax,dfdrac, dfdzax,dfdzac,
90
* dfdrbx,dfdrbc, dfdzbx,dfdzbc,
92
& za, zb, hE_x, hE_c, tol_rho)
93
if(ldew) func(n)=func(n)+hE_c*cfac+hE_x*xfac
94
Ec=Ec+hE_c*qwght(n)*cfac
95
Ex=Ex+hE_x*qwght(n)*xfac
96
Amat(n,1) = Amat(n,1)+dfdrac*cfac+dfdrax*xfac
97
Amat(n,2) = Amat(n,2)+dfdrbc*cfac+dfdrbx*xfac
98
dfdza=dfdzac*cfac+dfdzax*xfac
99
dfdzb=dfdzbc*cfac+dfdzbx*xfac
100
Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dfdza*0.5d0
101
Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + dfdzb*0.5d0
108
Subroutine xc_hcth_d2(tol_rho, xfac, lxfac, nlxfac,
109
& cfac, lcfac, nlcfac,rho, delrho,
110
& Amat, Amat2, Cmat, Cmat2, nq, ipol,
111
& Ex, Ec, qwght, ldew,func,funcname)
113
c$Id: xc_hcth.F 20247 2011-04-28 18:58:49Z d3y133 $
116
#include "errquit.fh"
118
#include "dft2drv.fh"
119
#include "mafdecls.fh"
121
logical ldew ! [input]
122
logical lcfac, nlcfac, lxfac, nlxfac ! [input]
123
double precision func(*) ![input/output]
124
double precision cfac, xfac ![input]
125
character*4 funcname ! functional name [input]
127
integer ipol ! no. of spin states [input]
128
integer nq ! no. of quadrature pts [input]
129
double precision tol_rho! [input]!threshold on density
130
double precision Ec ! Correlation energy [input/output]
131
double precision Ex ! Exchange energy [input/output]
132
double precision rho(nq,ipol*(ipol+1)/2)! Charge Density [input]
133
double precision delrho(nq,3,ipol) ! Charge Density Gradient[input]
134
double precision qwght(nq) ! Quadrature Weights [input]
135
double precision Amat(nq,ipol) !Sampling Matrices for the XC [output]
136
double precision Cmat(nq,*)!Potential & Energy [output]
137
double precision Amat2(nq,NCOL_AMAT2) ! XC functional seconds [output]
138
double precision Cmat2(nq,NCOL_CMAT2) ! XC functional seconds [output]
142
integer l_storage, i_prho, i_pdelrho, i_pAmat, i_pCmat, i_pfunc,
143
& i_qwght_copy, npert
144
double precision ExDum, EcDum
146
c First get the functional and first derivative values
148
call xc_hcth(tol_rho, xfac, lxfac, nlxfac, cfac, lcfac, nlcfac,
149
& rho, delrho, Amat, Cmat, nq, ipol, Ex, Ec, qwght, ldew, func,
152
c Compute the second derivative values by finite difference
154
call xc_setup_fd(tol_rho, rho, delrho, qwght, nq, ipol, .true.,
155
& l_storage, i_prho, i_pdelrho, i_pAmat, i_pCmat, i_pfunc,
158
c Compute functional first derivatives at perturbed density parameter
159
c values - note that the number of points is nq*2*npert and that the
160
c routine is called as unrestricted
163
call xc_hcth(tol_rho, xfac, lxfac, nlxfac, cfac, lcfac, nlcfac,
164
& dbl_mb(i_prho), dbl_mb(i_pdelrho), dbl_mb(i_pAmat),
165
& dbl_mb(i_pCmat), nq*2*npert, ipol, ExDum, EcDum,
166
& dbl_mb(i_qwght_copy), ldew, dbl_mb(i_pfunc), funcname)
167
call xc_make_fd(Amat2, Cmat2, nq, .true., dbl_mb(i_pAmat),
170
c Free temporary storage allocated by xc_setup_fd
172
if (.not.ma_free_heap(l_storage))
173
& call errquit('xc_hcth_d2: cannot pop stack',0, MA_ERR)
178
SUBROUTINE hcth(ipol,functional,
179
* dfdrax,dfdrac, dfdzax,dfdzac,
180
* dfdrbx,dfdrbc, dfdzbx,dfdzbc,
181
1 rhoa, rhob, za, zb, hE_x, hE_c,
184
c subroutine supplied by Fred Hamprecht fah@igc.phys.chem.ethz.ch
185
C SUPPLIED TO THE ROUTINE:
187
C rhoa -- value of rhoalpha at a given grid point
188
C rhob -- value of rhobeta at a given grid point
189
C za -- zeta_alpha, as defined in the TH1 paper (JCP 108 2545),
190
C that is mod(grad(rhoalpha)), a scalar quantity.
191
C zb -- mod(grad(rhobeta))
192
C zab -- zeta_{alpha beta} as defined in the TH1 paper, that is
193
C grad(rhoalpha).grad(rhobeta)
194
C energy -- a boolean variable deciding whether to compute the energy
195
C contribution at the point in space (true) or the
196
C appropriate derivatives (false) needed for the KS matrix
197
C _and_ the energy contribution.
199
C RETURNED FROM THE ROUTINE:
201
C hE_x -- the contribution to the energy at this point in space.
202
C hE_c -- the contribution to the energy at this point in space.
203
C dfdra -- partial functional derivative of F_xc with respect to
205
C dfdrb -- partial functional derivative of F_xc with respect to
207
C dfdza -- partial functional derivative of F_xc with respect to
208
C mod(grad(rhoalpha)), divided by za !!!!!!!!!
212
C This is a consequence of the Cadpac implementation
213
C dfdzb -- partial functional derivative of F_xc with respect to
214
C mod(grad(rhobeta)), divided by zb !!!!!!!!!
218
#include "errquit.fh"
220
double precision rhoa ![input]
221
double precision rhob ![input]
222
double precision za ![input]
223
double precision zb ![input]
224
integer ipol ![input]
225
double precision hE_x ![output]
226
double precision hE_c ![output]
227
double precision tol_rho ![input]
228
double precision dfdrax ![output]
229
double precision dfdrac ![output]
230
double precision dfdrbx ![output]
231
double precision dfdrbc ![output]
232
double precision dfdzax ![output]
233
double precision dfdzac ![output]
234
double precision dfdzbx ![output]
235
double precision dfdzbc ![output]
236
character*4 functional
238
PARAMETER (PI=3.1415926535898D0)
241
PARAMETER(limpow = 4)
242
Cfah limpow is equivalent to "m" in the Becke V paper, that is, the greatest
243
Cfah power of u appearing in the power expansion.
244
Cfah The "Becke V" paper is:
245
Cfah Becke A. D. Density-functional thermochemistry. V.
246
Cfah Systematic optimization of exchange-correlation functionals,
247
Cfah J. Chem. Phys., 1997, 107, 8554-8560
249
c variables passed to hcderiv
251
integer numfunc,nofunc,max_pow_u
252
parameter(numfunc=14)
253
double precision sol((limpow+1)*3), F((limpow+1)*3,4),
254
C & FF((limpow+1)*3,5,4),
257
Cfah sol -- contains the coefficients of the terms in F_xc
258
Cfah convention: sol(1) = c_{x alpha, 0}, c_{x beta, 0}
259
Cfah sol(2) = c_{c alpha alpha, 0}, c_{c beta beta, 0}
260
Cfah sol(3) = c_{c alpha beta, 0}
261
Cfah sol(4) = c_{x alpha, 1}, c_{x beta, 1}
262
Cfah sol(5) = c_{c alpha alpha, 1}, c_{c beta beta, 1}
263
Cfah sol(6) = c_{c alpha beta, 1}
267
Cfah f(5) -- contains the partial first functional derivatives of F_xc with
269
Cfah the four quantities (IN THIS ORDER): ra, rb, za, zb
271
Cfah ff(5,5) contains the second derivatives with
272
Cfah respect to the same five quantities
274
Cfah F_xa -- contains the alpha exchange bit containing the various powers
275
Cfah of u_{x alpha} (eq. (18) of Becke V paper)
278
Cfah F_caa -- contains the alpha parallel spin correlation bit with the powers
279
Cfah of u_{c alpha alpha}
282
Cfah F_cab -- contains the anti-parallel spin correlation bit with the powers
283
Cfah of u_{c alpha beta}
285
Cfah these transformed variables u will be defined and given short-cut names
287
double precision coeffs(3*(limpow+1),numfunc)
288
character*4 funcnam(numfunc)
289
integer maxpow(numfunc)
290
data maxpow / 2, 2 , 2, 4, 4, 4, 4,
291
, 2 , 4 , 4 , 2 , 4 , 4 , 2/
293
data funcnam/'b970','b980','b971','hcth','hcta','h120','h147',
294
, 'b97g','h407','hp14','b972','407p','b973','b97d'/
295
C B97 B98 B97-1 HCTH HCTH-A
298
data (coeffs(1,n),n=1,numfunc)/
299
/ +0.80940d+00,0.790194d0, +0.789518d+00,+0.109320d+01,
300
, +0.109878d+01,1.09163d0, 1.09025d0, 1.1068d0, 1.08184d0,
301
, +0.103161d+01,+0.827642D+00,+1.08018D0,+7.334648D-01,
305
data (coeffs(2,n),n=1,numfunc)/
306
, +0.17370d+00,-0.120163d0,+0.820011d-01,+0.222601d+00,
307
, +0.136823d-01, 0.48951d0, 0.56258d0, 0.4883d0, 1.18777d0,
308
, +0.282414d+01,+0.585808D+00,+0.80302D0,+5.623649D-01,
311
data (coeffs(3,n),n=1,numfunc)/
312
, +0.94540d+00,0.934715d0,+0.955689d+00,+0.729974d+00,
313
, +0.836897d+00,0.51473d0, 0.54235d0, 0.7961d0, 0.58908d0,
314
, +0.821827d-01,+0.999849D+00,+0.73604D0,+1.133830D+00,
317
data (coeffs(4,n),n=1,numfunc)/
318
, +0.50730d+00,0.400271d0,+0.573805d+00,-0.744056d+00,
319
, -0.251173d+01,-0.74720d0, -0.79920d0, -0.8765d0, -0.5183d0,
320
, -0.360781d+00,+0.478400D-01,-0.4117D0,+2.925270D-01,
323
data (coeffs(5,n),n=1,numfunc)/
324
, +0.23487d+01,2.82332d0,+0.271681d+01,-0.338622d-01,
325
, +0.268920d+00,-0.26070d0, -0.01710d0, -2.117d0, -2.4029d0,
326
, +0.318843d-01,-0.691682D+00,-1.0479D0,-1.322980D+00,
329
data (coeffs(6,n),n=1,numfunc)/
330
, +0.74710d+00,1.14105d0,+0.788552d+00,+0.335287d+01,
331
, +0.172051d+01,6.92980d0, 7.01460d0, 5.7060d0, 4.4237d0,
332
, +0.456466d+01,+0.140626D+01,+3.0270D0,-2.811967D+00,
335
data (coeffs(7,n),n=1,numfunc)/
336
, +0.74810d+00,0.832857d0,+0.660975d+00,+0.559920d+01,
337
, +0.156233d-01,5.07830d0, 5.57210d0, 4.2639d0, 3.4256d0,
338
, +0.351994d+01,+0.176125D+01,+2.4368D0,+3.338789D+00,
341
data (coeffs(8,n),n=1,numfunc)/
342
, -0.24868d+01,-2.59412d0,-0.287103d+01,-0.125170d-01,
343
, -0.550769d+00,0.43290d0, -1.30640d0, 2.3235d0, 5.6174d0,
344
, -0.178512d+01,+0.394796D+00,+4.9807D0,+6.359191D+00,
347
data (coeffs(9,n),n=1,numfunc)/
348
, -0.45961d+01,-5.33398d0,-0.547869d+01,-0.115430d+02,
349
, -0.278498d+01,-24.7070d0, -28.3820d0,-14.9820d0,-19.222d0,
350
, -0.135529d+02,-0.744060D+01,-10.075D0,+7.431302D+00,
353
data (coeffs(10,n),n=1,numfunc)/
354
, 0.0000000d0,0.0d000000, 0.00000000d0,-0.678549d+01,
355
, 0.00000000d0,-4.10750d0, -5.86760d0 ,0d0 , -2.6290d0,
356
, -0.495944d+01,0.d0000000000,+1.3890D0,-1.051158D+01,
359
data (coeffs(11,n),n=1,numfunc)/
360
, 0.0000000d0,0.0d000000, 0.00000000d0,-0.802496d+00,
361
, +0.103947d+01,-1.99250d0, 1.05750d0,0d0 , -9.1792d0,
362
, +0.239795d+01,0.d0000000000,-12.890D0,-7.464002D+00,
365
data (coeffs(12,n),n=1,numfunc)/
366
, 0.0000000d0,0.0d000000, 0.00000000d0,+0.808564d+01,
367
, -0.457504d+01, 23.1100d0, 35.0330d0,0d0 , 42.572d0,
368
, +0.133820d+02,0.000000d0,+20.611D0,-1.969342D+00,
371
data (coeffs(13,n),n=1,numfunc)/
372
, 0.0000000d0,0.0d000000, 0.00000000d0,+0.449357d+01,
373
, 0.00000000d0, 1.17170d0, 3.04540d0, 0d0 , 2.2886d0,
374
, +0.241165d+01,0.000000d0,-1.3529D0,+1.060907D+01,
377
data (coeffs(14,n),n=1,numfunc)/
378
, 0.0000000d0,0.0d000000, 0.00000000d0,+0.155396d+01,
379
, 0.00000000d0, 2.48530d0, 0.88540d0, 0d0 , 6.2480d0,
380
, -0.876909d+00,0.000000d0,9.6446D0,+1.827082D+00,
383
data (coeffs(15,n),n=1,numfunc)/
384
, 0.0000000d0,0.0d000000, 0.00000000d0,-0.447857d+01,
385
, 0.00000000d0,-11.3230d0, -20.4280d0, 0d0, -42.005d0,
386
, -0.317493d+01,0.000000d0,-29.418D0,-1.174423D+01,
389
c X +0.1! coeffs for HF exchange
391
C , +0.19430d+00,+0.210000d+00, 0.00000000d0, 0.00000000d0/
406
IF (rhoa .LT. tol_rho.and.rhob.lt.tol_rho) RETURN
407
Cfah numerical cutoff: if the density is too low, its contribution is
409
nofunc = -1 ! take care of compiler warnings
411
if(functional.eq.funcnam(n)) nofunc=n
413
if(nofunc.eq.-1) call errquit('xchcth: cant pair funcname ',0,
415
max_pow_u=maxpow(nofunc)
417
sol(n)=coeffs(n,nofunc)
420
C please refer to these coeffs as THCH1/iterate-e750-g500-v1-m4-n4
422
c sol( 1) = 0.109320D+01
423
c sol( 2) = 0.222601D+00
424
c sol( 3) = 0.729974D+00
425
c sol( 4) = -0.744056D+00
426
c sol( 5) = -0.338622D-01
427
c sol( 6) = 0.335287D+01
428
c sol( 7) = 0.559920D+01
429
c sol( 8) = -0.125170D-01
430
c sol( 9) = -0.115430D+02
431
c sol(10) = -0.678549D+01
432
c sol(11) = -0.802496D+00
433
c sol(12) = 0.808564D+01
434
c sol(13) = 0.449357D+01
435
c sol(14) = 0.155396D+01
436
c sol(15) = -0.447857D+01
438
CALL hcderiv(max_pow_u,ipol,
446
c DO n = 1, (max_pow_u+1)*3
447
c dfdra = dfdra + F(n,1)
448
c dfdrb = dfdrb + F(n,2)
449
c if(za.gt.tol_rho) dfdza = dfdza + F(n,3) / za
450
c if(zb.gt.tol_rho) dfdzb = dfdzb + F(n,4) / zb
453
c DO n = 1, (max_pow_u+1)*3
454
c dfdra = dfdra + F(n,1)
456
c if(za.gt.tol_rho) then
457
c DO n = 1, (max_pow_u+1)*3
458
c dfdza = dfdza + F(n,3) / za
462
Cfah big thanks to NCH: cadpac requires df/(za * dza), NOT
465
hE_x = hE_x + F_xc (n*3 + 1)
466
hE_c = hE_c + F_xc (n*3 + 2) + F_xc (n*3 + 3)
467
dfdrax = dfdrax + F(n*3+1,1)
468
dfdrac = dfdrac + F(n*3+2,1) + F(n*3+3,1)
469
if(za.gt.tol_rho) then
470
dfdzax = dfdzax + F(n*3+1,3) / za
471
dfdzac = dfdzac + (F(n*3+2,3)+F(n*3+3,3)) / za
474
dfdrbx = dfdrbx + F(n*3+1,2)
475
dfdrbc = dfdrbc + F(n*3+2,2) + F(n*3+3,2)
476
if(zb.gt.tol_rho) then
477
dfdzbx = dfdzbx + F(n*3+1,4) / zb
478
dfdzbc = dfdzbc + (F(n*3+2,4)+F(n*3+3,4)) / zb
484
SUBROUTINE hcderiv(max_pow_u,ipol,
491
c subroutine supplied by Fred Hamprecht fah@igc.phys.chem.ethz.ch
493
INTEGER max_pow_u,ipol
496
double precision f_xc((limpow + 1)*3)
497
double precision f((limpow+1)*3,4)
498
C, ff((limpow+1)*3,5,4)
499
double precision rhoa, rhob, za, zb,tol_rho
501
Cfah COMMON/special/h_atom
503
DOUBLE PRECISION sol((limpow+1)*3)
506
DOUBLE PRECISION dF_xa(4)
507
DOUBLE PRECISION dF_xb(4)
508
DOUBLE PRECISION dF_caa(4)
509
DOUBLE PRECISION dF_cbb(4)
510
DOUBLE PRECISION dF_cab(4)
511
Cfah these are the first derivatives of the terms of F_xc with respect to
512
Cfah the 4 quantities. the index
513
Cfah runs over the particular partial derivatives of each term.
514
Cfah More explicitly: these are the partial functional derivatives of
515
Cfah F_XXX with respect to rhoa, rhob, za and zb.
517
c DOUBLE PRECISION d2F_xa(4,4)
518
c DOUBLE PRECISION d2F_xb(4,4)
519
c DOUBLE PRECISION d2F_caa(4,4)
520
c DOUBLE PRECISION d2F_cbb(4,4)
521
c DOUBLE PRECISION d2F_cab(4,4)
523
Cfah these are the first derivatives of the different transformed variables
524
Cfah u with respect to rhoa, rhob, za and zb. These different derivatives
525
Cfah with respect to these 4 quantities named above are stored in these
529
PARAMETER (Pi = 3.1415926535898D0)
531
DOUBLE PRECISION s_a2, s_b2, s_avg2, u_caa, u_cbb, u_cab
532
DOUBLE PRECISION du_caa_by_drhoa, du_caa_by_dza, du_cbb_by_drhob
533
DOUBLE PRECISION du_cbb_by_dzb, du_cab_by_drhoa, du_cab_by_drhob
534
DOUBLE PRECISION du_cab_by_dza, du_cab_by_dzb
535
C, du_caa_by_drhoa_dza
536
C DOUBLE PRECISION du_caa_by_dza_dza, du_cbb_by_dzb_dzb
537
C DOUBLE PRECISION du_cbb_by_drhob_dzb, du_cab_by_drhoa_dza
538
C DOUBLE PRECISION du_cab_by_drhoa_dzb, du_cab_by_drhob_dza
539
C DOUBLE PRECISION du_cab_by_drhob_dzb, du_cab_by_dza_dza,
541
C DOUBLE PRECISION du_cab_by_dzb_dzb
542
DOUBLE PRECISION rsa, rsa12, rsa32, rsa21, rsb,
544
DOUBLE PRECISION rsab, rsab12, rsab32, rsab21
545
DOUBLE PRECISION drsa_by_drhoa, drsb_by_drhob, drsab_by_drhoa
546
DOUBLE PRECISION drsab_by_drhob
547
DOUBLE PRECISION zeta, dzeta_by_drhoa, dzeta_by_drhob
548
DOUBLE PRECISION fzeta, dfzeta_by_dzeta,
550
DOUBLE PRECISION e_crsab1, e_crsab0, a_crsab
551
DOUBLE PRECISION e_crsabzeta, de_crsa1_by_drsa, de_crsb1_by_drsb
552
DOUBLE PRECISION da_crsab_by_drsab, de_crsab0_by_drsab
553
DOUBLE PRECISION de_crsab1_by_drsab, de_crsabzeta_by_drsab
554
DOUBLE PRECISION de_crsabzeta_by_dzeta, e_caa, e_cbb, e_cab,
555
& de_caa_by_drhoa, de_cbb_by_drhob, de_cab_by_drhoa,
557
& c_naa, c_nbb, c_nab
558
DOUBLE PRECISION F_xs,F_xs0 ! this is a function which is called.
559
DOUBLE PRECISION dF_xs_by_drhos, dF_xs_by_dzs,
560
,dF_xs_by_drhos0,dF_xs_by_drhos1, dF_xs_by_dzs1
563
double precision x1,x2,x3,x4
564
double precision e_crs1
565
double precision drsbydrh
566
double precision decrsdrs
568
parameter (eps=1d-19)
570
C F_xs computes HCTH contribution to exchange Energy
571
C using Dirac functional as LDA part
573
C F_xs(n, sol(), rhoa, za)
575
F_xs(n1, x1, x2, x3) =
576
= (-3.D0*x1*(3.D0/Pi)**(1.D0/3.D0)*x2**(4.D0/3.D0)*
577
- ((0.004D0*x3**2.D0)/(0.004D0*x3**2.D0 +
578
- x2**(8.D0/3.D0)))**n1)/(2.D0*2.D0**(2.D0/3.D0))
580
= (-3.D0*x1*(3.D0/Pi)**(1.D0/3.D0)*x2**(4.D0/3.D0)
581
- )/(2.D0*2.D0**(2.D0/3.D0))
583
C dF_xs_by_drhos computes dE_x/drho derivative
584
Cfah computes the derivative of the term with u^n of the exchange part of
585
Cfah F_xc with respect to rho of the same spin.
586
Cfah n -- the power of u involved in this term
587
Cfah c_xs -- the coefficient c_xs(n) of the term of spin s with the
588
Cfah power n of u; is NOT passed over as an array.
589
Cfah rhos -- rhosigma, that is, either rhoalpha or rhobeta
590
Cfah zs -- mod(grad(rhosigma)), again for alpha or beta
591
C usage dF_xs_by_drhos(n, c_xs, rhos, zs)
593
dF_xs_by_drhos(n1, x1, x2, x3) = -(x1*(6.D0/Pi)**
594
* (1.D0/3.D0)*x2**(1.D0/3.D0)*((0.004D0*x3*x3)/
595
/ (x2**(8.D0/3.D0)+0.004D0*x3*x3))**n1)+
596
+ (x1*0.008D0*n1*(6.D0/Pi)**(1.D0/3.D0)*
597
* x2**3.D0*x3*x3*((0.004D0*x3*x3)/
598
/ (x2**(8.D0/3.D0)+0.004D0*x3*x3))**(-1 + n1))/
599
/ (x2**(8.D0/3.D0)+0.004D0*x3*x3)**2
600
dF_xs_by_drhos0(x1, x2, x3) = -(x1*(6.D0/Pi)**
601
* (1.D0/3.D0)*x2**(1.D0/3.D0))
602
dF_xs_by_drhos1(x1, x2, x3) = -(x1*(6.D0/Pi)**
603
* (1.D0/3.D0)*x2**(1.D0/3.D0)*((0.004D0*x3*x3)/
604
/ (x2**(8.D0/3.D0)+0.004D0*x3*x3)))+
605
+ (x1*0.008D0*(6.D0/Pi)**(1.D0/3.D0)*
607
/ (x2**(8.D0/3.D0)+0.004D0*x3*x3)**2
609
C F_xc with respect to zs
610
Cfah see above (function dF_xs_by_drhos) for definition of the
612
C usage dF_xs_by_dzs (n, c_xs, rhos, zs)
614
dF_xs_by_dzs(n1, x1, x2, x3) =
615
= (-3.d0*x1*n1*(3.D0/Pi)**(1.D0/3.D0)*
616
* x2**(4.D0/3.D0)*((0.004D0*x3*x3)/
617
/ (x2**(8.D0/3.D0)+ 0.004D0*x3*x3))**(-1+n1)*
618
* ((-0.000032D0*x3*x3*x3)/(x2**(8.D0/3.D0)+
619
+ 0.004D0*x3*x3)**2+(0.008D0*x3)/
620
/ (x2**(8.D0/3.D0)+0.004D0*x3*x3)))/
621
/ (2.D0**(5.D0/3.D0))
622
dF_xs_by_dzs1(x1, x2, x3) =
623
= (-3.d0*x1*(3.D0/Pi)**(1.D0/3.D0)*
625
* ((-0.000032D0*x3*x3*x3)/(x2**(8.D0/3.D0)+
626
+ 0.004D0*x3*x3)**2+(0.008D0*x3)/
627
/ (x2**(8.D0/3.D0)+0.004D0*x3*x3)))/
628
/ (2.D0**(5.D0/3.D0))
629
c dF_xs_by_dzs0(x1, x2, x3) = 0d0
632
c e_crsa1 = e_crs1(rsa12,rsa,rsa32,rsa21)
634
e_crs1(x1,x2,x3,x4) = -0.03108999999999999d0*
635
* dlog(1.d0 + 32.16468317787069d0/
636
/ (14.1189d0*x1+6.1977d0*x2 + 3.3662d0*x3 +
637
+ 0.6251699999999999d0*x4))*(1.d0 + 0.20548d0*x2)
638
drsbydrh(x1) = -((1.d0/x1)**(4.D0/3.D0)/
639
- (6.d0**(2.D0/3.D0)*Pi**(1.D0/3.D0)))
640
c usage decrsdrs(rsa,rsa12,rsa21,rsa32)
641
decrsdrs(x1,x2,x3,x4) = ((1.d0 + 0.20548d0*x1)*
642
- (6.1977d0 + 7.05945d0/x2 + 1.25034d0*x1+5.0493d0*x2))/
643
- ((6.1977d0*x1+14.1189d0*x2+0.6251699999999999d0*x3 +
644
+ 3.3662d0*x4)**2d0*(1.d0 + 32.16468317787069d0/
645
- (6.1977d0*x1+14.1189d0*x2+0.6251699999999999d0*x3+
646
+ 3.3662d0*x4))) - 0.006388373199999999d0*
647
- dlog(1.d0 + 32.16468317787069d0/(6.1977d0*x1 + 14.1189d0*x2 +
648
- 0.6251699999999999d0*x3 + 3.3662d0*x4))
651
DO n = 1, (max_pow_u+1)*3
653
Cfah later on, n has a different meaning: n as power of u, not
654
Cfah as number of the coefficient.
662
C DO n = 1, (max_pow_u+1)*3
667
C d2F_caa(j,k) = 0.D0
668
C d2F_cbb(j,k) = 0.D0
669
C d2F_cab(j,k) = 0.D0
672
DO j = 1, (max_pow_u+1)*3
676
Cfah --------------------------------------------------------------
678
Cfah call the expensive correlation parts here just once, and store their
679
Cfah values in a temporary variable. Then compute the actual F_c derivatives
680
Cfah with the various powers of u.
685
if(za.gt.tol_rho) s_a2 = za**2.D0 / rhoa**(8.D0/3.D0)
687
if(zb.gt.tol_rho) s_b2 = zb**2.D0 / rhob**(8.D0/3.D0)
688
s_avg2 = 0.5D0*(s_a2 + s_b2)
690
u_caa = 0.2D0*s_a2/(1.D0+0.2D0*s_a2)
691
u_cbb = 0.2D0*s_b2/(1.D0+0.2D0*s_b2)
693
u_cab = 0.006D0*s_avg2/(1.d0+0.006D0*s_avg2)
694
if(rhoa.gt.tol_rho) then
695
rsa = ((3.d0/Pi)**(1.D0/3.D0)*
696
- (1.d0/rhoa)**(1.D0/3.D0))/
698
rsa12 = rsa**(1.D0/2.D0)
699
rsa32 = rsa**(3.D0/2.D0)
708
if(rhob.gt.tol_rho) then
709
rsb = ((3.d0/Pi)**(1.D0/3.D0)*
710
- (1.d0/rhob)**(1.D0/3.D0))/
712
rsb12 = rsb**(1.D0/2.D0)
713
rsb32 = rsb**(3.D0/2.D0)
718
if(rhob.gt.tol_rho) then
719
e_crsb1 = e_crs1(rsb12,rsb,rsb32,rsb21)
720
de_crsb1_by_drsb = decrsdrs(rsb,rsb12,rsb21,rsb32)
723
de_crsb1_by_drsb = 0d0
726
du_cbb_by_drhob = (-1.6D0*zb**2*rhob**(5.D0/3.D0))/
727
- (3.d0*(0.2D0*zb**2 +
728
- rhob**(8.D0/3.D0))**2)
729
du_cbb_by_dzb = (2*0.2D0*zb*rhob**(8.D0/3.D0))/
730
- (0.2D0*zb**2 + rhob**(8.D0/3.D0))**2
732
du_cab_by_drhoa = (-16*0.006D0*za*za*rhoa**(5.D0/3.D0)*
733
- rhob**(16.D0/3.D0))/
734
- (3.d0*(0.006D0*zb**2*rhoa**(8.D0/3.D0) +
735
- 0.006D0*za*za*rhob**(8.D0/3.D0) +
736
- 2.d0*rhoa**(8.D0/3.D0)*
737
- rhob**(8.D0/3.D0))**2)
738
du_cab_by_dza = (4*0.006D0*za*rhoa**(8.D0/3.D0)*
739
- rhob**(16.D0/3.D0))/
740
- (0.006D0*zb**2*rhoa**(8.D0/3.D0) +
741
- 0.006D0*za*za*rhob**(8.D0/3.D0) +
742
- 2.d0*rhoa**(8.D0/3.D0)*
743
- rhob**(8.D0/3.D0))**2
744
du_cab_by_drhob = (-16*0.006D0*zb**2*rhob**(5.D0/3.D0)*
745
- rhoa**(16.D0/3.D0))/
746
- (3.d0*(0.006D0*za*za*rhob**(8.D0/3.D0) +
747
- 0.006D0*zb**2*rhoa**(8.D0/3.D0) +
748
- 2.d0*rhob**(8.D0/3.D0)*
749
- rhoa**(8.D0/3.D0))**2)
750
du_cab_by_dzb = (4*0.006D0*zb*rhoa**(16.D0/3.D0)*
751
- rhob**(8.D0/3.D0))/
752
- (0.006D0*zb**2*rhoa**(8.D0/3.D0) +
753
- 0.006D0*za*za*rhob**(8.D0/3.D0) +
754
- 2.d0*rhoa**(8.D0/3.D0)*
755
- rhob**(8.D0/3.D0))**2
757
du_cab_by_drhoa = 0d0
759
du_cab_by_drhob = 0d0
762
drsb_by_drhob = drsbydrh(rhob)
765
de_crsb1_by_drsb = 0d0
766
du_cbb_by_drhob = 0d0
768
du_cab_by_drhoa = 0d0
769
du_cab_by_drhob = 0d0
775
rsab = ((3.d0/Pi)**(1.D0/3.D0)*
776
- (1.d0/rho)**(1.D0/3.D0))/
778
rsab12 = rsab**(1.D0/2.D0)
779
rsab32 = rsab**(3.D0/2.D0)
782
zeta = (rhoa-rhob)/rho
783
if(zeta.lt.-1d0) zeta=-1d0
784
if(zeta.gt.1d0) zeta=1d0
786
if(abs(1d0-zeta).gt.eps) then
787
fzeta = (-2.d0 + sign(1d0,1.d0 - zeta)*(abs(1.d0 - zeta))**
789
- (1.d0 + zeta)**(4.D0/3.D0))/
790
- (-2.d0 + 2.d0*2.d0**(1.D0/3.D0))
798
if(rhoa.gt.tol_rho) then
799
e_crsa1 = e_crs1(rsa12,rsa,rsa32,rsa21)
803
if(rho.gt.tol_rho) then
804
e_crsab1 = e_crs1(rsab12,rsab,rsab32,rsab21)
805
e_crsab0 = -0.062182d0*dlog(1.d0 +
806
- 16.0818243221511d0/
807
- (7.595699999999999d0*rsab12 +
810
- 0.49294d0*rsab21))*
811
- (1.d0 + 0.2137d0*rsab)
813
a_crsab = 0.03377399999999999d0*
814
- dlog(1.d0 + 29.60857464321667/
815
- (10.35699999999999d0*rsab12 +
816
- 3.623099999999999d0*rsab +
818
- 0.49671d0*rsab21))*
819
- (1.d0 + 0.11125d0*rsab)
826
e_crsabzeta = e_crsab0+a_crsab*fzeta*(1.d0-zeta**4)/1.709921D0+
827
- (e_crsab1-e_crsab0)*fzeta*zeta**4
831
e_cab = rho*e_crsabzeta - rhoa*e_crsa1 - rhob*e_crsb1
832
if(rhoa.gt.tol_rho) then
833
du_caa_by_drhoa = (-1.6D0*za*za*rhoa**(5.D0/3.D0))/
835
- rhoa**(8.D0/3.D0))**2)
836
du_caa_by_dza = (2*0.2D0*za*rhoa**(8.D0/3.D0))/
837
- (0.2D0*za*za + rhoa**(8.D0/3.D0))**2
839
du_caa_by_drhoa = 0d0
845
Cfah Second derivatives are not required by cadpac.
846
Cfah du_caa_by_drhoa_dza = (16*0.2D0*za*rhoa**(5.D0/3.D0)*
847
Cfah - (0.2D0*za**2 - rhoa**(8.D0/3.D0)))/
848
Cfah - (3.*(0.2D0*za**2 +
849
Cfah - rhoa**(8.D0/3.D0))**3)
850
Cfah du_cbb_by_drhob_dzb = (16*0.2D0*zb*rhob**(5.D0/3.D0)*
851
Cfah - (0.2D0*zb**2 - rhob**(8.D0/3.D0)))/
852
Cfah - (3.*(0.2D0*zb**2 +
853
Cfah - rhob**(8.D0/3.D0))**3)
855
Cfah du_caa_by_dza_dza = (2*0.2D0*rhoa**(8.D0/3.D0)*
856
Cfah - (-3*0.2D0*za**2 + rhoa**(8.D0/3.D0))
858
Cfah - (0.2D0*za**2 + rhoa**(8.D0/3.D0))**3
859
Cfah du_cbb_by_dzb_dzb = (2*0.2D0*rhob**(8.D0/3.D0)*
860
Cfah - (-3*0.2D0*zb**2 + rhob**(8.D0/3.D0))
862
Cfah - (0.2D0*zb**2 + rhob**(8.D0/3.D0))**3
864
Cfah du_cab_by_drhoa_dza = (-32*0.006D0*rhoa**(5.D0/3.D0)*
865
Cfah - (0.006D0*za*zb**2*
866
Cfah - rhoa**(8.D0/3.D0)*
867
Cfah - rhob**(16.D0/3.D0) -
868
Cfah - 0.006D0*za**3*rhob**8 +
869
Cfah - 2*za*rhoa**(8.D0/3.D0)*rhob**8))/
870
Cfah - (3.*(0.006D0*zb**2*rhoa**(8.D0/3.D0) +
871
Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) +
872
Cfah - 2*rhoa**(8.D0/3.D0)*
873
Cfah - rhob**(8.D0/3.D0))**3)
874
Cfah du_cab_by_drhoa_dzb = (64*0.006D0**2*za**2*zb*
875
Cfah - rhoa**(13.D0/3.D0)*
876
Cfah - rhob**(16.D0/3.D0))/
877
Cfah - (3.*(0.006D0*zb**2*rhoa**(8.D0/3.D0) +
878
Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) +
879
Cfah - 2*rhoa**(8.D0/3.D0)*
880
Cfah - rhob**(8.D0/3.D0))**3)
881
Cfah du_cab_by_drhob_dza = (64*0.006D0**2*za*zb**2*
882
Cfah - rhoa**(16.D0/3.D0)*
883
Cfah - rhob**(13.D0/3.D0))/
884
Cfah - (3.*(0.006D0*zb**2*rhoa**(8.D0/3.D0) +
885
Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) +
886
Cfah - 2*rhoa**(8.D0/3.D0)*
887
Cfah - rhob**(8.D0/3.D0))**3)
888
Cfah du_cab_by_drhob_dzb = (-32*0.006D0*rhob**(5.D0/3.D0)*
889
Cfah - (-(0.006D0*zb**3*rhoa**8) +
890
Cfah - 0.006D0*za**2*zb*
891
Cfah - rhoa**(16.D0/3.D0)*
892
Cfah - rhob**(8.D0/3.D0) +
893
Cfah - 2*zb*rhoa**8*rhob**(8.D0/3.D0)))/
894
Cfah - (3.*(0.006D0*zb**2*rhoa**(8.D0/3.D0) +
895
Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) +
896
Cfah - 2*rhoa**(8.D0/3.D0)*
897
Cfah - rhob**(8.D0/3.D0))**3)
898
Cfah du_cab_by_dza_dza = (4*0.006D0*(0.006D0*zb**2*
899
Cfah - rhoa**(16.D0/3.D0)*
900
Cfah - rhob**(16.D0/3.D0) -
901
Cfah - 3*0.006D0*za**2*rhoa**(8.D0/3.D0)*
902
Cfah - rhob**8 + 2*rhoa**(16.D0/3.D0)*rhob**8
904
Cfah - (0.006D0*zb**2*rhoa**(8.D0/3.D0) +
905
Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) +
906
Cfah - 2*rhoa**(8.D0/3.D0)*
907
Cfah - rhob**(8.D0/3.D0))**3
908
Cfah du_cab_by_dza_dzb = (-16*0.006D0**2*za*zb*
909
Cfah - rhoa**(16.D0/3.D0)*
910
Cfah - rhob**(16.D0/3.D0))/
911
Cfah - (0.006D0*zb**2*rhoa**(8.D0/3.D0) +
912
Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) +
913
Cfah - 2*rhoa**(8.D0/3.D0)*
914
Cfah - rhob**(8.D0/3.D0))**3
915
Cfah du_cab_by_dzb_dzb = (4*0.006D0*rhoa**(16.D0/3.D0)*
916
Cfah - rhob**(8.D0/3.D0)*
917
Cfah - (-3*0.006D0*zb**2*rhoa**(8.D0/3.D0) +
918
Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) +
919
Cfah - 2*rhoa**(8.D0/3.D0)*
920
Cfah - rhob**(8.D0/3.D0)))/
921
Cfah - (0.006D0*zb**2*rhoa**(8.D0/3.D0) +
922
Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) +
923
Cfah - 2*rhoa**(8.D0/3.D0)*
924
Cfah - rhob**(8.D0/3.D0))**3
926
if(rhoa.gt.tol_rho) then
927
drsa_by_drhoa = drsbydrh(rhoa)
931
if(rho.gt.tol_rho) then
932
drsab_by_drhoa = drsbydrh(rho)
933
drsab_by_drhob = drsab_by_drhoa
939
dzeta_by_drhoa = 2.d0*rhob/rho**2
940
dzeta_by_drhob = -2.d0*rhoa/rho**2
942
dfzeta_by_dzeta = ((-4.d0*sign(1d0,1.d0 - zeta)*
943
* (abs(1.d0 - zeta))**(1.D0/3.D0))/
944
- 3.d0 + (4.d0*(1.d0 + zeta)**
946
- (-2.d0 + 2d0*2**(1.D0/3.D0))
949
da_crsab_by_drsab = (-1.d0*(1.d0 + 0.11125d0*rsab)*
950
- (3.623099999999999d0 +
951
- 5.178499999999999d0/rsab12 +
952
- 0.99342d0*rsab + 1.32039d0*rsab12))/
953
- ((1.d0 + 29.60857464321667d0/
954
- (3.623099999999999d0*rsab +
955
- 10.35699999999999d0*rsab12 +
956
- 0.49671d0*rsab21 + 0.88026d0*rsab32))*
957
- (3.623099999999999d0*rsab +
958
- 10.35699999999999d0*rsab12 +
959
- 0.49671d0*rsab21 + 0.880260*rsab32)**2) +
960
- 0.003757357499999999d0*
961
- dlog(1.d0 + 29.60857464321667/
962
- (3.623099999999999d0*rsab +
963
- 10.35699999999999d0*rsab12 +
964
- 0.49671d0*rsab21 + 0.88026d0*rsab32))
966
de_crsab0_by_drsab = (1.d0*(1.d0 + 0.2137d0*rsab)*
967
- (3.5876d0 + 3.797849999999999d0/rsab12 +
968
- 0.98588d0*rsab + 2.4573*rsab12))/
969
- ((3.5876d0*rsab + 7.595699999999999d0*rsab12 +
970
- 0.49294d0*rsab21 + 1.6382d0*rsab32)**2*
971
- (1.d0 + 16.0818243221511d0/
972
- (3.5876d0*rsab + 7.595699999999999d0*rsab12 +
973
- 0.49294d0*rsab21 + 1.6382d0*rsab32))) -
974
- 0.01328829339999999d0*
975
- dlog(1.d0 + 16.0818243221511d0/
976
- (3.5876*rsab + 7.595699999999999d0*rsab12 +
977
- 0.49294d0*rsab21 + 1.6382d0*rsab32))
978
if(rhoa.gt.tol_rho) then
979
de_crsa1_by_drsa = decrsdrs(rsa,rsa12,rsa21,rsa32)
980
de_crsab1_by_drsab = decrsdrs(rsab,rsab12,rsab21,rsab32)
982
de_crsa1_by_drsa = 0d0
983
de_crsab1_by_drsab = 0d0
987
de_crsabzeta_by_drsab = 1.124999956683108D0*(1.d0 - zeta**4)*
988
- (-2.d0+sign(1d0,1.d0-zeta)*(abs(1.d0-zeta))**(4.D0/3.D0)+
989
- (1.d0 + zeta)**(4.D0/3.D0))*
990
- da_crsab_by_drsab +
991
- de_crsab0_by_drsab +
993
- (- de_crsab0_by_drsab +
994
- de_crsab1_by_drsab )
996
de_crsabzeta_by_dzeta = 1.499999942244144D0*(-1.d0 + zeta**4)*
997
- (sign(1d0,1.d0 - zeta)*(abs(1.d0 - zeta))**(1.D0/3.D0) -
998
- (1.d0 + zeta)**(1.D0/3.D0))*a_crsab -
999
- 4.499999826732434D0*zeta**3*
1000
- (-2.d0 + sign(1d0,1.d0-zeta)*abs(1.d0-zeta)**(4.D0/3.D0)+
1001
- (1.d0 + zeta)**(4.D0/3.D0))*a_crsab +
1002
- (2*zeta**4*sign(1d0,1.d0-zeta)*(abs(1.d0-zeta)**(1.D0/3.D0)-
1003
- (1.d0 + zeta)**(1.D0/3.D0))*
1004
- (e_crsab0 - e_crsab1))/
1005
- (3.*(-1.d0 + 2**(1.D0/3.D0))) +
1006
- 4*fzeta*(-e_crsab0 +
1009
Cfah this is with application of the chain rule; I keep it that general
1010
Cfah because this way, I only have to define one "G".
1011
de_caa_by_drhoa = e_crsa1 + rhoa*de_crsa1_by_drsa*drsa_by_drhoa
1012
de_cbb_by_drhob = e_crsb1 + rhob*de_crsb1_by_drsb*drsb_by_drhob
1014
de_cab_by_drhoa = -e_crsa1 +
1016
- rhoa*de_crsa1_by_drsa*
1018
- rho*(de_crsabzeta_by_drsab*
1020
- de_crsabzeta_by_dzeta*
1023
de_cab_by_drhob = -e_crsb1 +
1025
- rhob*de_crsb1_by_drsb*
1027
- rho*(de_crsabzeta_by_drsab*
1029
- de_crsabzeta_by_dzeta*
1034
Cfah Here starts the big outer loop over the powers u
1036
c_naa = sol((n*3) + 2)
1038
c_nab = sol((n*3) + 3)
1040
Cfah construction of the F_xc itself
1041
Cfah -------------------------------
1042
IF (rhoa.GT.tol_rho) THEN
1044
F_xc(1) = F_xs0 (sol(1), rhoa, za)
1046
F_xc(n*3+1) = F_xs (n, sol((n*3) + 1), rhoa, za)
1049
if(u_caa.gt.tol_rho)then
1051
F_xc(2) = e_caa*c_naa
1053
F_xc(n*3+2) = e_caa*u_caa**n*c_naa
1057
IF (rhob.GT.tol_rho) THEN
1059
F_xc(1) = F_xc(1)+F_xs0(sol(1), rhob, zb)
1061
F_xc(n*3+1) = F_xc(n*3+1)+F_xs(n, sol((n*3) + 1), rhob, zb)
1064
if(u_cbb.gt.tol_rho) then
1066
F_xc(2) = F_xc(2)+e_cbb*c_nbb
1068
F_xc(n*3+2) = F_xc(n*3+2)+e_cbb*u_cbb**n*c_nbb
1072
if(u_cab.gt.tol_rho) then
1074
F_xc(3) = e_cab*c_nab
1076
F_xc(n*3+3) = e_cab*u_cab**n*c_nab
1080
Cfah print*, 'in deriv:', e_cab, u_cab, c_nab
1082
Cfah First Derivatives
1083
Cfah ---------------------
1085
if(za.gt.tol_rho)then
1087
dF_xa(1) = dF_xs_by_drhos0 ( sol(1), rhoa, za)
1090
dF_xa(1) = dF_xs_by_drhos1 (sol(4), rhoa, za)
1091
dF_xa(3) = dF_xs_by_dzs1 ( sol(4), rhoa, za)
1093
dF_xa(1) = dF_xs_by_drhos (n, sol((n*3) + 1), rhoa, za)
1094
dF_xa(3) = dF_xs_by_dzs (n, sol((n*3) + 1), rhoa, za)
1098
if(zb.gt.tol_rho) then
1100
dF_xb(2) = dF_xs_by_drhos0(sol(1), rhob, zb)
1103
dF_xb(2) = dF_xs_by_drhos1(sol(4), rhob, zb)
1104
dF_xb(4) = dF_xs_by_dzs1 ( sol(4), rhob, zb)
1106
dF_xb(2) = dF_xs_by_drhos (n, sol((n*3) + 1), rhob, zb)
1107
dF_xb(4) = dF_xs_by_dzs (n, sol((n*3) + 1), rhob, zb)
1111
if(u_caa.gt.tol_rho) then
1113
dF_caa(1) = c_naa*de_caa_by_drhoa
1116
dF_caa(1) = c_naa*u_caa*
1117
* de_caa_by_drhoa+c_naa*e_caa*du_caa_by_drhoa
1118
dF_caa(3) = c_naa*e_caa*du_caa_by_dza
1120
dF_caa(1) = c_naa*u_caa**n*
1121
* de_caa_by_drhoa+c_naa*n*e_caa*u_caa**(-1+n)*
1123
dF_caa(3) = c_naa*n*e_caa*u_caa**(-1+n)*du_caa_by_dza
1127
if(u_cbb.gt.tol_rho) then
1129
dF_cbb(2) = c_nbb*de_cbb_by_drhob
1132
dF_cbb(2) = c_nbb*u_cbb*de_cbb_by_drhob +
1133
- c_nbb*e_cbb*du_cbb_by_drhob
1134
dF_cbb(4) = c_nbb*e_cbb*du_cbb_by_dzb
1136
dF_cbb(2) = c_nbb*u_cbb**n*de_cbb_by_drhob +
1137
- c_nbb*n*e_cbb*u_cbb**(-1 + n)*du_cbb_by_drhob
1138
dF_cbb(4) = c_nbb*n*e_cbb*u_cbb**(-1+n)*du_cbb_by_dzb
1143
if(u_cab.gt.tol_rho) then
1145
dF_cab(1) = c_nab*de_cab_by_drhoa
1146
dF_cab(2) = c_nab*de_cab_by_drhob
1150
dF_cab(1) = c_nab*u_cab*
1151
* de_cab_by_drhoa+c_nab*n*e_cab*du_cab_by_drhoa
1152
dF_cab(2) = c_nab*u_cab*
1153
- de_cab_by_drhob+c_nab*n*e_cab*du_cab_by_drhob
1154
dF_cab(3) = c_nab*e_cab*du_cab_by_dza
1155
dF_cab(4) = c_nab*n*e_cab*du_cab_by_dzb
1157
dF_cab(1) = c_nab*u_cab**n*
1158
* de_cab_by_drhoa+c_nab*n*e_cab*u_cab**(-1+n)*
1160
dF_cab(2) = c_nab*u_cab**n*
1161
- de_cab_by_drhob+c_nab*n*e_cab*u_cab**(-1+n)*
1163
dF_cab(3) = c_nab*n*e_cab*
1164
- u_cab**(-1+n)*du_cab_by_dza
1165
dF_cab(4) = c_nab*n*e_cab*
1166
- u_cab**(-1+n)*du_cab_by_dzb
1170
Cfah Second Derivatives
1171
Cfah ------------------
1173
Cfah d2F_xa(1,1) = d2F_xs_by_drhos_drhos (n, sol((n*3) + 1),
1175
Cfah see comment below, for the (2,2) term.
1176
Cfah d2F_xa(1,2) = 0
1177
Cfah d2F_xa(1,3) = d2F_xs_by_drhos_dzs (n, sol((n*3) + 1), rhoa, za)
1178
Cfah d2F_xa(1,4) = 0
1179
Cfah d2F_xa(2,2) = 0
1180
Cfah d2F_xa(2,3) = 0
1181
Cfah d2F_xa(2,4) = 0
1182
Cfah d2F_xa(3,3) = d2F_xs_by_dzs_dzs (n, sol((n*3) + 1), rhoa, za)
1183
Cfah d2F_xa(3,4) = 0
1184
Cfah d2F_xa(4,4) = 0
1186
Cfah for alpha spin, elements are non-zero when both indices are odd;
1187
Cfah for beta spin, elements are non-zero when both indices are even.
1188
Cfah the matrix is symmetric, and the upper triangle contains the
1189
Cfah 10 elements given above and below.
1191
Cfah d2F_xb(1,1) = 0
1192
Cfah d2F_xb(1,2) = 0
1193
Cfah d2F_xb(1,3) = 0
1194
Cfah d2F_xb(1,4) = 0
1195
Cfah d2F_xb(2,2) = d2F_xs_by_drhos_drhos (n, sol((n*3) + 1),
1197
Cfah this term is NOT zero, but needs not be evaluated since we don't
1198
Cfah need it for the construction of v (cf. routine "va" in the fit
1200
Cfah d2F_xb(2,3) = 0
1201
Cfah d2F_xb(2,4) = d2F_xs_by_drhos_dzs (n, sol((n*3) + 1), rhob, zb)
1202
Cfah d2F_xb(3,3) = 0
1203
Cfah d2F_xb(3,4) = 0
1204
Cfah d2F_xb(4,4) = d2F_xs_by_dzs_dzs (n, sol((n*3) + 1), rhob, zb)
1207
Cfah d2F_caa(1,1) = !=0, but not needed
1208
Cfah d2F_caa(1,2) = 0.D0 (not needed)
1209
Cfah d2F_caa(1,3) = c_naa*n*u_caa**(-1 + n)*
1210
Cfah - de_caa_by_drhoa*
1211
Cfah - du_caa_by_dza +
1212
Cfah - c_naa*(-1 + n)*n*e_caa*
1213
Cfah - u_caa**(-2 + n)*
1214
Cfah - du_caa_by_dza*
1215
Cfah - du_caa_by_drhoa +
1216
Cfah - c_naa*n*e_caa*u_caa**(-1 + n)*
1217
Cfah - du_caa_by_drhoa_dza
1218
Cfah d2F_caa(1,4) = 0.D0
1219
Cfah d2F_caa(2,2) = 0.D0 (not needed)
1220
Cfah d2F_caa(2,3) = 0.D0
1221
Cfah d2F_caa(2,4) = 0.D0
1222
Cfah d2F_caa(3,3) = c_naa*n*e_caa*u_caa**(-2 + n)*
1223
Cfah - ((-1 + n)*du_caa_by_dza**
1225
Cfah - du_caa_by_dza_dza)
1226
Cfah d2F_caa(3,4) = 0.D0
1227
Cfah d2F_caa(4,4) = 0.D0
1230
Cfah d2F_cbb(1,1) = 0.D0 (not needed)
1231
Cfah d2F_cbb(1,2) = 0.D0 (not needed)
1232
Cfah d2F_cbb(1,3) = 0.D0
1233
Cfah d2F_cbb(1,4) = 0.D0
1234
Cfah d2F_cbb(2,2) = !=0, but not needed
1235
Cfah d2F_cbb(2,3) = 0.D0
1236
Cfah d2F_cbb(2,4) = c_nbb*n*u_cbb**(-1 + n)*
1237
Cfah - de_cbb_by_drhob*
1238
Cfah - du_cbb_by_dzb +
1239
Cfah - c_nbb*(-1 + n)*n*e_cbb*
1240
Cfah - u_cbb**(-2 + n)*
1241
Cfah - du_cbb_by_dzb*
1242
Cfah - du_cbb_by_drhob +
1243
Cfah - c_nbb*n*e_cbb*u_cbb**(-1 + n)*
1244
Cfah - du_cbb_by_drhob_dzb
1245
Cfah d2F_cbb(3,3) = 0.D0
1246
Cfah d2F_cbb(3,4) = 0.D0
1247
Cfah d2F_cbb(4,4) = c_nbb*n*e_cbb*u_cbb**(-2 + n)*
1248
Cfah - ((-1 + n)*du_cbb_by_dzb**
1250
Cfah - du_cbb_by_dzb_dzb)
1252
Cfah d2F_cab(1,1) = not needed
1253
Cfah d2F_cab(1,2) = not needed
1254
Cfah d2F_cab(1,3) = c_nab*n*u_cab**(-2 + n)*
1255
Cfah - ((-1 + n)*e_cab*
1256
Cfah - du_cab_by_dza
1257
Cfah - *du_cab_by_drhoa +
1258
Cfah - u_cab*(de_cab_by_drhoa*
1259
Cfah - du_cab_by_dza + e_cab*
1260
Cfah - du_cab_by_drhoa_dza
1262
Cfah d2F_cab(1,4) = c_nab*n*u_cab**(-2 + n)*
1263
Cfah - ( (-1 + n)*e_cab*
1264
Cfah - du_cab_by_dzb
1265
Cfah - *du_cab_by_drhoa +
1267
Cfah - (de_cab_by_drhoa*
1268
Cfah - du_cab_by_dzb +
1270
Cfah - du_cab_by_drhoa_dzb))
1271
Cfah d2F_cab(2,2) = not needed
1272
Cfah d2F_cab(2,3) = c_nab*n*u_cab**(-2 + n)*
1273
Cfah - ((-1 + n)*e_cab*
1274
Cfah - du_cab_by_dza
1275
Cfah - *du_cab_by_drhob +
1277
Cfah - (de_cab_by_drhob*
1278
Cfah - du_cab_by_dza +
1280
Cfah - du_cab_by_drhob_dza))
1281
Cfah d2F_cab(2,4) = c_nab*n*u_cab**(-2 + n)*
1282
Cfah - ((-1 + n)*e_cab*
1283
Cfah - du_cab_by_dzb
1284
Cfah - *du_cab_by_drhob +
1285
Cfah - u_cab*(de_cab_by_drhob*
1286
Cfah - du_cab_by_dzb + e_cab*
1287
Cfah - du_cab_by_drhob_dzb ))
1288
Cfah d2F_cab(3,3) = c_nab*n*e_cab*
1289
Cfah - u_cab**(-2 + n)*
1290
Cfah - ((-1 + n)*du_cab_by_dza**2 +
1292
Cfah - du_cab_by_dza_dza)
1293
Cfah d2F_cab(3,4) = c_nab*n*e_cab*
1294
Cfah - u_cab**(-2 + n)*
1295
Cfah - ((-1 + n)*du_cab_by_dzb*
1296
Cfah - du_cab_by_dza
1298
Cfah - du_cab_by_dza_dzb)
1299
Cfah d2F_cab(4,4) = c_nab*n*e_cab*
1300
Cfah - u_cab**(-2 + n)*
1301
Cfah - ((-1 + n)*du_cab_by_dzb**2 +
1303
Cfah - du_cab_by_dzb_dzb)
1305
Cfah here, the second derivatives are completed (Schwartz's rule:
1306
Cfah df/(dadb) = df/(dbda)
1309
Cfah d2F_xa(j,i) = d2F_xa(i,j)
1310
Cfah d2F_xb(j,i) = d2F_xb(i,j)
1311
Cfah d2F_caa(j,i) = d2F_caa(i,j)
1312
Cfah d2F_cbb(j,i) = d2F_cbb(i,j)
1313
Cfah d2F_cab(j,i) = d2F_cab(i,j)
1317
Cfah test for zero densities (as in beta part of H atom):
1318
IF (rhob.LT.tol_rho) THEN
1324
Cfah d2F_xb(i,j) = 0.D0
1325
Cfah d2F_cbb(i,j) = 0.D0
1326
Cfah d2F_cab(i,j) = 0.D0
1331
IF (rhoa.LT.tol_rho) THEN
1337
Cfah d2F_xa(i,j) = 0.D0
1338
Cfah d2F_caa(i,j) = 0.D0
1339
Cfah d2F_cab(i,j) = 0.D0
1345
Cfah Sum up all the partial derivatives with respect to the same function
1346
Cfah of terms containing different powers of u with the help of the big outer
1349
Cfah have the partial derivative
1352
F(n*3+1,i) = dF_xa(i) + dF_xb(i)
1353
F(n*3+2,i) = dF_caa(i) +dF_cbb(i)
1354
F(n*3+3,i) = dF_cab(i)
1356
Cfah FF(n*3+1,i,j) = d2F_xa(i,j) + d2F_xb(i,j)
1357
Cfah FF(n*3+2,i,j) = d2F_caa(i,j) + d2F_cbb(i,j)
1358
Cfah FF(n*3+3,i,j) = d2F_cab(i,j)
1362
Cfah these partial derivatives have not been computed because they are
1363
Cfah zero since we don't have a gradrhoagradrhob term in the Becke V functional
1368
Cfah FF(n*3+1,i,5) = 0
1369
Cfah FF(n*3+2,i,5) = 0
1370
Cfah FF(n*3+3,i,5) = 0
1371
Cfah FF(n*3+1,5,i) = 0
1372
Cfah FF(n*3+2,5,i) = 0
1373
Cfah FF(n*3+3,5,i) = 0
1382
Cfah-----------------------------------------------------------
added
removed
.pc/09_backported_6.1.1_fixes.patch/src/nwdft/xc/xc_hcth.F
