5
C> \file nwxc_x_hse08.F
6
C> The range separated HSE exchange enhancement factor
14
C> \brief Evaluate the range separated HSE exchange enhancement factor
16
C> This routine calculates the HSE exchange enhancement factor [1,2].
20
C> [1] J. Heyd, G.E. Scuceria, M. Ernzerhof, "Hybrid functionals based
21
C> on a screened Coulomb potential", J. Chem. Phys. <b>118</b>,
22
C> 8207 (2003), DOI: <a href="http://dx.doi.org/10.1063/1.1564060">
23
C> 10.1063/1.1564060</a>.
25
C> [2] J. Heyd, G.E. Scuceria, M. Ernzerhof, "Erratum: Hybrid
26
C> functionals based on a screened Coulomb potential", J. Chem. Phys.
27
C> <b>124</b>, 219906 (2006), DOI:
28
C> <a href="http://dx.doi.org/10.1063/1.2204597">
29
C> 10.1063/1.2204597</a>.
32
Subroutine nwxc_x_HSE08(cam_omega,ipol,rho,s,Fxhse,
35
Subroutine nwxc_x_HSE08_d2(cam_omega,ipol,rho,s,Fxhse,
36
& d10Fxhse,d01Fxhse,d20Fxhse,d02Fxhse,d11Fxhse)
42
c HSE evaluates the Heyd et al. Screened Coulomb
45
c Calculates the enhancement factor
48
double precision cam_omega
49
double precision rho,s,Fxhse,d10Fxhse,d01Fxhse
52
double precision d20Fxhse,d02Fxhse,d11Fxhse
55
double precision A,B,C,D,E
56
double precision ha2,ha3,ha4,ha5,ha6,ha7
57
double precision hb1,hb2,hb3,hb4,hb5,hb6,hb7,hb8,hb9
58
double precision smax,strans,sconst
60
double precision zero,one,two,three,four,five,six,seven,eight
61
double precision nine,ten
62
double precision fifteen,sixteen
64
double precision H,hnum,hden
65
double precision d1H,d1hnum,d1hden
66
double precision s2,s3,s4,s5,s6,s7,s8,s9
67
double precision Fs, d1Fs
68
double precision zeta, lambda, eta, kf, nu, chi, lambda2
69
double precision d1zeta,d1lambda,d1eta,d1nu,d1chi,d1lambda2
70
double precision EGs,d1EGs
71
double precision nu2,L2,L3,nu3,nu4,nu5,nu6
72
double precision Js,Ks,Ms,Ns
73
double precision d1Js,d1Ks,d1Ms,d1Ns
75
double precision tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8
76
double precision tmp9,tmp10,tmp11,tmp12,tmp13,tmp14,tmp15
77
double precision Fxhse1,Fxhse2,Fxhse3,Fxhse4,Fxhse5,Fxhse6
78
double precision d1Fxhse1,d1Fxhse2,d1Fxhse3,d1Fxhse4,d1Fxhse5
79
double precision d1Fxhse6,d1Fxhse7
81
double precision r42,r27,r12,r15,r14,r18,r20,r30,r56,r72
82
double precision r16,r32,r24,r48,r11,r64,r35
83
double precision pi,pi2,srpi,s02
84
double precision f12,f13,f32,f52,f72,f92
86
double precision d2H,d2hnum,d2hden
88
double precision d2zeta,d2lambda,d2eta,d2nu,d2chi,d2lambda2
89
double precision d2EGs
90
double precision d2Js,d2Ks,d2Ms,d2Ns
91
double precision d2Fxhse1,d2Fxhse2,d2Fxhse3,d2Fxhse4,d2Fxhse5
92
double precision d2Fxhse6,d2Fxhse7
93
double precision d11Fxhse1,d11Fxhse2,d11Fxhse3,d11Fxhse4
94
double precision d11Fxhse5,d11Fxhse6,d11Fxhse7
99
c Constants for HJS hole
102
& / 7.57211D-1,-1.06364D-1,-1.18649D-1,
103
& 6.09650D-1,-4.77963D-2 /
105
c Constants for fit of H(s) (PBE hole)
106
c Taken from JCTC_5_754 (2009)
108
Data ha2,ha3,ha4,ha5,ha6,ha7
109
& / 1.59941D-2,8.52995D-2,-1.60368D-1,1.52645D-1,
110
& -9.71263D-2,4.22061D-2 /
112
Data hb1,hb2,hb3,hb4,hb5,hb6,hb7,hb8,hb9
113
& / 5.33319D0,-12.4780D0,11.0988D0,-5.11013D0,
114
& 1.71468D0,-6.10380D-1,3.07555D-1,-7.70547D-2,
118
c Whole numbers used during evaluation
120
Data zero,one,two,three,four,five,six,seven,eight,nine,ten
121
& / 0D0,1D0,2D0,3D0,4D0,5D0,6D0,7D0,8D0,9D0,10D0 /
123
Data r11,r12,r14,r15,r16,r18,r20,r24,r27,r30,r32
124
& / 11D0,12D0,14D0,15D0,16D0,18D0,20D0,24D0,27d0,30D0,32D0 /
126
Data r35,r42,r48,r56,r64,r72
127
& / 35D0,42D0,48D0,56D0,64D0,72D0 /
129
c Fractions used during evaluation
146
c Calculate prelim variables
159
c Calculate H(s) the model exhange hole
161
hnum = ha2*s2 + ha3*s3 + ha4*s4 + ha5*s5 + ha6*s6 + ha7*s7
162
hden = one + hb1*s + hb2*s2 + hb3*s3 + hb4*s4 + hb5*s5 +
163
& hb6*s6 + hb7*s7 + hb8*s8 + hb9*s9
167
c Calculate helper variables
173
kf = (three*pi2*rho)**f13
175
kf = (six*pi2*rho)**f13
178
chi = nu/dsqrt(lambda+nu**two)
179
lambda2 = (one+chi)*(lambda+nu**two)
182
c Calculate F(H(s)) for the model exhange hole
184
Fs = one-s2/(r27*C*(one+s02))-zeta/(two*C)
189
EGs = -(two/five)*C*Fs*lambda - (four/r15)*B*lambda**two -
190
& (six/five)*A*lambda**three -
191
& (four/five)*srpi*lambda**(seven/two) -
192
& (r12/five)*(lambda**(seven/two))*(dsqrt(zeta)-dsqrt(eta))
195
c Calculate the denominators needed
199
Js = (dsqrt(zeta+nu2)+dsqrt(eta+nu2))*(dsqrt(zeta+nu2)+nu)
200
Ks = (dsqrt(zeta+nu2)+dsqrt(eta+nu2))*(dsqrt(eta+nu2)+nu)
201
Ms = (dsqrt(zeta+nu2)+dsqrt(lambda+nu2))*(dsqrt(lambda+nu2)+nu)
202
Ns = (dsqrt(eta+nu2)+dsqrt(lambda+nu2))*(dsqrt(lambda+nu2)+nu)
205
c The final value for the enhancement factor is
208
tmp2 = one + (nine/eight)*chi + (three/eight)*chi**two
209
Fxhse1 = A*(zeta/Js + eta/Ks)
210
Fxhse2 = -(four/nine)*B/lambda2
211
Fxhse3 = -(four/nine)*C*Fs*tmp1/lambda2**two
212
Fxhse4 = -(eight/nine)*EGs*tmp2/lambda2**three
213
Fxhse5 = two*zeta*dlog(one -D/Ms)
214
Fxhse6 = -two*eta*dlog(one -(D-A)/Ns)
216
Fxhse = Fxhse1+Fxhse2+Fxhse3+Fxhse4+Fxhse5+Fxhse6
218
c Calculate the first derivative of H with respect to the
219
c reduced density gradient, s.
221
d1hnum = two*ha2*s + three*ha3*s2 + four*ha4*s3 +
222
& five*ha5*s4 + six*ha6*s5 + seven*ha7*s6
224
d1hden = hb1 + two*hb2*s +three*hb3*s2 + four*hb4*s3 +
225
& five*hb5*s4 + six*hb6*s5 + seven*hb7*s6 +
226
& eight*hb8*s7 + nine*hb9*s8
227
d1H = (hden*d1hnum -hnum*d1hden)/hden**two
230
c calculate first derivative of variables needed with respect to s
232
d1zeta = two*s*H + s2*d1H
235
d1chi = -f12*nu*d1zeta/(lambda + nu2)**f32
236
d1lambda2 = d1chi*(lambda + nu**two) + (one+chi)*d1lambda
237
!d1lambda2 = (d1lambda*(one-chi)+lambda*d1chi)/(one-chi)**two
240
c calculate the first derivative of Fs with respect to s
242
d1Fs = -two*s/(r27*C*(one+s02)**two) - d1zeta/(two*C)
245
c Calculate the first derivate of EGs with respect to s
247
d1EGs = -(two/five)*C*(d1Fs*lambda + Fs*d1lambda) -
248
& (eight/r15)*B*lambda*d1lambda -
249
& (r18/five)*A*lambda*lambda*d1lambda -
250
& (r14/five)*srpi*d1lambda*lambda**f52 -
251
& (r42/five)*(lambda**f52)*
252
& d1lambda*(dsqrt(zeta)-dsqrt(eta))-
253
& (six/five)*(lambda**(seven/two))*
254
& (d1zeta/dsqrt(zeta)-d1eta/dsqrt(eta))
257
c Calculate the first derivate of denominators needed with respect
260
tmp1 = (dsqrt(zeta+nu2)+nu)/(dsqrt(eta+nu2))
261
tmp2 = (dsqrt(eta+nu2)+nu)/(dsqrt(zeta+nu2))
263
d1Js = f12*d1zeta*(two+tmp1+tmp2)
266
tmp3 = (dsqrt(zeta+nu2)+nu)/(dsqrt(lambda+nu2))
267
tmp4 = (dsqrt(lambda+nu2)+nu)/(dsqrt(zeta+nu2))
268
d1Ms = f12*d1zeta*(two +tmp3+tmp4)
270
tmp5 = (dsqrt(lambda+nu2)+nu)/(dsqrt(eta+nu2))
271
tmp6 = (dsqrt(eta+nu2)+nu)/(dsqrt(lambda+nu2))
272
d1Ns = f12*d1zeta*(two + tmp5+tmp6)
275
c Calculate the derivative of the 08-Fxhse with respect to s
278
L3 = lambda2*lambda2*lambda2
279
d1Fxhse1 = A*( (Js*d1zeta - zeta*d1Js)/(Js*Js) +
280
& (Ks*d1zeta - eta*d1Ks)/(Ks*Ks) )
282
d1Fxhse2 = (four/nine)*B*d1lambda2/L2
284
tmp9 = d1lambda2/lambda2
285
tmp7 = d1Fs - two*Fs*tmp9
289
d1Fxhse3 = -(four*C/(nine*L2))*(tmp7*tmp8+tmp10)
291
c d1Fxhse3 = -(four/nine)*(C/(L2*L2))*
292
c & (L2*d1Fs - two*lambda2*Fs*d1lambda2 +
293
c & f12*(L2*(d1Fs*chi+Fs*d1chi) -
294
c & two*lambda2*chi*Fs*d1chi))
296
c tmp4 = d1chi/(one+chi) + d1lambda/(lambda+nu2)
297
c tmp4 = (one-chi)*d1lambda2/lambda
298
c tmp1 = (eight/three+three*chi+chi*chi)*tmp4
299
c tmp2 = d1chi + (two/three)*chi*d1chi
300
c tmp3 = (eight/nine + chi + f13*chi*chi)
301
c d1FXhse4 = ((tmp1-tmp2)*EGs-tmp3*d1EGs)/L3
303
tmp7 = one + (nine/eight)*chi+(three/eight)*chi*chi
304
tmp8 = (nine/eight)*d1chi + (six/eight)*chi*d1chi
306
d1Fxhse4 = -(eight/(nine*L3))*((d1EGs-three*EGs*tmp9)*tmp7
308
c d1Fxhse4 = -(eight/nine)*(L3*d1EGs -
309
c & three*EGs*L2*d1lambda2)/(L3*L3) -
310
c & (L3*(d1EGs*chi + EGs*d1chi) -
311
c & three*EGs*chi*L2*d1lambda2)/(L3*L3)-
312
c & (L3*(d1EGs*chi*chi+two*EGs*chi*d1chi)-
313
c & three*EGs*chi*chi*L2*d1lambda2)/(three*L3*L3)
315
d1Fxhse5 = two*d1zeta*dlog(one-D/Ms) +
316
& two*zeta*D*d1Ms/(Ms*Ms*(one-D/Ms))
318
d1Fxhse6 = -two*d1eta*dlog(one- (D-A)/Ns) -
319
& two*eta*(D-A)*d1Ns/(Ns*Ns*(one-(D-A)/Ns))
321
d10Fxhse = d1Fxhse1+d1Fxhse2+d1Fxhse3+d1Fxhse4+d1Fxhse5+d1Fxhse6
323
c Calculate the derivative of 08-Fxhse with respect to nu
327
d1Fxhse1 = -((A*(nu + dsqrt(eta + nu2))*zeta)/
328
& (dsqrt(eta + nu2)*dsqrt(nu2 + zeta)*
329
& (nu + dsqrt(nu2 + zeta))*
330
& (dsqrt(eta + nu2) + dsqrt(nu2 + zeta))))
332
d1Fxhse2 = -((A*eta*(nu/dsqrt(eta + nu2) + nu/
333
& dsqrt(nu2 + zeta)))/
334
& ((nu + dsqrt(eta + nu2))*
335
& (dsqrt(eta + nu2) + dsqrt(nu2 + zeta))**two)) -
336
& (A*eta*(one + nu/dsqrt(eta + nu2)))/
337
& ((nu + dsqrt(eta + nu2))**two*
338
& (dsqrt(eta + nu2) + dsqrt(nu2 + zeta)))
340
d1Fxhse3 = (four*B)/(nine*(lambda + nu2)**(f32))
342
d1Fxhse4 = (two*C*Fs)/(three*(lambda + nu2)**(f52))
344
d1Fxhse5 = (five*EGs*(lambda**two + four*nu3*
345
& (nu + dsqrt(lambda + nu2)) +
346
& lambda*nu*(five*nu + three*dsqrt(lambda + nu2))))/
347
& (three*(lambda + nu2)**four*(nu + dsqrt(lambda + nu2))**three)
349
d1Fxhse6 = (two*D*zeta*(nu + dsqrt(nu2 + zeta)))/
350
& (dsqrt(lambda + nu2)*dsqrt(nu2 + zeta)*
351
& (-D + lambda + (nu + dsqrt(lambda + nu2))*
352
& (nu + dsqrt(nu2 + zeta))))
354
d1Fxhse7 = (two*(A - D)*eta*(nu + dsqrt(eta + nu2)))/
355
& (dsqrt(eta + nu2)*dsqrt(lambda + nu2)*
356
& (A - D + lambda + nu2 + nu*dsqrt(eta + nu2) +
357
& nu*dsqrt(lambda + nu2) +
358
& dsqrt(eta + nu2)*dsqrt(lambda + nu2)))
361
c d1Fxhse1 = -((A*zeta*(nu/dsqrt(eta + nu2) +
362
c & nu/dsqrt(nu2 + zeta)))/((nu + dsqrt(nu2 + zeta))*
363
c & (dsqrt(eta + nu2) + dsqrt(nu2 + zeta))**two)) -
364
c & (A*zeta*(one + nu/dsqrt(nu2 + zeta)))/
365
c & ((nu + dsqrt(nu2 + zeta))**two*(dsqrt(eta + nu2) +
366
c & dsqrt(nu2 + zeta)))
368
c d1Fxhse2 = -((A*eta*(nu/dsqrt(eta + nu2) +
369
c & nu/dsqrt(nu2 + zeta)))/((nu + dsqrt(eta + nu2))*
370
c & (dsqrt(eta + nu2) + dsqrt(nu2 + zeta))**two)) -
371
c & (A*eta*(one + nu/dsqrt(eta + nu2)))/
372
c & ((nu + dsqrt(eta + nu2))**two*(dsqrt(eta + nu2) +
373
c & dsqrt(nu2 + zeta)))
375
c d1Fxhse3 = (four*B*(-(nu2/(lambda + nu2)**(three/two)) +
376
c & one/dsqrt(lambda + nu2)))/(nine*(lambda + nu2)*
377
c & (one + nu/dsqrt(lambda + nu2))**two) +
378
c & (eight*B*nu)/(nine*(lambda + nu2)**two*
379
c & (one + nu/dsqrt(lambda + nu2)))
381
c d1Fxhse4 = (eight*C*Fs*(-(nu2/(lambda + nu2)**(three/two)) +
382
c & one/dsqrt(lambda + nu2))*
383
c & (one + nu/(two*dsqrt(lambda + nu2))))/(nine*(lambda + nu2)**two*
384
c & (one + nu/dsqrt(lambda + nu2))**three) -
385
c & (four*C*Fs*(-nu2/(two*(lambda + nu2)**(three/two)) +
386
c & one/(two*dsqrt(lambda + nu2))))/(nine*(lambda + nu2)**two*
387
c & (one + nu/dsqrt(lambda + nu2))**two) +
388
c & (r16*C*Fs*nu*(one + nu/(two*dsqrt(lambda + nu2))))/
389
c & (nine*(lambda + nu2)**three*
390
c & (one + nu/dsqrt(lambda + nu2))**two)
392
c d1Fxhse5 = (-eight*EGs*((-three*nu3)/(four*(lambda + nu2)**two) -
393
c & (nine*nu2)/(eight*(lambda + nu2)**(three/two)) +
394
c & (three*nu)/(four*(lambda + nu2)) +
395
c & nine/(eight*dsqrt(lambda + nu2))))/(nine*(lambda + nu2)**three*
396
c & (one + nu/dsqrt(lambda + nu2))**three) +
397
c & (eight*EGs*(-(nu2/(lambda + nu2)**(three/two)) +
398
c & one/dsqrt(lambda + nu2))*
399
c & (one + (three*nu2)/(eight*(lambda + nu2)) +
400
c & (nine*nu)/(eight*dsqrt(lambda + nu2))))/
401
c & (three*(lambda + nu2)**three*
402
c & (one + nu/dsqrt(lambda + nu2))**four)
403
c & + (r16*EGs*nu*(one + (three*nu2)/(eight*(lambda + nu2)) +
404
c & (nine*nu)/(eight*dsqrt(lambda + nu2))))/
405
c & (three*(lambda + nu2)**four*
406
c & (one + nu/dsqrt(lambda + nu2))**three)
408
c d1Fxhse6 = (two*zeta*((D*(nu/dsqrt(lambda + nu2) +
409
c & nu/dsqrt(nu2 + zeta)))/
410
c & ((nu + dsqrt(lambda + nu2))*(dsqrt(lambda + nu2) +
411
c & dsqrt(nu2 + zeta))**two) +
412
c & (D*(one + nu/dsqrt(lambda + nu2)))/
413
c & ((nu + dsqrt(lambda + nu2))**two*(dsqrt(lambda + nu2) +
414
c & dsqrt(nu2 + zeta)))))/
415
c & (one - D/((nu + dsqrt(lambda + nu2))*
416
c & (dsqrt(lambda + nu2) + dsqrt(nu2 + zeta))))
418
c d1Fxhse7 = (-two*eta*(((-A + D)*(nu/dsqrt(eta + nu2) +
419
c & nu/dsqrt(lambda + nu2)))/
420
c & ((nu + dsqrt(lambda + nu2))*(dsqrt(eta + nu2) +
421
c & dsqrt(lambda + nu2))**two) +
422
c & ((-A + D)*(one + nu/dsqrt(lambda + nu2)))/
423
c & ((nu + dsqrt(lambda + nu2))**two*(dsqrt(eta + nu2) +
424
c & dsqrt(lambda + nu2)))))/
425
c & (one - (-A + D)/((nu + dsqrt(lambda + nu2))*
426
c & (dsqrt(eta + nu2) + dsqrt(lambda + nu2))))
429
d01Fxhse = d1Fxhse1+d1Fxhse2+d1Fxhse3+d1Fxhse4+d1Fxhse5+d1Fxhse6+
434
c Calculate the second derivative of H
436
d2hnum = two*ha2+six*ha3*s+r12*ha4*s2+r20*ha5*s3+
437
& r30*ha6*s4 + r42*ha7*s5
439
d2hden = two*hb2+six*hb3*s+r12*hb4*s2+r20*hb5*s3+
440
& r30*hb6*s4+r42*hb7*s5+r56*hb8*s6+r72*hb9*s7
442
d2H = (hden*(hden*d2hnum-hnum*d2hden) -
443
& two*(hden*d1hnum-hnum*d1hden)*d1hden)/hden**three
446
c calculate second derivative of variables needed
448
d2zeta = two*H +four*s*d1H + s2*d2H
452
d2chi = -f12*(nu/(lambda+nu2)**(f32))*
453
& (-f32*d1zeta*d1zeta/(lambda+nu2) +d2zeta)
455
d2lambda2 =(one-chi)*(d2lambda-d2lambda*chi+lambda*d2chi)+
456
& two*d1chi*(d1lambda-d1lambda*chi+lambda*d1chi)
457
d2lambda2 = d2lambda2/(one-chi)**three
460
c calculate the second derivative of Fs
462
d2Fs = -(two/(r27*C))*(one-three*s02)/
463
& ((one+s02)**three)-d2zeta/(two*C)
466
c Calculate the second derivate of EGs
468
tmp7 = -(two/five)*C*(d2Fs*lambda+two*d1Fs*d1lambda+Fs*d2lambda)
469
tmp8 = -(eight/r15)*B*(d1lambda*d1lambda+lambda*d2lambda)
470
tmp9 =-(r18/five)*A*lambda*(two*d1lambda*d1lambda+lambda*d2lambda)
471
tmp10 = -(r14/five)*srpi*(f52*(lambda**f32)*d1lambda*d1lambda
472
& +(lambda**f52)*d2lambda)
473
tmp11 = -(r42/five)*(f52*(lambda**f32)*d1lambda*d1lambda
474
& +(lambda**f52)*d2lambda)*(dsqrt(zeta)-dsqrt(eta))
475
tmp12 = -(r42/five)*(lambda**f52)*d1lambda*
476
& (d1zeta/dsqrt(zeta)-d1eta/dsqrt(eta))
477
tmp13 = -(six/five)*(lambda**(seven/two))*
478
& (-d1zeta*d1zeta/(two*zeta**f32)+d2zeta/dsqrt(zeta)
479
& +d1eta*d1eta/(two*eta**f32)-d2eta/dsqrt(eta))
481
d2EGs = tmp7+tmp8+tmp9+tmp10+tmp11+tmp12+tmp13
483
c Calculate the second derivate of denominators needed
485
tmp7 = two/(dsqrt(zeta+nu2)*dsqrt(eta+nu2))
486
tmp8 = (dsqrt(eta+nu2)+nu)/(zeta+nu2)**f32
487
tmp9 = (dsqrt(zeta+nu2)+nu)/(eta+nu2)**f32
489
d2Js = f12*(d2zeta*(two+tmp1+tmp2) +
490
& f12*d1zeta*d1zeta*(tmp7-tmp8-tmp9))
494
tmp10 = two/(dsqrt(zeta+nu2)*dsqrt(lambda+nu2))
495
tmp11 = (dsqrt(lambda+nu2)+nu)/(zeta+nu2)**f32
496
tmp12 = (dsqrt(zeta+nu2)+nu)/(lambda+nu2)**f32
497
d2Ms = f12*(d2zeta*(two+tmp3+tmp4) +
498
& f12*d1zeta*d1zeta*(tmp10-tmp11-tmp12))
500
tmp13 = two/(dsqrt(eta+nu2)*dsqrt(lambda+nu2))
501
tmp14 = (dsqrt(lambda+nu2)+nu)/(eta+nu2)**f32
502
tmp15 = (dsqrt(eta+nu2)+nu)/(lambda+nu2)**f32
503
d2Ns = f12*(d2zeta*(two+tmp5+tmp6) +
504
& f12*d1zeta*d1zeta*(tmp13-tmp14-tmp15))
508
c Calculate the second derivative of the Fxhse with respect to the
509
c reduced density gradient, s.
511
tmp1 = A*(Js*Js*d2zeta -two*Js*d1Js*d1zeta
512
& + two*d1Js*d1Js*zeta-Js*d2Js*zeta)/
514
tmp2 = A*(Ks*Ks*d2zeta -two*Ks*d1Ks*d1zeta
515
& + two*d1Ks*d1Ks*eta-Ks*d2Ks*eta)/
517
d2Fxhse1 = tmp1 +tmp2
519
d2Fxhse2 = (four/nine)*B*(L2*d2lambda2-
520
& two*lambda2*d1lambda2*d1lambda2)/(L2*L2)
522
tmp4 = d1lambda2/lambda2
523
tmp5 = d2lambda2/lambda2
524
d2Fxhse3 = -((four*C)/(nine*L2))*(d2Fs +
525
& six*Fs*tmp4**two - two*Fs*tmp5 -
527
& f12*((d2Fs*chi+two*d1Fs*d1chi+Fs*d2chi)-
528
& four*tmp4*(d1Fs*chi+Fs*d1chi) +
529
& six*chi*Fs*tmp4*tmp4)-chi*Fs*tmp5)
531
tmp1 = -d2chi -(two/three)*(d1chi*d1chi+chi*d2chi)
532
& + (six*d1chi+four*chi*d1chi
533
& -tmp4*(r32/three+r12*chi+four*chi*chi))*tmp4
534
& +((eight/three) +three*chi +chi*chi)*tmp5
536
tmp2 = -(two +(four/three)*chi)*d1chi
537
& +(two*chi*chi+(r16/three)+six*chi)*tmp4
539
tmp3 = -((eight/nine)+chi+(one/three)*chi*chi)
541
d2Fxhse4 = (EGs*tmp1 + d1EGs*tmp2+ d2EGs*tmp3)/L3
545
d2Fxhse5 = two*d2zeta*dlog(tmp1)+four*d1zeta*D*d1Ms/(Ms*Ms*tmp1)-
546
& two*zeta*(D*d1Ms/(Ms*Ms*tmp1))**two +
547
& two*zeta*D*(Ms*Ms*d2Ms-two*Ms*d1Ms*d1Ms)/(tmp1*Ms**four)
549
tmp1 = (one-(D-A)/Ns)
550
d2Fxhse6 = -two*d2eta*dlog(tmp1)-
551
& four*d1eta*(D-A)*d1Ns/(Ns*Ns*tmp1)+
552
& two*eta*((D-A)*d1Ns/(Ns*Ns*tmp1))**two -
553
& two*eta*(D-A)*(Ns*Ns*d2Ns-two*Ns*d1Ns*d1Ns)/(tmp1*Ns**four)
555
d20Fxhse = d2Fxhse1+d2Fxhse2+d2Fxhse3+d2Fxhse4+d2Fxhse5+d2Fxhse6
557
c Calculate the second derivative of Fxhse with respect to the
558
c parameter nu (cam_omega/kf).
563
d2Fxhse1 = A*zeta*((two*(one + nu/dsqrt(nu2 + zeta))*
564
& (nu/dsqrt(eta + nu2) + nu/dsqrt(nu2 + zeta)))/
565
& ((nu + dsqrt(nu2 + zeta))**two*(dsqrt(eta + nu2) +
566
& dsqrt(nu2 + zeta))**two) +
567
& ((two*(one + nu/dsqrt(nu2 + zeta))**two)/
568
& (nu + dsqrt(nu2 + zeta))**three -
569
& (-(nu2/(nu2 + zeta)**(three/two)) + one/dsqrt(nu2 + zeta))/
570
& (nu + dsqrt(nu2 + zeta))**two)/
571
& (dsqrt(eta + nu2) + dsqrt(nu2 + zeta)) +
572
& ((two*(nu/dsqrt(eta + nu2) + nu/dsqrt(nu2 + zeta))**two)/
573
& (dsqrt(eta + nu**two) + dsqrt(nu2 + zeta))**three -
574
& (-(nu2/(eta + nu2)**(three/two)) + one/dsqrt(eta + nu2) -
575
& nu2/(nu2 + zeta)**(three/two) + one/dsqrt(nu2 + zeta))/
576
& (dsqrt(eta + nu2) + dsqrt(nu2 + zeta))**two)/
577
& (nu + dsqrt(nu2 + zeta)))
579
d2Fxhse2 = A*eta*((two*(one + nu/dsqrt(eta + nu2))*
580
& (nu/dsqrt(eta + nu2) + nu/dsqrt(nu2 + zeta)))/
581
& ((nu + dsqrt(eta + nu2))**two*(dsqrt(eta + nu2) +
582
& dsqrt(nu2 + zeta))**two) +
583
& ((two*(one + nu/dsqrt(eta + nu2))**two)/
584
& (nu + dsqrt(eta + nu2))**three -
585
& (-(nu2/(eta + nu2)**(three/two)) + one/dsqrt(eta + nu2))/
586
& (nu + dsqrt(eta + nu2))**two)/(dsqrt(eta + nu2) +
587
& dsqrt(nu2 + zeta)) +
588
& ((two*(nu/dsqrt(eta + nu2) + nu/dsqrt(nu2 + zeta))**two)/
589
& (dsqrt(eta + nu2) + dsqrt(nu2 + zeta))**three -
590
& (-(nu2/(eta + nu2)**(three/two)) + one/dsqrt(eta + nu2) -
591
& nu2/(nu2 + zeta)**(three/two) + one/dsqrt(nu2 + zeta))/
592
& (dsqrt(eta + nu2) + dsqrt(nu2 + zeta))**two)/
593
& (nu + dsqrt(eta + nu2)))
595
d2Fxhse3 = (-four*B*((four*nu*(-(nu2/(lambda + nu2)**(three/two))
596
& + one/dsqrt(lambda + nu2)))/((lambda + nu2)**two*
597
& (one + nu/dsqrt(lambda + nu2))**two) +
598
& ((eight*nu2)/(lambda + nu2)**three -
599
& two/(lambda + nu2)**two)/(one + nu/dsqrt(lambda + nu2)) +
600
& ((two*(-(nu2/(lambda + nu2)**(three/two)) +
601
& one/dsqrt(lambda + nu2))**two)/
602
& (one + nu/dsqrt(lambda + nu2))**three -
603
& ((three*nu3)/(lambda + nu2)**(five/two) -
604
& (three*nu)/(lambda + nu2)**(three/two))/
605
& (one + nu/dsqrt(lambda + nu2))**two)/(lambda + nu2)))/nine
607
d2Fxhse4 =(-four*C*Fs*((-four*(-(nu2/(lambda + nu2)**(three/two))
608
& + one/dsqrt(lambda + nu2))*
609
& ((-nu2/(two*(lambda + nu2)**(three/two)) +
610
& one/(two*dsqrt(lambda + nu2)))/(lambda + nu2)**two -
611
& (four*nu*(one + nu/(two*dsqrt(lambda + nu2))))/
612
& (lambda + nu2)**three))/(one + nu/dsqrt(lambda + nu2))**three +
613
& ((one + nu/(two*dsqrt(lambda + nu2)))*
614
& ((six*(-(nu2/(lambda + nu2)**(three/two)) +
615
& one/dsqrt(lambda + nu2))**two)/
616
& (one + nu/dsqrt(lambda + nu2))**four -
617
& (two*((three*nu3)/(lambda + nu2)**(five/two) -
618
& (three*nu)/(lambda + nu2)**(three/two)))/
619
& (one + nu/dsqrt(lambda + nu2))**three))/
620
& (lambda + nu2)**two +
621
& ((-eight*nu*(-nu2/(two*(lambda + nu2)**(three/two)) +
622
& one/(two*dsqrt(lambda + nu2))))/(lambda + nu2)**three +
623
& ((r24*nu2)/(lambda + nu2)**four - four/(lambda + nu2)**three)*
624
& (one + nu/(two*dsqrt(lambda + nu2))) +
625
& (-(nu/(lambda + nu2)**(three/two)) +
626
& (nu*((three*nu2)/(lambda + nu2)**(five/two) -
627
& (lambda + nu2)**(-three/two)))/two)/(lambda + nu2)**two)/
628
& (one + nu/dsqrt(lambda + nu2))**two))/nine
630
d2Fxhse5 = (-eight*EGs*((r48*nu2)/(lambda + nu2)**five -
631
& six/(lambda + nu2)**four)*
632
& (one + (three*nu2)/(eight*(lambda + nu2)) +
633
& (nine*nu)/(eight*dsqrt(lambda + nu2))))/
634
& (nine*(one + nu/dsqrt(lambda + nu2))**three) +
635
& (r32*EGs*nu*(((-three*nu3)/(four*(lambda + nu2)**two) -
636
& (nine*nu2)/(eight*(lambda + nu2)**(three/two)) +
637
& (three*nu)/(four*(lambda + nu2)) +
638
& nine/(eight*dsqrt(lambda + nu2)))/
639
& (one + nu/dsqrt(lambda + nu2))**three -
640
& (three*(-(nu2/(lambda + nu2)**(three/two)) +
641
& one/dsqrt(lambda + nu2))*
642
& (one + (three*nu2)/(eight*(lambda + nu2)) +
643
& (nine*nu)/(eight*dsqrt(lambda + nu2))))/
644
& (one + nu/dsqrt(lambda + nu2))**four))/
645
& (three*(lambda + nu2)**four) -
646
& (eight*EGs*((-six*(-(nu2/(lambda + nu2)**(three/two)) +
647
& one/dsqrt(lambda + nu2))*
648
& ((-three*nu3)/(four*(lambda + nu2)**two) -
649
& (nine*nu2)/(eight*(lambda + nu2)**(three/two)) +
650
& (three*nu)/(four*(lambda + nu2)) +
651
& nine/(eight*dsqrt(lambda + nu2))))/
652
& (one + nu/dsqrt(lambda + nu2))**four +
653
& ((-three*nu2)/(lambda + nu2)**two -
654
& (nine*nu)/(four*(lambda + nu2)**(three/two)) +
655
& three/(four*(lambda + nu2)) +
656
& (three*nu2*((eight*nu2)/(lambda + nu2)**three -
657
& two/(lambda + nu2)**two))/eight +
658
& (nine*nu*((three*nu2)/(lambda + nu2)**(five/two) -
659
& (lambda + nu2)**(-three/two)))/eight)/
660
& (one + nu/dsqrt(lambda + nu2))**three +
661
& (one + (three*nu2)/(eight*(lambda + nu2)) +
662
& (nine*nu)/(eight*dsqrt(lambda + nu2)))*
663
& ((r12*(-(nu2/(lambda + nu2)**(three/two)) +
664
& one/dsqrt(lambda + nu2))**two)/
665
& (one + nu/dsqrt(lambda + nu2))**five -
666
& (three*((three*nu3)/(lambda + nu2)**(five/two) -
667
& (three*nu)/(lambda + nu2)**(three/two)))/
668
& (one + nu/dsqrt(lambda + nu2))**four)))/
669
& (nine*(lambda + nu2)**three)
671
d2Fxhse6 = two*zeta*(-(((D*(nu/dsqrt(lambda + nu2) +
672
& nu/dsqrt(nu2 + zeta)))/
673
& ((nu + dsqrt(lambda + nu2))*(dsqrt(lambda + nu2) +
674
& dsqrt(nu2 + zeta))**two) +
675
& (D*(one + nu/dsqrt(lambda + nu2)))/
676
& ((nu + dsqrt(lambda + nu2))**two*(dsqrt(lambda + nu2) +
677
& dsqrt(nu2 + zeta))))**two/
678
& (one - D/((nu + dsqrt(lambda + nu2))*
679
& (dsqrt(lambda + nu2) + dsqrt(nu2 + zeta))))**two) +
680
& ((-two*D*(nu/dsqrt(lambda + nu2) + nu/dsqrt(nu2 + zeta))**two)/
681
& ((nu + dsqrt(lambda + nu2))*(dsqrt(lambda + nu2) +
682
& dsqrt(nu2 + zeta))**three) +
683
& (D*(-(nu2/(lambda + nu2)**(three/two)) +
684
& one/dsqrt(lambda + nu2) - nu2/(nu2 + zeta)**(three/two) +
685
& one/dsqrt(nu2 + zeta)))/
686
& ((nu + dsqrt(lambda + nu2))*(dsqrt(lambda + nu2) +
687
& dsqrt(nu2 + zeta))**two) -
688
& (two*D*(one + nu/dsqrt(lambda + nu2))*
689
& (nu/dsqrt(lambda + nu2) + nu/dsqrt(nu2 + zeta)))/
690
& ((nu + dsqrt(lambda + nu2))**two*(dsqrt(lambda + nu2) +
691
& dsqrt(nu2 + zeta))**two) -
692
& (two*D*(one + nu/dsqrt(lambda + nu2))**two)/
693
& ((nu + dsqrt(lambda + nu2))**three*(dsqrt(lambda + nu2) +
694
& dsqrt(nu2 + zeta))) +
695
& (D*(-(nu2/(lambda + nu2)**(three/two)) +
696
& one/dsqrt(lambda + nu2)))/
697
& ((nu + dsqrt(lambda + nu2))**two*(dsqrt(lambda + nu2) +
698
& dsqrt(nu2 + zeta))))/
699
& (one - D/((nu + dsqrt(lambda + nu2))*
700
& (dsqrt(lambda + nu2) + dsqrt(nu2 + zeta)))))
702
d2Fxhse7 = -two*eta*(-((((-A + D)*(nu/dsqrt(eta + nu2) +
703
& nu/dsqrt(lambda + nu2)))/
704
& ((nu + dsqrt(lambda + nu2))*(dsqrt(eta + nu2) +
705
& dsqrt(lambda + nu2))**two) +
706
& ((-A + D)*(one + nu/dsqrt(lambda + nu2)))/
707
& ((nu + dsqrt(lambda + nu2))**two*(dsqrt(eta + nu2) +
708
& dsqrt(lambda + nu2))))**two/
709
& (one - (-A + D)/((nu + dsqrt(lambda + nu2))*(dsqrt(eta + nu2) +
710
& dsqrt(lambda + nu2))))**two) +
711
& ((-two*(-A + D)*(nu/dsqrt(eta + nu2) +
712
& nu/dsqrt(lambda + nu2))**two)/
713
& ((nu + dsqrt(lambda + nu2))*(dsqrt(eta + nu2) +
714
& dsqrt(lambda + nu2))**three) -
715
& (two*(-A + D)*(one + nu/dsqrt(lambda + nu2))*
716
& (nu/dsqrt(eta + nu2) + nu/dsqrt(lambda + nu2)))/
717
& ((nu + dsqrt(lambda + nu2))**two*(dsqrt(eta + nu2) +
718
& dsqrt(lambda + nu2))**two) +
719
& ((-A + D)*(-(nu2/(eta + nu2)**(three/two)) +
720
& one/dsqrt(eta + nu2) - nu2/(lambda + nu2)**(three/two) +
721
& one/dsqrt(lambda + nu2)))/
722
& ((nu + dsqrt(lambda + nu2))*(dsqrt(eta + nu2) +
723
& dsqrt(lambda + nu2))**two) -
724
& (two*(-A + D)*(one + nu/dsqrt(lambda + nu2))**two)/
725
& ((nu + dsqrt(lambda + nu2))**three*(dsqrt(eta + nu2) +
726
& dsqrt(lambda + nu2))) +
727
& ((-A + D)*(-(nu2/(lambda + nu2)**(three/two)) +
728
& one/dsqrt(lambda + nu2)))/
729
& ((nu + dsqrt(lambda + nu2))**two*(dsqrt(eta + nu2) +
730
& dsqrt(lambda + nu2))))/
731
& (one - (-A + D)/((nu + dsqrt(lambda + nu2))*
732
& (dsqrt(eta + nu2) + dsqrt(lambda + nu2)))))
734
d02Fxhse = d2Fxhse1+d2Fxhse2+d2Fxhse3+d2Fxhse4+d2Fxhse5+d2Fxhse6+
738
c Calculate the mixed partial derivative of the enhancement factor
739
c with respect to both the reduced density gradient, s, and the
740
c parameter nu (cam_omega/kF). Note that the enhancement factor
741
c satisfies continuity requirements and therefore it does not matter
742
c what order the derivatives are taken in.
746
d11Fxhse1 = (A*(nu*zeta**three*d1eta - two*nu2*(nu2 + eta)*
747
& (eta + two*nu*(nu + dsqrt(nu2 + eta)))*
748
& (nu + dsqrt(nu2 + zeta))*d1zeta +
749
& nu*(eta + two*nu*(nu + dsqrt(nu2 + eta)))*zeta*
750
& (nu*(nu + dsqrt(nu2 + zeta))*d1eta - (nu2 + eta)*d1zeta) +
751
& zeta**two*((eta*(nu + dsqrt(nu2 + zeta)) +
752
& nu*(three*nu2 + two*dsqrt(nu2 + eta)*dsqrt(nu2 + zeta) +
753
& two*nu*(dsqrt(nu2 + eta) + dsqrt(nu2 + zeta))))*
754
& d1eta + (nu2 + eta)*(nu + dsqrt(nu2 + eta))*d1zeta)))/
755
& (two*(nu2 + eta)**(f32)*(nu2 + zeta)**(f32)*
756
& (nu + dsqrt(nu2 + zeta))**two*
757
& (dsqrt(nu2 + eta) + dsqrt(nu2 + zeta))**two)
759
d11Fxhse2 = (A*((-nu2 - zeta)*
760
& (-(eta**two*(nu + dsqrt(nu2 + zeta))) +
761
& nu*eta*(zeta + two*nu*(nu + dsqrt(nu2 + zeta))) +
762
& two*nu2*(nu + dsqrt(nu2 + eta))*(zeta +
763
& two*nu*(nu + dsqrt(nu2 + zeta))))*d1eta +
764
& eta*(nu2 + eta)*(nu*eta + (nu + dsqrt(nu2 + eta))*
765
& (zeta + two*nu*(nu + dsqrt(nu2 + zeta))))*d1zeta))/
766
& (two*(nu2 + eta)**(f32)*(nu + dsqrt(nu2 + eta))**two*
767
& (nu2 + zeta)**(f32)*
768
& (dsqrt(nu2 + eta) + dsqrt(nu2 + zeta))**two)
770
d11Fxhse3 = (-two*B*(two*nu2*(nu + dsqrt(nu2 + lambda)) +
771
& lambda*(two*nu + dsqrt(nu2 + lambda)))*d1lambda)/
772
& (three*(nu2 + lambda)**three*
773
& (nu + dsqrt(nu2 + lambda))**two)
775
d11Fxhse4 = (C*(two*(nu2 + lambda)*d1Fs - five*Fs*d1lambda))/
776
& (three*(nu2 + lambda)**(f72))
778
d11Fxhse5 = (five*(eight*nu4*(nu + dsqrt(nu2 + lambda)) +
779
& lambda**two*(four*nu + dsqrt(nu2 + lambda)) +
780
& four*nu2*lambda*(three*nu + two*dsqrt(nu2 + lambda)))*
781
& (two*(nu2 + lambda)*d1EGs - seven*EGs*d1lambda))/
782
& (six*(nu2 + lambda)**five*(nu + dsqrt(nu2 + lambda))**four)
784
d11Fxhse6 = (D*((-nu - two*dsqrt(nu2 + lambda))*zeta**three*
786
& two*nu2*(nu2 + lambda)*(-D + lambda + two*nu*
787
& (nu + dsqrt(nu2 + lambda)))*(nu + dsqrt(nu2 + zeta))*
788
& d1zeta + zeta**two*((D*(nu + dsqrt(nu2 + zeta)) -
789
& three*lambda*(nu + dsqrt(nu2 + zeta)) -
790
& nu*(five*nu2 + six*nu*dsqrt(nu2 + lambda) +
791
& four*nu*dsqrt(nu2 + zeta) + four*dsqrt(nu2 + lambda)*
792
& dsqrt(nu2 + zeta)))*
793
& d1lambda + (nu2 + lambda)*(nu + dsqrt(nu2 + lambda))*d1zeta) +
794
& zeta*(-(nu2*(-D + three*lambda + four*nu*
795
& (nu + dsqrt(nu2 + lambda)))*(nu + dsqrt(nu2 + zeta))*d1lambda) +
796
& (nu2 + lambda)*(two*nu*(nu + dsqrt(nu2 + lambda))*
797
& (two*nu + dsqrt(nu2 + zeta)) -
798
& D*(nu + two*dsqrt(nu2 + zeta)) +
799
& lambda*(nu + two*dsqrt(nu2 + zeta)))*d1zeta)))/
800
& ((nu2 + lambda)**(f32)*(nu2 + zeta)**(f32)*
801
& (-D + lambda + (nu + dsqrt(nu2 + lambda))*
802
& (nu + dsqrt(nu2 + zeta)))**two)
804
d11Fxhse7 = ((A - D)*((nu2 + lambda)*(eta**two*
805
& (nu + dsqrt(nu2 + lambda)) +
806
& two*nu2*(nu + dsqrt(nu2 + eta))*
807
& (A - D + two*nu2 + lambda + two*nu*dsqrt(nu2 + lambda)) +
808
& eta*(A*nu - D*nu + four*nu3 + two*A*dsqrt(nu2 + eta) -
809
& two*D*dsqrt(nu2 + eta) + two*nu2*dsqrt(nu2 + eta) +
810
& (nu + two*dsqrt(nu2 + eta))*lambda +
811
& four*nu2*dsqrt(nu2 + lambda) +
812
& two*nu*dsqrt(nu2 + eta)*dsqrt(nu2 + lambda)))*d1eta -
814
& (eta*(nu + two*dsqrt(nu2 + lambda)) +
815
& (nu + dsqrt(nu2 + eta))*
816
& (A - D + four*nu2 + three*lambda +
817
& four*nu*dsqrt(nu2 + lambda)))*d1lambda))/
818
& ((nu2 + eta)**(f32)*(nu2 + lambda)**(f32)*
819
& (A - D + nu2 + nu*dsqrt(nu2 + eta) +
820
& lambda + nu*dsqrt(nu2 + lambda) +
821
& dsqrt(nu2 + eta)*dsqrt(nu2 + lambda))**two)
824
c d11Fxhse1 = (A*zeta*(nu/dsqrt(nu2 + eta) +
825
c & nu/dsqrt(nu2 + zeta))*d1zeta)/
826
c & (two*dsqrt(nu2 + zeta)*(nu + dsqrt(nu2 + zeta))**two*
827
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + zeta))**two) -
828
c & (A*(nu/dsqrt(nu2 + eta) + nu/dsqrt(nu2 + zeta))*d1zeta)/
829
c & ((nu + dsqrt(nu2 + zeta))*(dsqrt(nu2 + eta) +
830
c & dsqrt(nu2 + zeta))**two) +
831
c & (A*zeta*(one + nu/dsqrt(nu2 + zeta))*d1zeta)/
832
c & (dsqrt(nu2 + zeta)*(nu + dsqrt(nu2 + zeta))**three*
833
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + zeta))) +
834
c & (A*nu*zeta*d1zeta)/
835
c & (two*(nu2 + zeta)**(three/two)*(nu + dsqrt(nu2 + zeta))**two*
836
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + zeta))) -
837
c & (A*(one + nu/dsqrt(nu2 + zeta))*d1zeta)/
838
c & ((nu + dsqrt(nu2 + zeta))**two*(dsqrt(nu2 + eta) +
839
c & dsqrt(nu2 + zeta))) -
840
c & (A*zeta*(-(nu*d1eta)/(two*(nu2 + eta)**(three/two)) -
841
c & (nu*d1zeta)/(two*(nu2 + zeta)**(three/two))))/
842
c & ((nu + dsqrt(nu2 + zeta))*(dsqrt(nu2 + eta) +
843
c & dsqrt(nu2 + zeta))**two) +
844
c & (two*A*zeta*(nu/dsqrt(nu2 + eta) + nu/dsqrt(nu2 + zeta))*
845
c & (d1eta/(two*dsqrt(nu2 + eta)) +
846
c & d1zeta/(two*dsqrt(nu2 + zeta))))/
847
c & ((nu + dsqrt(nu2 + zeta))*(dsqrt(nu2 + eta) +
848
c & dsqrt(nu2 + zeta))**three) +
849
c & (A*zeta*(one + nu/dsqrt(nu2 + zeta))*
850
c & (d1eta/(two*dsqrt(nu2 + eta)) +
851
c & d1zeta/(two*dsqrt(nu2 + zeta))))/
852
c & ((nu + dsqrt(nu2 + zeta))**two*(dsqrt(nu2 + eta) +
853
c & dsqrt(nu2 + zeta))**two)
855
c d11Fxhse2 = (A*eta*(nu/dsqrt(nu2 + eta) +
856
c & nu/dsqrt(nu2 + zeta))*d1eta)/
857
c & (two*dsqrt(nu2 + eta)*(nu + dsqrt(nu2 + eta))**two*
858
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + zeta))**two) -
859
c & (A*(nu/dsqrt(nu2 + eta) + nu/dsqrt(nu2 + zeta))*d1eta)/
860
c & ((nu + dsqrt(nu2 + eta))*(dsqrt(nu2 + eta) +
861
c & dsqrt(nu2 + zeta))**two) +
862
c & (A*eta*(one + nu/dsqrt(nu2 + eta))*d1eta)/
863
c & (dsqrt(nu2 + eta)*(nu + dsqrt(nu2 + eta))**three*
864
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + zeta))) +
865
c & (A*nu*eta*d1eta)/
866
c & (two*(nu2 + eta)**(three/two)*(nu + dsqrt(nu2 + eta))**two*
867
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + zeta))) -
868
c & (A*(one + nu/dsqrt(nu2 + eta))*d1eta)/
869
c & ((nu + dsqrt(nu2 + eta))**two*(dsqrt(nu2 + eta) +
870
c & dsqrt(nu2 + zeta))) -
871
c & (A*eta*(-(nu*d1eta)/(two*(nu2 + eta)**(three/two)) -
872
c & (nu*d1zeta)/(two*(nu2 + zeta)**(three/two))))/
873
c & ((nu + dsqrt(nu2 + eta))*(dsqrt(nu2 + eta) +
874
c & dsqrt(nu2 + zeta))**two) +
875
c & (two*A*eta*(nu/dsqrt(nu2 + eta) + nu/dsqrt(nu2 + zeta))*
876
c & (d1eta/(two*dsqrt(nu2 + eta)) +
877
c & d1zeta/(two*dsqrt(nu2 + zeta))))/
878
c & ((nu + dsqrt(nu2 + eta))*(dsqrt(nu2 + eta) +
879
c & dsqrt(nu2 + zeta))**three) +
880
c & (A*eta*(one + nu/dsqrt(nu2 + eta))*
881
c & (d1eta/(two*dsqrt(nu2 + eta)) +
882
c & d1zeta/(two*dsqrt(nu2 + zeta))))/
883
c & ((nu + dsqrt(nu2 + eta))**two*(dsqrt(nu2 + eta) +
884
c & dsqrt(nu2 + zeta))**two)
886
c d11Fxhse3 = (four*B*nu*(-(nu2/(nu2 + lambda)**(three/two)) +
887
c & one/dsqrt(nu2 + lambda))*d1lambda)/
888
c & (nine*(nu2 + lambda)**(five/two)*
889
c & (one + nu/dsqrt(nu2 + lambda))**three) +
890
c & (four*B*nu2*d1lambda)/(nine*(nu2 + lambda)**(seven/two)*
891
c & (one + nu/dsqrt(nu2 + lambda))**two) -
892
c & (four*B*(-(nu2/(nu2 + lambda)**(three/two)) +
893
c & one/dsqrt(nu2 + lambda))*d1lambda)/
894
c & (nine*(nu2 + lambda)**two*
895
c & (one + nu/dsqrt(nu2 + lambda))**two) -
896
c & (r16*B*nu*d1lambda)/(nine*(nu2 + lambda)**three*
897
c & (one + nu/dsqrt(nu2 + lambda))) +
898
c & (four*B*((three*nu2*d1lambda)/
899
c & (two*(nu2 + lambda)**(five/two)) -
900
c & d1lambda/(two*(nu2 + lambda)**(three/two))))/
901
c & (nine*(nu2 + lambda)*(one + nu/dsqrt(nu2 + lambda))**two)
903
c d11Fxhse4 = (eight*C*(-(nu2/(nu2 + lambda)**(three/two)) +
904
c & one/dsqrt(nu2 + lambda))*
905
c & (one + nu/(two*dsqrt(nu2 + lambda)))*d1Fs)/
906
c & (nine*(nu2 + lambda)**two*
907
c & (one + nu/dsqrt(nu2 + lambda))**three) -
908
c & (four*C*(-nu2/(two*(nu2 + lambda)**(three/two)) +
909
c & one/(two*dsqrt(nu2 + lambda)))*d1Fs)/
910
c & (nine*(nu2 + lambda)**two*
911
c & (one + nu/dsqrt(nu2 + lambda))**two) +
912
c & (r16*C*nu*(one + nu/(two*dsqrt(nu2 + lambda)))*d1Fs)/
913
c & (nine*(nu2 + lambda)**three*
914
c & (one + nu/dsqrt(nu2 + lambda))**two) +
915
c & (four*C*nu*Fs*(-(nu2/(nu2 + lambda)**(three/two)) +
916
c & one/dsqrt(nu2 + lambda))*(one + nu/
917
c & (two*dsqrt(nu2 + lambda)))*
918
c & d1lambda)/(three*(nu2 + lambda)**(seven/two)*
919
c & (one + nu/dsqrt(nu2 + lambda))**four) -
920
c & (four*C*nu*Fs*(-nu2/(two*(nu2 + lambda)**(three/two)) +
921
c & one/(two*dsqrt(nu2 + lambda)))*d1lambda)/
922
c & (nine*(nu2 + lambda)**(seven/two)*
923
c & (one + nu/dsqrt(nu2 + lambda))**three) -
924
c & (two*C*nu*Fs*(-(nu2/(nu2 + lambda)**(three/two)) +
925
c & one/dsqrt(nu2 + lambda))*d1lambda)/
926
c & (nine*(nu2 + lambda)**(seven/two)*
927
c & (one + nu/dsqrt(nu2 + lambda))**three) +
928
c & (r16*C*nu2*Fs*(one + nu/(two*dsqrt(nu2 + lambda)))*d1lambda)/
929
c & (nine*(nu2 + lambda)**(nine/two)*
930
c & (one + nu/dsqrt(nu2 + lambda))**three) -
931
c & (r16*C*Fs*(-(nu2/(nu2 + lambda)**(three/two)) +
932
c & one/dsqrt(nu2 + lambda))*(one + nu/(two*dsqrt(nu2 + lambda)))*
933
c & d1lambda)/(nine*(nu2 + lambda)**three*
934
c & (one + nu/dsqrt(nu2 + lambda))**three) -
935
c & (four*C*nu2*Fs*d1lambda)/(nine*(nu2 + lambda)**(nine/two)*
936
c & (one + nu/dsqrt(nu2 + lambda))**two) +
937
c & (eight*C*Fs*(-nu2/(two*(nu2 + lambda)**(three/two)) +
938
c & one/(two*dsqrt(nu2 + lambda)))*d1lambda)/
939
c & (nine*(nu2 + lambda)**three*
940
c & (one + nu/dsqrt(nu2 + lambda))**two) -
941
c & (r16*C*nu*Fs*(one + nu/(two*dsqrt(nu2 + lambda)))*d1lambda)/
942
c & (three*(nu2 + lambda)**four*
943
c & (one + nu/dsqrt(nu2 + lambda))**two) +
944
c & (eight*C*Fs*(one + nu/(two*dsqrt(nu2 + lambda)))*
945
c & ((three*nu2*d1lambda)/(two*(nu2 + lambda)**(five/two)) -
946
c & d1lambda/(two*(nu2 + lambda)**(three/two))))/
947
c & (nine*(nu2 + lambda)**two*
948
c & (one + nu/dsqrt(nu2 + lambda))**three) -
949
c & (four*C*Fs*((three*nu2*d1lambda)/
950
c & (four*(nu2 + lambda)**(five/two)) -
951
c & d1lambda/(four*(nu2 + lambda)**(three/two))))/
952
c & (nine*(nu2 + lambda)**two*
953
c & (one + nu/dsqrt(nu2 + lambda))**two)
955
c d11Fxhse5 = (-eight*((-three*nu3)/(four*(nu2 + lambda)**two) -
956
c & (nine*nu2)/(eight*(nu2 + lambda)**(three/two)) +
957
c & (three*nu)/(four*(nu2 + lambda)) +
958
c & nine/(eight*dsqrt(nu2 + lambda)))*d1EGs)/
959
c & (nine*(nu2 + lambda)**three*
960
c & (one + nu/dsqrt(nu2 + lambda))**three) +
961
c & (eight*(-(nu2/(nu2 + lambda)**(three/two)) +
962
c & one/dsqrt(nu2 + lambda))*
963
c & (one + (three*nu2)/(eight*(nu2 + lambda)) +
964
c & (nine*nu)/(eight*dsqrt(nu2 + lambda)))*d1EGs)/
965
c & (three*(nu2 + lambda)**three*
966
c & (one + nu/dsqrt(nu2 + lambda))**four) +
967
c & (r16*nu*(one + (three*nu2)/(eight*(nu2 + lambda)) +
968
c & (nine*nu)/(eight*dsqrt(nu2 + lambda)))*d1EGs)/
969
c & (three*(nu2 + lambda)**four*
970
c & (one + nu/dsqrt(nu2 + lambda))**three) -
971
c & (four*nu*EGs*((-three*nu3)/(four*(nu2 + lambda)**two) -
972
c & (nine*nu2)/(eight*(nu2 + lambda)**(three/two)) +
973
c & (three*nu)/(four*(nu2 + lambda)) +
974
c & nine/(eight*dsqrt(nu2 + lambda)))*d1lambda)/
975
c & (three*(nu2 + lambda)**(nine/two)*
976
c & (one + nu/dsqrt(nu2 + lambda))**four) +
977
c & (eight*EGs*((-three*nu3)/(four*(nu2 + lambda)**two) -
978
c & (nine*nu2)/(eight*(nu2 + lambda)**(three/two)) +
979
c & (three*nu)/(four*(nu2 + lambda)) +
980
c & nine/(eight*dsqrt(nu2 + lambda)))*d1lambda)/
981
c & (three*(nu2 + lambda)**four*
982
c & (one + nu/dsqrt(nu2 + lambda))**three) +
983
c & (r16*nu*EGs*(-(nu2/(nu2 + lambda)**(three/two)) +
984
c & one/dsqrt(nu2 + lambda))*
985
c & (one + (three*nu2)/(eight*(nu2 + lambda)) +
986
c & (nine*nu)/(eight*dsqrt(nu2 + lambda)))*d1lambda)/
987
c & (three*(nu2 + lambda)**(nine/two)*
988
c & (one + nu/dsqrt(nu2 + lambda))**five) +
989
c & (eight*nu2*EGs*(one + (three*nu2)/(eight*(nu2 + lambda)) +
990
c & (nine*nu)/(eight*dsqrt(nu2 + lambda)))*d1lambda)/
991
c & ((nu2 + lambda)**(r11/two)*
992
c & (one + nu/dsqrt(nu2 + lambda))**four) -
993
c & (eight*EGs*(-(nu2/(nu2 + lambda)**(three/two)) +
994
c & one/dsqrt(nu2 + lambda))*
995
c & (one + (three*nu2)/(eight*(nu2 + lambda)) +
996
c & (nine*nu)/(eight*dsqrt(nu2 + lambda)))*d1lambda)/
997
c & ((nu2 + lambda)**four*(one + nu/dsqrt(nu2 + lambda))**four) -
998
c & (r64*nu*EGs*(one + (three*nu2)/(eight*(nu2 + lambda)) +
999
c & (nine*nu)/(eight*dsqrt(nu2 + lambda)))*d1lambda)/
1000
c & (three*(nu2 + lambda)**five*
1001
c & (one + nu/dsqrt(nu2 + lambda))**three) -
1002
c & (eight*EGs*((three*nu3*d1lambda)/
1003
c & (two*(nu2 + lambda)**three) +
1004
c & (r27*nu2*d1lambda)/(r16*(nu2 + lambda)**(five/two)) -
1005
c & (three*nu*d1lambda)/(four*(nu2 + lambda)**two) -
1006
c & (nine*d1lambda)/(r16*(nu2 + lambda)**(three/two))))/
1007
c & (nine*(nu2 + lambda)**three*
1008
c & (one + nu/dsqrt(nu2 + lambda))**three) +
1009
c & (eight*EGs*(one + (three*nu2)/(eight*(nu2 + lambda)) +
1010
c & (nine*nu)/(eight*dsqrt(nu2 + lambda)))*
1011
c & ((three*nu2*d1lambda)/(two*(nu2 + lambda)**(five/two)) -
1012
c & d1lambda/(two*(nu2 + lambda)**(three/two))))/
1013
c & (three*(nu2 + lambda)**three*
1014
c & (one + nu/dsqrt(nu2 + lambda))**four) +
1015
c & (r16*nu*EGs*((-three*nu2*d1lambda)/
1016
c & (eight*(nu2 + lambda)**two) -
1017
c & (nine*nu*d1lambda)/(r16*(nu2 + lambda)**(three/two))))/
1018
c & (three*(nu2 + lambda)**four*
1019
c & (one + nu/dsqrt(nu2 + lambda))**three)
1021
c d11Fxhse6 = (two*((D*(nu/dsqrt(nu2 + lambda) +
1022
c & nu/dsqrt(nu2 + zeta)))/
1023
c & ((nu + dsqrt(nu2 + lambda))*(dsqrt(nu2 + lambda) +
1024
c & dsqrt(nu2 + zeta))**two) +
1025
c & (D*(one + nu/dsqrt(nu2 + lambda)))/
1026
c & ((nu + dsqrt(nu2 + lambda))**two*(dsqrt(nu2 + lambda) +
1027
c & dsqrt(nu2 + zeta))))*
1028
c & d1zeta)/(one - D/((nu + dsqrt(nu2 + lambda))*
1029
c & (dsqrt(nu2 + lambda) + dsqrt(nu2 + zeta)))) +
1030
c & (two*zeta*(-(D*(nu/dsqrt(nu2 + lambda) +
1031
c & nu/dsqrt(nu2 + zeta))*d1lambda)/
1032
c & (two*dsqrt(nu2 + lambda)*(nu + dsqrt(nu2 + lambda))**two*
1033
c & (dsqrt(nu2 + lambda) + dsqrt(nu2 + zeta))**two) -
1034
c & (D*(one + nu/dsqrt(nu2 + lambda))*d1lambda)/
1035
c & (dsqrt(nu2 + lambda)*(nu + dsqrt(nu2 + lambda))**three*
1036
c & (dsqrt(nu2 + lambda) + dsqrt(nu2 + zeta))) -
1037
c & (D*nu*d1lambda)/
1038
c & (two*(nu2 + lambda)**(three/two)*
1039
c & (nu + dsqrt(nu2 + lambda))**two*(dsqrt(nu2 + lambda) +
1040
c & dsqrt(nu2 + zeta))) +
1041
c & (D*(-(nu*d1lambda)/(two*(nu2 + lambda)**(three/two)) -
1042
c & (nu*d1zeta)/(two*(nu2 + zeta)**(three/two))))/
1043
c & ((nu + dsqrt(nu2 + lambda))*(dsqrt(nu2 + lambda) +
1044
c & dsqrt(nu2 + zeta))**two) -
1045
c & (two*D*(nu/dsqrt(nu2 + lambda) + nu/dsqrt(nu2 + zeta))*
1046
c & (d1lambda/(two*dsqrt(nu2 + lambda)) +
1047
c & d1zeta/(two*dsqrt(nu2 + zeta))))/
1048
c & ((nu + dsqrt(nu2 + lambda))*(dsqrt(nu2 + lambda) +
1049
c & dsqrt(nu2 + zeta))**three) -
1050
c & (D*(one + nu/dsqrt(nu2 + lambda))*
1051
c & (d1lambda/(two*dsqrt(nu2 + lambda)) +
1052
c & d1zeta/(two*dsqrt(nu2 + zeta))))/
1053
c & ((nu + dsqrt(nu2 + lambda))**two*(dsqrt(nu2 + lambda) +
1054
c & dsqrt(nu2 + zeta))**two)))/
1055
c & (one - D/((nu + dsqrt(nu2 + lambda))*
1056
c & (dsqrt(nu2 + lambda) + dsqrt(nu2 + zeta)))) -
1057
c & (two*zeta*((D*(nu/dsqrt(nu2 + lambda) +
1058
c & nu/dsqrt(nu2 + zeta)))/
1059
c & ((nu + dsqrt(nu2 + lambda))*(dsqrt(nu2 + lambda) +
1060
c & dsqrt(nu2 + zeta))**two) +
1061
c & (D*(one + nu/dsqrt(nu2 + lambda)))/
1062
c & ((nu + dsqrt(nu2 + lambda))**two*(dsqrt(nu2 + lambda) +
1063
c & dsqrt(nu2 + zeta))))*
1065
c & (two*dsqrt(nu2 + lambda)*(nu + dsqrt(nu2 + lambda))**two*
1066
c & (dsqrt(nu2 + lambda) + dsqrt(nu2 + zeta))) +
1067
c & (D*(d1lambda/(two*dsqrt(nu2 + lambda)) +
1068
c & d1zeta/(two*dsqrt(nu2 + zeta))))/
1069
c & ((nu + dsqrt(nu2 + lambda))*(dsqrt(nu2 + lambda) +
1070
c & dsqrt(nu2 + zeta))**two)))/
1071
c & (one - D/((nu + dsqrt(nu2 + lambda))*
1072
c & (dsqrt(nu2 + lambda) + dsqrt(nu2 + zeta))))**two
1074
c d11Fxhse7 = (-two*(((-A + D)*(nu/dsqrt(nu2 + eta) +
1075
c & nu/dsqrt(nu2 + lambda)))/
1076
c & ((nu + dsqrt(nu2 + lambda))*
1077
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + lambda))**two) +
1078
c & ((-A + D)*(one + nu/dsqrt(nu2 + lambda)))/
1079
c & ((nu + dsqrt(nu2 + lambda))**two*(dsqrt(nu2 + eta) +
1080
c & dsqrt(nu2 + lambda))))*d1eta)/(one - (-A + D)/
1081
c & ((nu + dsqrt(nu2 + lambda))*
1082
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + lambda)))) -
1083
c & (two*eta*(-((-A + D)*(nu/dsqrt(nu2 + eta) +
1084
c & nu/dsqrt(nu2 + lambda))*d1lambda)/
1085
c & (two*dsqrt(nu2 + lambda)*(nu + dsqrt(nu2 + lambda))**two*
1086
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + lambda))**two) -
1087
c & ((-A + D)*(one + nu/dsqrt(nu2 + lambda))*d1lambda)/
1088
c & (dsqrt(nu2 + lambda)*(nu + dsqrt(nu2 + lambda))**three*
1089
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + lambda))) -
1090
c & ((-A + D)*nu*d1lambda)/
1091
c & (two*(nu2 + lambda)**(three/two)*
1092
c & (nu + dsqrt(nu2 + lambda))**two*
1093
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + lambda))) +
1094
c & ((-A + D)*(-(nu*d1eta)/(two*(nu2 + eta)**(three/two)) -
1095
c & (nu*d1lambda)/(two*(nu2 + lambda)**(three/two))))/
1096
c & ((nu + dsqrt(nu2 + lambda))*
1097
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + lambda))**two) -
1098
c & (two*(-A + D)*(nu/dsqrt(nu2 + eta) +
1099
c & nu/dsqrt(nu2 + lambda))*
1100
c & (d1eta/(two*dsqrt(nu2 + eta)) +
1101
c & d1lambda/(two*dsqrt(nu2 + lambda))))/
1102
c & ((nu + dsqrt(nu2 + lambda))*
1103
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + lambda))**three) -
1104
c & ((-A + D)*(one + nu/dsqrt(nu2 + lambda))*
1105
c & (d1eta/(two*dsqrt(nu2 + eta)) +
1106
c & d1lambda/(two*dsqrt(nu2 + lambda))))/
1107
c & ((nu + dsqrt(nu2 + lambda))**two*(dsqrt(nu2 + eta) +
1108
c & dsqrt(nu2 + lambda))**two)))/
1109
c & (one - (-A + D)/((nu + dsqrt(nu2 + lambda))*
1110
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + lambda)))) +
1111
c & (two*eta*(((-A + D)*(nu/dsqrt(nu2 + eta) +
1112
c & nu/dsqrt(nu2 + lambda)))/
1113
c & ((nu + dsqrt(nu2 + lambda))*(dsqrt(nu2 + eta) +
1114
c & dsqrt(nu2 + lambda))**two) +
1115
c & ((-A + D)*(one + nu/dsqrt(nu2 + lambda)))/
1116
c & ((nu + dsqrt(nu2 + lambda))**two*
1117
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + lambda))))*
1118
c & (((-A + D)*d1lambda)/
1119
c & (two*dsqrt(nu2 + lambda)*(nu + dsqrt(nu2 + lambda))**two*
1120
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + lambda))) +
1121
c & ((-A + D)*(d1eta/(two*dsqrt(nu2 + eta)) +
1122
c & d1lambda/(two*dsqrt(nu2 + lambda))))/
1123
c & ((nu + dsqrt(nu2 + lambda))*(dsqrt(nu2 + eta) +
1124
c & dsqrt(nu2 + lambda))**two)))/
1125
c & (one - (-A + D)/((nu + dsqrt(nu2 + lambda))*
1126
c & (dsqrt(nu2 + eta) + dsqrt(nu2 + lambda))))**two
1129
d11Fxhse = d11Fxhse1+d11Fxhse2+d11Fxhse3+d11Fxhse4+d11Fxhse5+
1130
& d11Fxhse6+d11Fxhse7
1136
#ifndef SECOND_DERIV
1137
#define SECOND_DERIV
1139
c Compile source again for the 2nd derivative case
1141
#include "nwxc_x_hse08.F"
1143
c $Id: nwxc_x_hse08.F 23648 2013-02-27 21:52:25Z d3y133 $
added
removed
src/nwxc/nwxc_x_hse08.F
