1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2 Subroutine HSE08Fx(ipol,rho,s,Fxhse,d10Fxhse,d01Fxhse) 3#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 4 Subroutine HSE08Fx_d2(ipol,rho,s,Fxhse,d10Fxhse,d01Fxhse, 5 & d20Fxhse,d02Fxhse,d11Fxhse) 6#else 7 Subroutine HSE08Fx_d3(ipol,rho,s,Fxhse,d10Fxhse,d01Fxhse, 8 & d20Fxhse,d02Fxhse,d11Fxhse,d30Fxhse, 9 & d21Fxhse,d12Fxhse,d03Fxhse) 10#endif 11c 12c$Id$ 13c 14 implicit none 15c 16ccase...start 17#include "case.fh" 18ccase...end 19c 20c HSE evaluates the Heyd et al. Screened Coulomb 21c Exchange Functional 22c 23c Calculates the enhancement factor 24c 25 integer ipol 26 double precision rho,s,Fxhse,d10Fxhse,d01Fxhse 27c 28#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 29 double precision d20Fxhse,d02Fxhse,d11Fxhse 30#endif 31c 32#ifdef THIRD_DERIV 33 double precision d30Fxhse,d03Fxhse,d21Fxhse,d12Fxhse 34#endif 35 36 double precision A,B,C,D,E 37 double precision ha2,ha3,ha4,ha5,ha6,ha7 38 double precision hb1,hb2,hb3,hb4,hb5,hb6,hb7,hb8,hb9 39 double precision smax,strans,sconst 40c 41 double precision zero,one,two,three,four,five,six,seven,eight 42 double precision nine,ten 43 double precision fifteen,sixteen 44 45 double precision H,hnum,hden 46 double precision d1H,d1hnum,d1hden 47 double precision s2,s3,s4,s5,s6,s7,s8,s9 48 double precision Fs, d1Fs 49 double precision zeta, lambda, eta, kf, nu, chi, lambda2 50 double precision d1zeta,d1lambda,d1eta,d1nu,d10chi,d10lambda2 51 double precision d01chi,d01lambda2 52 double precision EGs,d1EGs 53 double precision nu2,L2,L3,nu3,nu4,nu5,nu6,nu7,nu8 54 double precision Js,Ks,Ms,Ns 55 double precision Js2,Ks2,Ms2,Ns2 56 double precision d10Js,d10Ks,d10Ms,d10Ns 57 double precision d01Js,d01Ks,d01Ms,d01Ns 58 59 double precision tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8 60 double precision tmp9,tmp10,tmp11,tmp12,tmp13,tmp14,tmp15 61 double precision Fxhse1,Fxhse2,Fxhse3,Fxhse4,Fxhse5,Fxhse6 62 double precision d1Fxhse1,d1Fxhse2,d1Fxhse3,d1Fxhse4,d1Fxhse5 63 double precision d1Fxhse6,d1Fxhse7 64 65 double precision r42,r27,r12,r15,r14,r18,r20,r30,r56,r72 66 double precision r16,r32,r24,r48,r11,r64,r35,r60,r36,r63 67 double precision r189,r256,r120,r210,r336,r504,r21,r105,r126 68 double precision pi,pi2,srpi,s02 69 double precision f12,f13,f32,f52,f72,f92 70#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 71 double precision d2H,d2hnum,d2hden 72 double precision d2Fs 73 double precision d2zeta,d2lambda,d2eta,d2nu,d20chi,d20lambda2 74 double precision d02chi,d02lambda2,d11chi,d11lambda2 75 double precision d2EGs 76 double precision L5 77 double precision Js3,Ks3,Ms4,Ns4 78 double precision d20Js,d20Ks,d20Ms,d20Ns 79 double precision d02Js,d02Ks,d02Ms,d02Ns 80 double precision d11Js,d11Ks,d11Ms,d11Ns 81 double precision d2Fxhse1,d2Fxhse2,d2Fxhse3,d2Fxhse4,d2Fxhse5 82 double precision d2Fxhse6,d2Fxhse7 83 double precision d11Fxhse1,d11Fxhse2,d11Fxhse3,d11Fxhse4 84 double precision d11Fxhse5,d11Fxhse6,d11Fxhse7 85#endif 86c 87#ifdef THIRD_DERIV 88 double precision d3H,d3hnum,d3hden 89 double precision d3Fs 90 double precision d3zeta,d3lambda,d3eta,d3nu,d30chi,d30lambda2 91 double precision d03chi,d03lambda2,d21chi,d12chi 92 double precision d21lambda2,d12lambda2 93 double precision L4,L6 94 double precision Js4,Ks4,Ms3,Ns3 95 double precision d3EGs 96 double precision d30Js,d30Ks,d30Ms,d30Ns 97 double precision d03Js,d03Ks,d03Ms,d03Ns 98 double precision d21Js,d21Ks,d21Ms,d21Ns 99 double precision d12Js,d12Ks,d12Ms,d12Ns 100 double precision d3Fxhse1,d3Fxhse2,d3Fxhse3,d3Fxhse4,d3Fxhse5 101 double precision d3Fxhse6,d3Fxhse7 102 double precision d21Fxhse1,d21Fxhse2,d21Fxhse3,d21Fxhse4 103 double precision d21Fxhse5,d21Fxhse6,d21Fxhse7 104 double precision d12Fxhse1,d12Fxhse2,d12Fxhse3,d12Fxhse4 105 double precision d12Fxhse5,d12Fxhse6,d12Fxhse7 106#endif 107c 108c Constants for HJS hole 109c 110 Data A,B,C,D,E 111 1 / 7.57211D-1,-1.06364D-1,-1.18649D-1, 112 2 6.09650D-1,-4.77963D-2 / 113c 114c Constants for fit of H(s) (PBE hole) 115c Taken from JCTC_5_754 (2009) 116c 117 Data ha2,ha3,ha4,ha5,ha6,ha7 118 1 / 1.59941D-2,8.52995D-2,-1.60368D-1,1.52645D-1, 119 2 -9.71263D-2,4.22061D-2 / 120c 121 Data hb1,hb2,hb3,hb4,hb5,hb6,hb7,hb8,hb9 122 1 / 5.33319D0,-12.4780D0,11.0988D0,-5.11013D0, 123 2 1.71468D0,-6.10380D-1,3.07555D-1,-7.70547D-2, 124 3 3.34840D-2 / 125 126c 127c Whole numbers used during evaluation 128c 129 Data zero,one,two,three,four,five,six,seven,eight,nine,ten 130 1 / 0.D0,1.D0,2.D0,3.D0,4.D0,5.D0,6.D0,7.D0,8.D0,9.D0,10.D0 / 131 132 Data r11,r12,r14,r15,r16,r18,r20,r24,r27,r30,r32 133 1 / 11.D0,12.D0,14.D0,15.D0,16.D0,18.D0,20.D0,24.D0,27.d0, 134 2 30.D0,32.D0 / 135 136 Data r35,r42,r48,r56,r64,r72 137 1 / 35.D0,42.D0,48.D0,56.D0,64.D0,72.D0 / 138 139 Data r21,r36,r60,r63,r120,r189,r210,r256,r336,r504 140 1 / 21.D0,36.D0,60.D0,63.D0,120.D0,189.D0,210.D0,256.D0, 141 2 336.D0,504.D0 / 142 143 Data r105,r126 144 1 / 105.D0,126.D0 / 145c 146c Fractions used during evaluation 147c 148 Data f12 149 1 / 0.5D0 / 150c 151c General constants 152c 153 f13 = one/three 154 f32 = three/two 155 f52 = five/two 156 f72 = seven/two 157 f92 = nine/two 158 pi = ACos(-one) 159 pi2 = pi*pi 160 srpi = dsqrt(pi) 161c 162c 163c Calculate prelim variables 164c 165 s2 = s*s 166 s02 = s2/four 167 s3 = s2*s 168 s4 = s3*s 169 s5 = s4*s 170 s6 = s5*s 171 s7 = s6*s 172 s8 = s7*s 173 s9 = s8*s 174 175!********************************* 176! Calculate the Enhancement Factor 177!********************************* 178 179c 180c Calculate H(s) the model exhange hole 181c 182 hnum = ha2*s2 + ha3*s3 + ha4*s4 + ha5*s5 + ha6*s6 + ha7*s7 183 hden = one + hb1*s + hb2*s2 + hb3*s3 + hb4*s4 + hb5*s5 + 184 1 hb6*s6 + hb7*s7 + hb8*s8 + hb9*s9 185 H = hnum/hden 186 187c 188c Calculate helper variables 189c 190 zeta = s2*H 191 eta = A + zeta 192 lambda = D + zeta 193 if (ipol.eq.1) then 194 kf = (three*pi2*rho)**f13 195 else 196 kf = (six*pi2*rho)**f13 197 endif 198 nu = cam_omega/kf 199 chi = nu/dsqrt(lambda+nu**two) 200 lambda2 = (one+chi)*(lambda+nu**two) 201 202c 203c Calculate F(H(s)) for the model exhange hole 204c 205 Fs = one-s2/(r27*C*(one+s02))-zeta/(two*C) 206 207c 208c Calculate EGs 209c 210 EGs = -(two/five)*C*Fs*lambda - (four/r15)*B*lambda*lambda - 211 1 (six/five)*A*lambda*lambda*lambda - 212 2 (four/five)*srpi*lambda**(seven/two) - 213 3 (r12/five)*(lambda**(seven/two))*(dsqrt(zeta)-dsqrt(eta)) 214 215c 216c Calculate the denominators needed 217c 218 219 nu2 = nu*nu 220 Js = (dsqrt(zeta+nu2)+dsqrt(eta+nu2))*(dsqrt(zeta+nu2)+nu) 221 Ks = (dsqrt(zeta+nu2)+dsqrt(eta+nu2))*(dsqrt(eta+nu2)+nu) 222 Ms = (dsqrt(zeta+nu2)+dsqrt(lambda+nu2))*(dsqrt(lambda+nu2)+nu) 223 Ns = (dsqrt(eta+nu2)+dsqrt(lambda+nu2))*(dsqrt(lambda+nu2)+nu) 224 225c 226c The final value for the enhancement factor is 227c 228 tmp1 = one + f12*chi 229 tmp2 = one + (nine/eight)*chi + (three/eight)*chi**two 230 Fxhse1 = A*(zeta/Js + eta/Ks) 231 Fxhse2 = -(four/nine)*B/lambda2 232 Fxhse3 = -(four/nine)*C*Fs*tmp1/lambda2**two 233 Fxhse4 = -(eight/nine)*EGs*tmp2/lambda2**three 234 Fxhse5 = two*zeta*dlog(one -D/Ms) 235 Fxhse6 = -two*eta*dlog(one -(D-A)/Ns) 236 237 Fxhse = Fxhse1+Fxhse2+Fxhse3+Fxhse4+Fxhse5+Fxhse6 238 239!********************************************************* 240! Calculate the First Derivative of the Enhancement Factor 241!********************************************************* 242 243c 244c Calculate the first derivative of H with respect to the 245c reduced density gradient, s. 246c 247 d1hnum = two*ha2*s + three*ha3*s2 + four*ha4*s3 + 248 1 five*ha5*s4 + six*ha6*s5 + seven*ha7*s6 249 250 d1hden = hb1 + two*hb2*s +three*hb3*s2 + four*hb4*s3 + 251 1 five*hb5*s4 + six*hb6*s5 + seven*hb7*s6 + 252 2 eight*hb8*s7 + nine*hb9*s8 253 d1H = (hden*d1hnum -hnum*d1hden)/hden**two 254 255c 256c calculate first derivative of variables needed with respect to s 257c 258 d1zeta = two*s*H + s2*d1H 259 d1eta = d1zeta 260 d1lambda = d1zeta 261 d10chi = -f12*nu*d1zeta/(lambda + nu2)**f32 262 d10lambda2 = d10chi*(lambda + nu**two) + (one+chi)*d1lambda 263 264c 265c calculate the first derivative of Fs with respect to s 266c 267 d1Fs = -two*s/(r27*C*(one+s02)**two) - d1zeta/(two*C) 268 269c 270c Calculate the first derivate of EGs with respect to s 271c 272 d1EGs = -(two/five)*C*(d1Fs*lambda + Fs*d1lambda) - 273 1 (eight/r15)*B*lambda*d1lambda - 274 2 (r18/five)*A*lambda*lambda*d1lambda - 275 3 (r14/five)*srpi*d1lambda*lambda**f52 - 276 4 (r42/five)*(lambda**f52)* 277 5 d1lambda*(dsqrt(zeta)-dsqrt(eta))- 278 6 (six/five)*(lambda**(seven/two))* 279 7 (d1zeta/dsqrt(zeta)-d1eta/dsqrt(eta)) 280 281c 282c Calculate the first derivate of denominators needed with respect 283c to s 284c 285 tmp1 = (dsqrt(zeta+nu2)+nu)/(dsqrt(eta+nu2)) 286 tmp2 = (dsqrt(eta+nu2)+nu)/(dsqrt(zeta+nu2)) 287 288 d10Js = f12*d1zeta*(two+tmp1+tmp2) 289 d10Ks = d10Js 290 291 tmp3 = (dsqrt(zeta+nu2)+nu)/(dsqrt(lambda+nu2)) 292 tmp4 = (dsqrt(lambda+nu2)+nu)/(dsqrt(zeta+nu2)) 293 d10Ms = f12*d1zeta*(two +tmp3+tmp4) 294 295 tmp5 = (dsqrt(lambda+nu2)+nu)/(dsqrt(eta+nu2)) 296 tmp6 = (dsqrt(eta+nu2)+nu)/(dsqrt(lambda+nu2)) 297 d10Ns = f12*d1zeta*(two + tmp5+tmp6) 298c 299c 300c Calculate the derivative of the 08-Fxhse with respect to s 301c 302 L2 = lambda2*lambda2 303 L3 = lambda2*lambda2*lambda2 304 Js2 = Js*Js 305 Ks2 = Ks*Ks 306 Ms2 = Ms*Ms 307 Ns2 = Ns*Ns 308c 309 d1Fxhse1 = A*( (Js*d1zeta - zeta*d10Js)/(Js2) + 310 1 (Ks*d1zeta - eta*d10Ks)/(Ks2) ) 311c 312 d1Fxhse2 = (four/nine)*B*d10lambda2/L2 313c 314 tmp9 = d10lambda2/lambda2 315 tmp7 = d1Fs - two*Fs*tmp9 316 tmp8 = one + f12*chi 317 tmp10 = f12*Fs*d10chi 318 319 d1Fxhse3 = -(four*C/(nine*L2))*(tmp7*tmp8+tmp10) 320c 321 tmp7 = one + (nine/eight)*chi+(three/eight)*chi*chi 322 tmp8 = (nine/eight)*d10chi + (six/eight)*chi*d10chi 323 324 d1Fxhse4 = -(eight/(nine*L3))*((d1EGs-three*EGs*tmp9)*tmp7 325 1 + EGs*tmp8) 326c 327 d1Fxhse5 = two*d1zeta*dlog(one-D/Ms) + 328 1 two*zeta*D*d10Ms/(Ms2*(one-D/Ms)) 329c 330 d1Fxhse6 = -two*d1eta*dlog(one- (D-A)/Ns) - 331 1 two*eta*(D-A)*d10Ns/(Ns2*(one-(D-A)/Ns)) 332c 333 d10Fxhse = d1Fxhse1+d1Fxhse2+d1Fxhse3+d1Fxhse4+d1Fxhse5+d1Fxhse6 334c 335c 336c Calculate the derivative of 08-Fxhse with respect to nu 337c 338 nu3 = nu2*nu 339c 340c Define all derivatives with respect to nu necessary to 341c evaluate enhancement factor derivatives. 342c 343 d01chi = lambda/(nu2 + lambda)**f32 344 d01lambda2 = (lambda*d01chi)/(one - chi)**two 345 d01Js = (eta*(nu + dsqrt(nu2 + zeta)) + 346 1 (nu + dsqrt(nu2 + eta))* 347 2 (zeta + two*nu*(nu + dsqrt(nu2 + zeta))))/ 348 3 (dsqrt(nu2 + eta)*dsqrt(nu2 + zeta)) 349 d01Ks = d01Js 350 d01Ms = (lambda*(nu + dsqrt(nu2 + zeta)) + 351 1 (nu + dsqrt(nu2 + lambda))* 352 2 (zeta + two*nu*(nu + dsqrt(nu2 + zeta))))/ 353 3 (dsqrt(nu2 + lambda)*dsqrt(nu2 + zeta)) 354 d01Ns = (eta*(nu + dsqrt(nu2 + lambda)) + 355 1 (nu + dsqrt(nu2 + eta))* 356 2 (lambda + two*nu*(nu + dsqrt(nu2 + lambda))))/ 357 3 (dsqrt(nu2 + eta)*dsqrt(nu2 + lambda)) 358c 359 d1Fxhse1 = -((A*(nu + dsqrt(eta + nu2))*zeta)/ 360 1 (dsqrt(eta + nu2)*dsqrt(nu2 + zeta)* 361 2 (nu + dsqrt(nu2 + zeta))* 362 3 (dsqrt(eta + nu2) + dsqrt(nu2 + zeta)))) 363c 364 d1Fxhse2 = -((A*eta*(nu/dsqrt(eta + nu2) + nu/ 365 1 dsqrt(nu2 + zeta)))/ 366 2 ((nu + dsqrt(eta + nu2))* 367 3 (dsqrt(eta + nu2) + dsqrt(nu2 + zeta))**two)) - 368 4 (A*eta*(one + nu/dsqrt(eta + nu2)))/ 369 5 ((nu + dsqrt(eta + nu2))**two* 370 6 (dsqrt(eta + nu2) + dsqrt(nu2 + zeta))) 371c 372 d1Fxhse3 = (four*B)/(nine*(lambda + nu2)**(f32)) 373c 374 d1Fxhse4 = (two*C*Fs)/(three*(lambda + nu2)**(f52)) 375c 376 d1Fxhse5 = (five*EGs*(lambda**two + four*nu3* 377 1 (nu + dsqrt(lambda + nu2)) + 378 2 lambda*nu*(five*nu + three*dsqrt(lambda + nu2))))/ 379 3 (three*(lambda + nu2)**four*(nu + dsqrt(lambda + nu2))**three) 380c 381 d1Fxhse6 = (two*D*zeta*(nu + dsqrt(nu2 + zeta)))/ 382 1 (dsqrt(lambda + nu2)*dsqrt(nu2 + zeta)* 383 2 (-D + lambda + (nu + dsqrt(lambda + nu2))* 384 3 (nu + dsqrt(nu2 + zeta)))) 385c 386 d1Fxhse7 = (two*(A - D)*eta*(nu + dsqrt(eta + nu2)))/ 387 1 (dsqrt(eta + nu2)*dsqrt(lambda + nu2)* 388 2 (A - D + lambda + nu2 + nu*dsqrt(eta + nu2) + 389 3 nu*dsqrt(lambda + nu2) + 390 4 dsqrt(eta + nu2)*dsqrt(lambda + nu2))) 391c 392 d01Fxhse = d1Fxhse1+d1Fxhse2+d1Fxhse3+d1Fxhse4+d1Fxhse5+d1Fxhse6+ 393 1 d1Fxhse7 394c 395 396!********************************************************** 397! Calculate the Second Derivative of the Enhancement Factor 398!********************************************************** 399#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 400c 401c Calculate the second derivative of H 402c 403 d2hnum = two*ha2+six*ha3*s+r12*ha4*s2+r20*ha5*s3+ 404 1 r30*ha6*s4 + r42*ha7*s5 405 406 d2hden = two*hb2+six*hb3*s+r12*hb4*s2+r20*hb5*s3+ 407 1 r30*hb6*s4+r42*hb7*s5+r56*hb8*s6+r72*hb9*s7 408 409 d2H = (hden*(hden*d2hnum-hnum*d2hden) - 410 1 two*(hden*d1hnum-hnum*d1hden)*d1hden)/hden**three 411 412c 413c calculate second derivative of variables needed 414c 415 d2zeta = two*H +four*s*d1H + s2*d2H 416 d2eta = d2zeta 417 d2lambda = d2zeta 418 419 d20chi = -f12*(nu/(lambda+nu2)**(f32))* 420 1 (-f32*d1zeta*d1zeta/(lambda+nu2) +d2zeta) 421 422 d20lambda2 =(one-chi)*(d2lambda-d2lambda*chi+lambda*d20chi)+ 423 1 two*d10chi*(d1lambda-d1lambda*chi+lambda*d10chi) 424 d20lambda2 = d20lambda2/(one-chi)**three 425c 426c calculate the second derivative of Fs 427c 428 d2Fs = -(two/(r27*C))*(one-three*s02)/ 429 1 ((one+s02)**three)-d2zeta/(two*C) 430 431c 432c Calculate the second derivative of EGs 433c 434 tmp7 = -(two/five)*C*(d2Fs*lambda+two*d1Fs*d1lambda+Fs*d2lambda) 435 tmp8 = -(eight/r15)*B*(d1lambda*d1lambda+lambda*d2lambda) 436 tmp9 =-(r18/five)*A*lambda*(two*d1lambda*d1lambda+lambda*d2lambda) 437 tmp10 = -(r14/five)*srpi*(f52*(lambda**f32)*d1lambda*d1lambda 438 1 +(lambda**f52)*d2lambda) 439 tmp11 = -(r42/five)*(f52*(lambda**f32)*d1lambda*d1lambda 440 1 +(lambda**f52)*d2lambda)*(dsqrt(zeta)-dsqrt(eta)) 441 tmp12 = -(r42/five)*(lambda**f52)*d1lambda* 442 1 (d1zeta/dsqrt(zeta)-d1eta/dsqrt(eta)) 443 tmp13 = -(six/five)*(lambda**(seven/two))* 444 1 (-d1zeta*d1zeta/(two*zeta**f32)+d2zeta/dsqrt(zeta) 445 2 +d1eta*d1eta/(two*eta**f32)-d2eta/dsqrt(eta)) 446 447 d2EGs = tmp7+tmp8+tmp9+tmp10+tmp11+tmp12+tmp13 448c 449c Calculate the second derivative of denominators needed 450c 451 tmp7 = two/(dsqrt(zeta+nu2)*dsqrt(eta+nu2)) 452 tmp8 = (dsqrt(eta+nu2)+nu)/(zeta+nu2)**f32 453 tmp9 = (dsqrt(zeta+nu2)+nu)/(eta+nu2)**f32 454 455 d20Js = f12*(d2zeta*(two+tmp1+tmp2) + 456 1 f12*d1zeta*d1zeta*(tmp7-tmp8-tmp9)) 457 458 d20Ks = d20Js 459 460 tmp10 = two/(dsqrt(zeta+nu2)*dsqrt(lambda+nu2)) 461 tmp11 = (dsqrt(lambda+nu2)+nu)/(zeta+nu2)**f32 462 tmp12 = (dsqrt(zeta+nu2)+nu)/(lambda+nu2)**f32 463 d20Ms = f12*(d2zeta*(two+tmp3+tmp4) + 464 1 f12*d1zeta*d1zeta*(tmp10-tmp11-tmp12)) 465 466 tmp13 = two/(dsqrt(eta+nu2)*dsqrt(lambda+nu2)) 467 tmp14 = (dsqrt(lambda+nu2)+nu)/(eta+nu2)**f32 468 tmp15 = (dsqrt(eta+nu2)+nu)/(lambda+nu2)**f32 469 d20Ns = f12*(d2zeta*(two+tmp5+tmp6) + 470 1 f12*d1zeta*d1zeta*(tmp13-tmp14-tmp15)) 471 472 473c 474c Calculate the second derivative of the Fxhse with respect to the 475c reduced density gradient, s. 476c 477 Ms4 = Ms2*Ms2 478 Ns4 = Ns2*Ns2 479c 480 tmp1 = A*(Js2*d2zeta -two*Js*d10Js*d1zeta 481 1 + two*d10Js*d10Js*zeta-Js*d20Js*zeta)/ 482 2 (Js2*Js) 483 tmp2 = A*(Ks2*d2zeta -two*Ks*d10Ks*d1zeta 484 1 + two*d10Ks*d10Ks*eta-Ks*d20Ks*eta)/ 485 2 (Ks2*Ks) 486 d2Fxhse1 = tmp1 +tmp2 487c 488 d2Fxhse2 = (four/nine)*B*(L2*d20lambda2- 489 1 two*lambda2*d10lambda2*d10lambda2)/(L2*L2) 490c 491 tmp4 = d10lambda2/lambda2 492 tmp5 = d20lambda2/lambda2 493 d2Fxhse3 = -((four*C)/(nine*L2))*(d2Fs + 494 1 six*Fs*tmp4**two - two*Fs*tmp5 - 495 2 four*d1Fs*tmp4 + 496 3 f12*((d2Fs*chi+two*d1Fs*d10chi+Fs*d20chi)- 497 4 four*tmp4*(d1Fs*chi+Fs*d10chi) + 498 5 six*chi*Fs*tmp4*tmp4)-chi*Fs*tmp5) 499c 500 tmp1 = -EGs*(eight/nine)/L3 501 tmp2 = (-d1EGs + three*EGs*d10lambda2/lambda2)*(eight/nine)/L3 502 tmp3 = (-d2EGs - r12*EGs*d10lambda2**two/L2 + 503 1 three*(two*d1EGs*d10lambda2 + EGs*d20lambda2)/lambda2) 504 2 *(eight/nine)/L3 505 tmp4 = (one + (nine/eight)*chi + (three/eight)*chi**two) 506 tmp5 = (three*(three + two*chi)*d10chi)/eight 507 tmp6 = (three*(two*d10chi**two + (three + two*chi)*d20chi))/eight 508 509 d2Fxhse4 = (tmp1*tmp6+two*tmp2*tmp5+tmp3*tmp4) 510c 511 tmp1 = (one-D/Ms) 512 513 d2Fxhse5 = two*d2zeta*dlog(tmp1)+four*d1zeta*D*d10Ms/(Ms2*tmp1) 514 1 - two*zeta*(D*d10Ms/(Ms2*tmp1))**two + 515 2 two*zeta*D*(Ms2*d20Ms-two*Ms*d10Ms*d10Ms)/ 516 3 (tmp1*Ms**four) 517c 518 tmp1 = (one-(D-A)/Ns) 519 d2Fxhse6 = -two*d2eta*dlog(tmp1)- 520 1 four*d1eta*(D-A)*d10Ns/(Ns2*tmp1)+ 521 2 two*eta*((D-A)*d10Ns/(Ns2*tmp1))**two - 522 3 two*eta*(D-A)*(Ns2*d20Ns-two*Ns*d10Ns*d10Ns)/ 523 4 (tmp1*Ns**four) 524c 525 d20Fxhse = d2Fxhse1+d2Fxhse2+d2Fxhse3+d2Fxhse4+d2Fxhse5+d2Fxhse6 526c 527c 528c Calculate the second derivative of Fxhse with respect to the 529c parameter nu (cam_omega/kf). 530c 531 nu5 = nu3*nu2 532 nu6 = nu5*nu 533c 534c Calculate second derivatives with respect to nu necessary to 535c simplify the appearance of the code. 536c 537 d02chi = (-three*nu*lambda)/(nu2 + lambda)**f52 538 d02lambda2 = -((lambda*(two*d01chi**two - 539 1 (-one + chi)*d02chi))/ 540 2 (-one + chi)**three) 541 d02Js = two*nu*(one/dsqrt(nu2 + eta) + one/dsqrt(nu2 + zeta))* 542 1 (one + nu/dsqrt(nu2 + zeta)) + 543 2 (zeta*(dsqrt(nu2 + eta) + 544 3 dsqrt(nu2 + zeta)))/(nu2 + zeta)**f32 + 545 4 (nu + dsqrt(nu2 + zeta))* 546 5 (one/dsqrt(nu2 + eta) + one/dsqrt(nu2 + zeta) + 547 6 nu2*(-(nu2 + eta)**(-f32) - 548 7 (nu2 + zeta)**(-f32))) 549 d02Ks = d02Js 550 d02Ms = two*nu*(one + nu/dsqrt(nu2 + lambda))* 551 1 (one/dsqrt(nu2 + lambda) + one/dsqrt(nu2 + zeta)) 552 2 + (lambda*(dsqrt(nu2 + lambda) + 553 3 dsqrt(nu2 + zeta)))/(nu2 + lambda)**f32 + 554 4 (nu + dsqrt(nu2 + lambda))* 555 5 (one/dsqrt(nu2 + lambda) + one/dsqrt(nu2 + zeta) + 556 6 nu2*(-(nu2 + lambda)**(-f32) - 557 7 (nu2 + zeta)**(-f32))) 558 d02Ns = two*nu*(one/dsqrt(nu2 + eta) + 559 1 one/dsqrt(nu2 + lambda))* 560 2 (one + nu/dsqrt(nu2 + lambda)) + 561 3 (lambda*(dsqrt(nu2 + eta) + 562 4 dsqrt(nu2 + lambda)))/(nu2 + lambda)**f32 563 5 + (nu + dsqrt(nu2 + lambda))* 564 6 (one/dsqrt(nu2 + eta) + one/dsqrt(nu2 + lambda) + 565 7 nu2*(-(nu2 + eta)**(-f32) - 566 8 (nu2 + lambda)**(-f32))) 567c 568 L5 = L2*L3 569 Js3 = Js2*Js 570 Ks3 = Ks2*Ks 571c 572 d2Fxhse1 = A*((zeta*(two*d01Js**two - Js*d02Js))/Js3 + 573 1 (eta*(two*d01Ks**two - Ks*d02Ks))/Ks3) 574c 575 d2Fxhse2 = (four*B*(-two*d01lambda2**two + lambda2*d02lambda2))/ 576 1 (nine*L3) 577c 578 d2Fxhse3 = -(ten*C*Fs*nu)/(three*(lambda+nu2)**f72) 579c 580 d2Fxhse4 = (EGs*(-four* 581 1 (eight + nine*chi + three*chi**two)*d01lambda2**two - 582 2 lambda2**two*(two*d01chi**two + (three + two*chi)*d02chi) + 583 3 lambda2*(six*(three + two*chi)*d01chi*d01lambda2 + 584 4 (eight + nine*chi + three*chi**two)*d02lambda2)))/ 585 5 (three*L5) 586c 587 d2Fxhse5 = (two*D*zeta*((D - two*Ms)*d01Ms**two + 588 1 Ms*(-D + Ms)*d02Ms))/((D - Ms)**two*Ms**two) 589c 590 d2Fxhse6 = (two*(A - D)*eta*((-A + D - two*Ns)*d01Ns**two + 591 1 Ns*(A - D + Ns)*d02Ns))/ 592 2 (Ns**two*(A - D + Ns)**two) 593c 594 d02Fxhse = d2Fxhse1 + d2Fxhse2 + d2Fxhse3 + d2Fxhse4 + 595 1 d2Fxhse5 + d2Fxhse6 596c 597c 598c Calculate the mixed partial derivative of the enhancement factor 599c with respect to both the reduced density gradient, s, and the 600c parameter nu (cam_omega/kF). Note that the enhancement factor 601c satisfies continuity requirements and therefore it does not matter 602c what order the derivatives are taken in. 603c 604 nu4 = nu3*nu 605c 606c Calculate mixed partial derivatives for the variables depending 607c on both s and nu. 608c 609 d11chi = (three*nu2*d1lambda)/ 610 1 (two*(nu2 + lambda)**f52) - 611 2 d1lambda/(two*(nu2 + lambda)**f32) 612 d11lambda2 = (d1lambda*d01chi)/ 613 1 (one - chi)**two + 614 2 (two*lambda*d01chi* 615 3 d10chi)/(one - chi)**three + 616 4 (lambda*d11chi)/ 617 5 (one - chi)**two 618 d11Js = (-(nu*eta**two*d1zeta) + 619 1 nu*zeta*((-nu2 - zeta)* 620 2 d1eta + 621 3 nu*(nu + dsqrt(nu2 + eta))* 622 4 d1zeta) + 623 5 eta*((nu2 + zeta)* 624 6 (nu + dsqrt(nu2 + zeta))*d1eta 625 7 + (-nu3 + (nu + dsqrt(nu2 + eta))*zeta)* 626 8 d1zeta))/ 627 9 (two*(nu2 + eta)**f32*(nu2 + zeta)**f32) 628 d11Ks = d11Js 629 d11Ms = (-(nu*lambda**two*d1zeta) + 630 1 nu*zeta*((-nu2 - zeta)* 631 2 d1lambda + 632 3 nu*(nu + dsqrt(nu2 + lambda))* 633 4 d1zeta) + 634 5 lambda*((nu2 + zeta)* 635 6 (nu + dsqrt(nu2 + zeta))* 636 7 d1lambda + 637 8 (-nu3 + (nu + dsqrt(nu2 + lambda))*zeta)* 638 9 d1zeta))/ 639 A (two*(nu2 + lambda)**f32*(nu2 + zeta)**f32) 640 d11Ns = (-(nu*eta**two*d1lambda) + 641 1 nu*lambda*((-nu2 - lambda)* 642 2 d1eta + 643 3 nu*(nu + dsqrt(nu2 + eta))* 644 4 d1lambda) + 645 5 eta*((nu2 + lambda)* 646 6 (nu + dsqrt(nu2 + lambda))* 647 7 d1eta + 648 8 (-nu3 + (nu + dsqrt(nu2 + eta))*lambda)* 649 9 d1lambda))/ 650 A (two*(nu2 + eta)**f32*(nu2 + lambda)**f32) 651c 652 d11Fxhse1 = -((A*(-two*Ks3*zeta*d01Js*d10Js + 653 1 Js*Ks3*(d1zeta*d01Js + zeta*d11Js) + 654 2 Js3*(-two*eta*d01Ks*d10Ks + 655 3 Ks*(d1eta*d01Ks + eta*d11Ks))))/(Js3*Ks3)) 656c 657 d11Fxhse2 = (-two*B*(two*nu2*(nu + dsqrt(nu2 + lambda)) + 658 1 lambda*(two*nu + dsqrt(nu2 + lambda)))*d1lambda)/ 659 2 (three*(nu2 + lambda)**three* 660 3 (nu + dsqrt(nu2 + lambda))**two) 661c 662 d11Fxhse3 = (C*(two*(nu2 + lambda)*d1Fs - five*Fs*d1lambda))/ 663 1 (three*(nu2 + lambda)**(f72)) 664c 665 d11Fxhse4 = (five*(eight*nu4*(nu + dsqrt(nu2 + lambda)) + 666 1 lambda**two*(four*nu + dsqrt(nu2 + lambda)) + 667 2 four*nu2*lambda*(three*nu + two*dsqrt(nu2 + lambda)))* 668 3 (two*(nu2 + lambda)*d1EGs - seven*EGs*d1lambda))/ 669 4 (six*(nu2 + lambda)**five*(nu + dsqrt(nu2 + lambda))**four) 670c 671 d11Fxhse5 = (two*D*(D*zeta*d01Ms*d10Ms + 672 1 Ms**two*(d1zeta*d01Ms + zeta*d11Ms) - 673 2 Ms*(D*d1zeta*d01Ms + 674 3 zeta*(two*d01Ms*d10Ms + D*d11Ms))))/((D - Ms)**two*Ms**two) 675c 676 d11Fxhse6 = (two*(A - D)*((-A + D)*eta*d01Ns*d10Ns + 677 1 Ns**two*(d1eta*d01Ns + eta*d11Ns) + 678 2 Ns*((A - D)*d1eta*d01Ns + 679 3 eta*(-two*d01Ns*d10Ns + (A - D)*d11Ns))))/ 680 4 (Ns**two*(A - D + Ns)**two) 681c 682 d11Fxhse = d11Fxhse1 + d11Fxhse2 + d11Fxhse3 + d11Fxhse4 + 683 1 d11Fxhse5 + d11Fxhse6 684c 685#endif 686!********************************************************* 687! Calculate the Third Derivative of the Enhancement Factor 688!********************************************************* 689#ifdef THIRD_DERIV 690c 691c Calculate the third order derivative of H with respect to s. 692c 693 d3hnum = six*ha3+r24*ha4*s+r60*ha5*s2+ 694 1 r120*ha6*s3 + r210*ha7*s4 695 696 d3hden = six*hb3+r24*hb4*s+r60*hb5*s2+ 697 1 r120*hb6*s3+r210*hb7*s4+r336*hb8*s5+r504*hb9*s6 698 699 d3H = (-six*hnum*d1hden**three + 700 1 six*hden*d1hden*(d1hnum*d1hden + hnum*d2hden) + 701 2 hden**three*d3hnum - hden**two*(three*d1hden*d2hnum + 702 3 three*d1hnum*d2hden + hnum*d3hden))/hden**four 703c 704c Calculate the third order derivatives of the "helper variables" 705c 706 d3zeta = six*d1H + six*s*d2H + s2*d3H 707 d3eta = d3zeta 708 d3lambda = d3zeta 709 d30chi = -(nu*(r15*d1lambda*d1lambda*d1lambda 710 1 - r18*(nu2 + lambda)*d1lambda* 711 2 d2lambda + four*(nu2 + lambda)**two*d3lambda))/ 712 3 (eight*(nu2 + lambda)**f72) 713 d30lambda2 = d3lambda/(one - chi) + 714 1 (three*d2lambda*d10chi)/ 715 2 (-one + chi)**two - 716 3 (three*d1lambda* 717 4 (two*d10chi**two - 718 5 (-one + chi)*d20chi))/ 719 6 (-one + chi)**three + 720 7 (lambda*(six*d10chi**three - 721 8 six*(-one + chi)*d10chi* 722 9 d20chi + 723 A (-one + chi)**two*d30chi))/ 724 B (-one + chi)**four 725c 726c Calculate the third order derivative of Fs 727c 728 d3Fs = -(r256*s*(s2 - four) + 729 1 nine*(s2 + four)**four*d3zeta)/ 730 2 (r18*C*(s2 + four)**four) 731 732c 733c Calculate the third order derivative of EGs 734c 735 tmp1 = (-two*C*(three*d1lambda*d2Fs + 736 1 three*d1Fs*d2lambda + lambda*d3Fs + 737 2 Fs*d3lambda))/five 738 tmp2 = (-four*B*(six*d1lambda*d2lambda + 739 1 two*lambda*d3lambda))/r15 740 tmp3 = (-six*A*(six*d1lambda**three + 741 1 r18*lambda*d1lambda*d2lambda + 742 2 three*lambda**two*d3lambda))/five 743 tmp4 = (-four*srpi*((r105*dsqrt(lambda)* 744 1 d1lambda**three)/eight + (r105*lambda**f32*d1lambda* 745 2 d2lambda)/four + (seven*lambda**f52*d3lambda)/two))/five 746 tmp5 = (-r36*(-d1eta/(two*dsqrt(eta)) + 747 1 d1zeta/(two*dsqrt(zeta)))* 748 2 ((r35*lambda**f32*d1lambda**two)/four + 749 3 (seven*lambda**f52*d2lambda)/two))/five 750 4 - (r126*lambda**f52*d1lambda*(d1eta**two/(four*eta**f32) - 751 5 d1zeta**two/(four*zeta**f32) - d2eta/(two*dsqrt(eta)) + 752 6 d2zeta/(two*dsqrt(zeta))))/five - 753 7 (r12*(-dsqrt(eta) + dsqrt(zeta))* 754 8 ((r105*dsqrt(lambda)*d1lambda**three)/ 755 9 eight + (r105*lambda**f32*d1lambda*d2lambda)/four + 756 A (seven*lambda**f52*d3lambda)/two))/five 757 B - (r12*lambda**f72*((-three*d1eta**three)/(eight*eta**f52) + 758 C (three*d1zeta**three)/(eight*zeta**f52) + 759 D (three*d1eta*d2eta)/(four*eta**f32) - 760 E (three*d1zeta*d2zeta)/(four*zeta**f32) - 761 F d3eta/(two*dsqrt(eta)) + d3zeta/(two*dsqrt(zeta))))/five 762 763 d3EGs = tmp1 + tmp2 + tmp3 + tmp4 + tmp5 764 765c 766c Calculate derivatives of denominators needed for third derivatives 767c 768 d30Js = (three*(nu + dsqrt(nu2 + zeta))* 769 1 d1eta**three - 770 2 (three*(nu2 + eta)*d1eta**two* 771 3 d1zeta)/dsqrt(nu2 + zeta) - 772 4 (three*(nu2 + eta)*d1eta* 773 5 ((nu2 + eta)*d1zeta**two + 774 6 two*(nu2 + zeta)* 775 7 ((zeta + nu*(nu + dsqrt(nu2 + zeta)))* 776 8 d2eta - 777 9 (nu2 + eta)*d2zeta)))/ 778 A (nu2 + zeta)**f32 + 779 B ((nu2 + eta)**two* 780 C (three*(eta + nu*(nu + dsqrt(nu2 + eta)))* 781 D d1zeta**three + 782 E six*(nu2 + zeta)*d1zeta* 783 F ((nu2 + zeta)*d2eta - 784 G (eta + nu*(nu + dsqrt(nu2 + eta)))* 785 H d2zeta) + 786 I four*(nu2 + zeta)**two* 787 J ((zeta + nu*(nu + dsqrt(nu2 + zeta)))* 788 K d3eta + 789 L (nu2 + eta + nu*dsqrt(nu2 + eta) + 790 M two*dsqrt(nu2 + eta)* 791 N dsqrt(nu2 + zeta))* 792 O d3zeta)))/ 793 P (nu2 + zeta)**f52)/(eight*(nu2 + eta)**f52) 794 d30Ks = d30Js 795 d30Ms = (three*(nu + dsqrt(nu2 + zeta))* 796 1 d1lambda**three - 797 2 (three*(nu2 + lambda)*d1lambda**two* 798 3 d1zeta)/dsqrt(nu2 + zeta) - 799 4 (three*(nu2 + lambda)*d1lambda* 800 5 ((nu2 + lambda)*d1zeta**two + 801 6 two*(nu2 + zeta)* 802 7 ((zeta + nu*(nu + dsqrt(nu2 + zeta)))* 803 8 d2lambda - 804 9 (nu2 + lambda)*d2zeta)))/ 805 A (nu2 + zeta)**f32 + 806 B ((nu2 + lambda)**two* 807 C (three*(lambda + nu*(nu + dsqrt(nu2 + lambda)))* 808 D d1zeta**three + 809 E six*(nu2 + zeta)*d1zeta* 810 F ((nu2 + zeta)*d2lambda - 811 G (lambda + 812 H nu*(nu + dsqrt(nu2 + lambda)))* 813 I d2zeta) + 814 J four*(nu2 + zeta)**two* 815 K ((nu2 + zeta + nu*dsqrt(nu2 + zeta) + 816 L two*dsqrt(nu2 + lambda)* 817 M dsqrt(nu2 + zeta))* 818 N d3lambda + 819 O (lambda + 820 P nu*(nu + dsqrt(nu2 + lambda)))* 821 Q d3zeta)))/ 822 R (nu2 + zeta)**f52)/(eight*(nu2 + lambda)**f52) 823 d30Ns = (three*(nu + dsqrt(nu2 + lambda))* 824 1 d1eta**three - 825 2 (three*(nu2 + eta)*d1eta**two* 826 3 d1lambda)/dsqrt(nu2 + lambda) 827 4 - (three*(nu2 + eta)*d1eta* 828 5 ((nu2 + eta)*d1lambda**two + 829 6 two*(nu2 + lambda)* 830 7 ((lambda + 831 8 nu*(nu + dsqrt(nu2 + lambda)))* 832 9 d2eta - 833 A (nu2 + eta)*d2lambda)))/ 834 B (nu2 + lambda)**f32 + 835 C ((nu2 + eta)**two* 836 D (three*(eta + nu*(nu + dsqrt(nu2 + eta)))* 837 E d1lambda**three + 838 F six*(nu2 + lambda)*d1lambda* 839 G ((nu2 + lambda)*d2eta - 840 H (eta + nu*(nu + dsqrt(nu2 + eta)))* 841 I d2lambda) + 842 J four*(nu2 + lambda)**two* 843 K ((lambda + 844 L nu*(nu + dsqrt(nu2 + lambda)))* 845 M d3eta + 846 N (nu2 + eta + nu*dsqrt(nu2 + eta) + 847 O two*dsqrt(nu2 + eta)* 848 P dsqrt(nu2 + lambda))* 849 Q d3lambda)))/ 850 R (nu2 + lambda)**f52)/(eight*(nu2 + eta)**f52) 851c 852c Calculate derivatives of the enhancement factor with respect 853c to s. 854c 855 L4 = L2*L2 856 L6 = L4*L2 857 Ms3 = Ms2*Ms 858 Ns3 = Ns2*Ns 859 Js4 = Js2*Js2 860 Ks4 = Ks2*Ks2 861c 862 d3Fxhse1 = (A*(Js3*Ks4*d3zeta - 863 1 six*Ks4*zeta*d10Js**three + 864 2 six*Js*Ks4*d10Js* 865 3 (d1zeta*d10Js + zeta*d20Js) - 866 4 Js**two*Ks4*(three*d2zeta*d10Js + 867 5 three*d1zeta*d20Js + 868 6 zeta*d30Js) + Js4* 869 7 (Ks3*d3eta - six*eta*d10Ks**three + 870 8 six*Ks*d10Ks*(d1eta*d10Ks + eta*d20Ks) - 871 9 Ks**two*(three*d2eta*d10Ks + three*d1eta*d20Ks + 872 A eta*d30Ks))))/(Js4*Ks4) 873c 874 d3Fxhse2 = (four*B*(six*d10lambda2**three - 875 1 six*lambda2*d10lambda2*d20lambda2 + 876 2 lambda2**two*d30lambda2))/(nine*L4) 877c 878 d3Fxhse3 = (two*C*(r24*(two + chi)*Fs*d10lambda2**three - 879 1 r18*lambda2*d10lambda2* 880 2 ((two + chi)*d1Fs*d10lambda2 + 881 3 Fs*(d10chi*d10lambda2 + (two + 882 4 chi)*d20lambda2)) - 883 5 L3*((two + chi)*d3Fs + three*d2Fs*d10chi + 884 6 three*d1Fs*d20chi + Fs*d30chi) + 885 7 two*L2*(three*(two + chi)*d2Fs*d10lambda2 + 886 8 three*d1Fs*(two*d10chi*d10lambda2 + 887 9 (two + chi)*d20lambda2) + 888 A Fs*(three*d10lambda2*d20chi + three*d10chi*d20lambda2 + 889 B (two + chi)*d30lambda2))))/(nine*L5) 890c 891 tmp1 = -EGs*(eight/nine)/L3 892 tmp2 = (-d1EGs + (three*EGs*d10lambda2)/(lambda2))*(eight/nine)/L3 893 tmp3 = (-d2EGs + (six*d1EGs*d10lambda2)/(lambda2) 894 1 - (EGs*((r12*d10lambda2**two)/L2 - 895 2 (three*d20lambda2)/lambda2)))*(eight/nine)/L3 896 tmp4 = (-d3EGs + (nine*d2EGs*d10lambda2)/lambda2 - 897 1 three*(d1EGs*((r12*d10lambda2**two)/L2 - 898 2 (three*d20lambda2)/lambda2)) - (EGs* 899 3 ((-r60*d10lambda2**three)/L3 + 900 4 (r36*d10lambda2*d20lambda2)/L2 901 5 - (three*d30lambda2)/lambda2)))*(eight/nine)/L3 902 tmp5 = (one + (nine/eight)*chi + (three/eight)*chi*chi) 903 tmp6 = (nine*d10chi)/eight + (three*chi*d10chi)/four 904 tmp7 = (nine*d20chi)/eight + (three*(two*d10chi**two + 905 1 two*chi*d20chi))/eight 906 tmp8 = (nine*d30chi)/eight + (three*(six*d10chi* 907 1 d20chi + two*chi*d30chi))/eight 908 909 910 d3Fxhse4 = (tmp1*tmp8+3d0*tmp2*tmp7+3d0*tmp3*tmp6+ 911 1 tmp4*tmp5) 912c 913 d3Fxhse5 = two*dlog(one - D/Ms)*d3zeta - 914 1 (six*D*d2zeta*d10Ms)/(D*Ms - Ms**two) + 915 2 (six*D*d1zeta*((D - two*Ms)*d10Ms**two + Ms*(-D + Ms)*d20Ms))/ 916 3 ((D - Ms)**two*Ms**two) - 917 4 (two*D*zeta* 918 5 (two*(D**two - three*D*Ms + three*Ms**two)*d10Ms**three - 919 6 three*Ms*(D**two - three*D*Ms + two*Ms**two)*d10Ms*d20Ms + 920 7 (D - Ms)**two*Ms**two*d30Ms))/((D - Ms)**three*Ms3) 921c 922 d3Fxhse6 = -two*dlog((A - D + Ns)/Ns)*d3eta + 923 1 (six*(A - D)*d2eta*d10Ns)/(Ns*(A - D + Ns)) + 924 2 (six*(A - D)*d1eta*((-A + D - two*Ns)*d10Ns**two + 925 3 Ns*(A - D + Ns)*d20Ns))/(Ns**two*(A - D + Ns)**two) + 926 4 (two*(A - D)*eta*(two*((A - D)**two + three*(A - D)*Ns 927 5 + three*Ns**two)*d10Ns**three - 928 6 three*Ns*((A - D)**two + three*(A - D)*Ns + 929 7 two*Ns**two)*d10Ns*d20Ns + 930 8 Ns**two*(A - D + Ns)**two*d30Ns))/ 931 9 (Ns3*(A - D + Ns)**three) 932c 933 d30Fxhse = d3Fxhse1 + d3Fxhse2 + d3Fxhse3 + d3Fxhse4 + 934 1 d3Fxhse5 + d3Fxhse6 935c 936c 937c Calculate the third derivative of Fxhse with respect to the 938c parameter nu (cam_omega/kf). 939c 940 nu7 = nu6*nu 941 nu8 = nu7*nu 942c 943c Calculate third derivatives with respect to nu necessary to 944c simplify the appearance of the code. 945c 946 d03chi = (three*(four*nu2 - lambda)*lambda)/(nu2 + lambda)**f72 947 d03lambda2 = (lambda*(six*d01chi**three - 948 1 six*(-one + chi)*d01chi* 949 2 d02chi + 950 3 (-one + chi)**two*d03chi))/ 951 4 (-one + chi)**four 952 d03Js = (three*(-(nu*eta**three*zeta) + 953 1 nu4*(nu + dsqrt(nu2 + eta))*zeta**two - 954 2 nu*eta*zeta* 955 3 (two*nu4 - two*nu*dsqrt(nu2 + eta)*zeta + 956 4 zeta**two) + 957 5 eta**two*(two*nu2*zeta*dsqrt(nu2 + zeta) + 958 6 nu4*(nu + dsqrt(nu2 + zeta)) + 959 7 zeta**two* 960 8 (two*nu + dsqrt(nu2 + eta) + 961 9 dsqrt(nu2 + zeta)))))/ 962 A ((nu2 + eta)**f52*(nu2 + zeta)**f52) 963 d03Ks = d03Js 964 d03Ms = (three*(-(nu*lambda**three*zeta) + 965 1 nu4*(nu + dsqrt(nu2 + lambda))*zeta**two - 966 2 nu*lambda*zeta* 967 3 (two*nu4 - two*nu*dsqrt(nu2 + lambda)*zeta + 968 4 zeta**two) + 969 5 lambda**two*(two*nu2*zeta* 970 6 dsqrt(nu2 + zeta) + 971 7 nu4*(nu + dsqrt(nu2 + zeta)) + 972 8 zeta**two* 973 9 (two*nu + dsqrt(nu2 + lambda) + 974 A dsqrt(nu2 + zeta)))))/ 975 B ((nu2 + lambda)**f52*(nu2 + zeta)**f52) 976 d03Ns = (three*(-(nu*eta**three*lambda) + 977 1 nu4*(nu + dsqrt(nu2 + eta))*lambda**two - 978 2 nu*eta*lambda* 979 3 (two*nu4 - two*nu*dsqrt(nu2 + eta)*lambda + 980 4 lambda**two) + 981 5 eta**two*(two*nu2*lambda* 982 6 dsqrt(nu2 + lambda) + 983 7 nu4*(nu + dsqrt(nu2 + lambda)) + 984 8 lambda**two* 985 9 (two*nu + dsqrt(nu2 + eta) + 986 A dsqrt(nu2 + lambda)))))/ 987 B ((nu2 + eta)**f52*(nu2 + lambda)**f52) 988c 989 d3Fxhse1 = A*(-((zeta*(six*d01Js**three - six*Js*d01Js*d02Js + 990 1 Js**two*d03Js))/Js4) - 991 2 (eta*(six*d01Ks**three - six*Ks*d01Ks*d02Ks + 992 3 Ks**two*d03Ks))/Ks4) 993c 994 d3Fxhse2 = (four*B*(six*d01lambda2**three - 995 1 six*lambda2*d01lambda2*d02lambda2 + 996 2 lambda2**two*d03lambda2))/(nine*L4) 997c 998 d3Fxhse3 = (two*C*Fs*(r24*(two + chi)* 999 1 d01lambda2**three - r18*lambda2*d01lambda2* 1000 2 (d01chi*d01lambda2 + (two + chi)*d02lambda2) 1001 3 - L3*d03chi + two*L2*(three*d01lambda2*d02chi + 1002 4 three*d01chi*d02lambda2 + 1003 5 (two + chi)*d03lambda2)))/(nine*L5) 1004c 1005 d3Fxhse4 = (EGs*(r20*(eight + nine*chi + three*chi**two)* 1006 1 d01lambda2**three - r12*lambda2*d01lambda2* 1007 2 (three*(three + two*chi)*d01chi*d01lambda2 + 1008 3 (eight + nine*chi + three*chi**two)*d02lambda2) - 1009 4 L3*(six*d01chi*d02chi + (three + two*chi)*d03chi) + 1010 5 L2*(r18*d01chi**two*d01lambda2 + nine*(three + two*chi)* 1011 6 d01lambda2*d02chi + nine*(three + two*chi)*d01chi* 1012 7 d02lambda2 + (eight + nine*chi + three*chi**two)* 1013 8 d03lambda2)))/(three*L6) 1014c 1015 d3Fxhse5 = (-two*D*zeta*(two*(D**two - three*D*Ms + 1016 1 three*Ms**two)*d01Ms**three - 1017 2 three*Ms*(D**two - three*D*Ms + two*Ms**two)*d01Ms*d02Ms + 1018 3 (D - Ms)**two*Ms**two*d03Ms))/((D - Ms)**three*Ms3) 1019c 1020 d3Fxhse6 = (two*(A - D)*eta*(two*((A - D)**two + 1021 1 three*(A - D)*Ns + three*Ns**two)*d01Ns**three - 1022 2 three*Ns*((A - D)**two + 1023 3 three*(A - D)*Ns + two*Ns**two)*d01Ns*d02Ns + 1024 4 Ns**two*(A - D + Ns)**two*d03Ns))/ 1025 5 (Ns3*(A - D + Ns)**three) 1026c 1027 d03Fxhse = d3Fxhse1 + d3Fxhse2 + d3Fxhse3 + d3Fxhse4 + 1028 1 d3Fxhse5 + d3Fxhse6 1029c 1030c Calculate the mixed third derivative of Fxhse (dnu ds^2) 1031c 1032c 1033c Calculate mixed partial derivatives for the variables depending 1034c on both s and nu. 1035c 1036 d21chi = (three*(-four*nu2 + lambda)*d1lambda**two + 1037 1 two*(two*nu4 + nu2*lambda - lambda**two)*d2lambda)/ 1038 2 (four*(nu2 + lambda)**f72) 1039 d21lambda2 = ((-one + chi)**two*d2lambda*d01chi + 1040 1 two*(-one + chi)*d1lambda*(-two*d01chi*d10chi + 1041 2 (-one + chi)*d11chi) + 1042 3 lambda*(two*d01chi*(three*d10chi**two - 1043 4 (-one + chi)*d20chi) + 1044 5 (-one + chi)*(-four*d10chi*d11chi + (-one + chi)*d21chi)))/ 1045 6 (-one + chi)**four 1046 d21Js = (-(nu*(nu2 + eta)*(nu2 + zeta)**two*d1eta**two) + 1047 1 three*nu2*(nu2 + zeta)**f52*d1eta**two - 1048 2 (nu2 + eta)*(nu2 + zeta)**f52*d1eta**two + 1049 3 three*nu*(nu2 + zeta)**three*d1eta**two - 1050 4 two*nu*(nu2 + eta)**two*(nu2 + zeta)*d1eta*d1zeta - 1051 5 two*nu*(nu2 + eta)*(nu2 + zeta)**two*d1eta*d1zeta + 1052 6 three*nu2*(nu2 + eta)**f52*d1zeta**two + 1053 7 three*nu*(nu2 + eta)**three*d1zeta**two - 1054 8 nu*(nu2 + eta)**two*(nu2 + zeta)*d1zeta**two - 1055 9 (nu2 + eta)**f52*(nu2 + zeta)*d1zeta**two + 1056 A two*nu*(nu2 + eta)**two*(nu2 + zeta)**two*d2eta - 1057 B two*nu2*(nu2 + eta)*(nu2 + zeta)**f52*d2eta + 1058 C two*(nu2 + eta)**two*(nu2 + zeta)**f52*d2eta - 1059 D two*nu*(nu2 + eta)*(nu2 + zeta)**three*d2eta - 1060 E two*nu2*(nu2 + eta)**f52*(nu2 + zeta)*d2zeta - 1061 F two*nu*(nu2 + eta)**three*(nu2 + zeta)*d2zeta + 1062 G two*nu*(nu2 + eta)**two*(nu2 + zeta)**two*d2zeta + 1063 H two*(nu2 + eta)**f52*(nu2 + zeta)**two*d2zeta)/ 1064 I (four*(nu2 + eta)**f52*(nu2 + zeta)**f52) 1065 d21Ks = d21Js 1066 d21Ms = -(nu*(nu2 + lambda)*(nu2 + zeta)**two*d1lambda**two - 1067 1 three*nu2*(nu2 + zeta)**f52*d1lambda**two + 1068 2 (nu2 + lambda)*(nu2 + zeta)**f52*d1lambda**two - 1069 3 three*nu*(nu2 + zeta)**three*d1lambda**two + 1070 4 two*nu*(nu2 + lambda)**two*(nu2 + zeta)*d1lambda*d1zeta + 1071 5 two*nu*(nu2 + lambda)*(nu2 + zeta)**two*d1lambda*d1zeta - 1072 6 three*nu2*(nu2 + lambda)**f52*d1zeta**two - 1073 7 three*nu*(nu2 + lambda)**three*d1zeta**two + 1074 8 nu*(nu2 + lambda)**two*(nu2 + zeta)*d1zeta**two + 1075 9 (nu2 + lambda)**f52*(nu2 + zeta)*d1zeta**two - 1076 A two*nu*(nu2 + lambda)**two*(nu2 + zeta)**two*d2lambda + 1077 B two*nu2*(nu2 + lambda)*(nu2 + zeta)**f52*d2lambda - 1078 C two*(nu2 + lambda)**two*(nu2 + zeta)**f52*d2lambda + 1079 D two*nu*(nu2 + lambda)*(nu2 + zeta)**three*d2lambda + 1080 E two*nu2*(nu2 + lambda)**f52*(nu2 + zeta)*d2zeta + 1081 F two*nu*(nu2 + lambda)**three*(nu2 + zeta)*d2zeta - 1082 G two*nu*(nu2 + lambda)**two*(nu2 + zeta)**two*d2zeta - 1083 H two*(nu2 + lambda)**f52*(nu2 + zeta)**two*d2zeta)/ 1084 I (four*(nu2 + lambda)**f52*(nu2 + zeta)**f52) 1085 d21Ns = -(nu*(nu2 + eta)*(nu2 + lambda)**two*d1eta**two - 1086 1 three*nu2*(nu2 + lambda)**f52*d1eta**two + 1087 2 (nu2 + eta)*(nu2 + lambda)**f52*d1eta**two - 1088 3 three*nu*(nu2 + lambda)**three*d1eta**two + 1089 4 two*nu*(nu2 + eta)**two*(nu2 + lambda)*d1eta*d1lambda + 1090 5 two*nu*(nu2 + eta)*(nu2 + lambda)**two*d1eta*d1lambda - 1091 6 three*nu2*(nu2 + eta)**f52*d1lambda**two - 1092 7 three*nu*(nu2 + eta)**three*d1lambda**two + 1093 8 nu*(nu2 + eta)**two*(nu2 + lambda)*d1lambda**two + 1094 9 (nu2 + eta)**f52*(nu2 + lambda)*d1lambda**two - 1095 A two*nu*(nu2 + eta)**two*(nu2 + lambda)**two*d2eta + 1096 B two*nu2*(nu2 + eta)*(nu2 + lambda)**f52*d2eta - 1097 C two*(nu2 + eta)**two*(nu2 + lambda)**f52*d2eta + 1098 D two*nu*(nu2 + eta)*(nu2 + lambda)**three*d2eta + 1099 E two*nu2*(nu2 + eta)**f52*(nu2 + lambda)*d2lambda + 1100 F two*nu*(nu2 + eta)**three*(nu2 + lambda)*d2lambda - 1101 G two*nu*(nu2 + eta)**two*(nu2 + lambda)**two*d2lambda - 1102 H two*(nu2 + eta)**f52*(nu2 + lambda)**two*d2lambda)/ 1103 I (four*(nu2 + eta)**f52*(nu2 + lambda)**f52) 1104c 1105 d21Fxhse1 = -((A*(six*Ks4*zeta*d01Js*d10Js**two - 1106 1 two*Js*Ks4*(two*d1zeta*d01Js*d10Js + 1107 2 zeta*(two*d10Js*d11Js + d01Js*d20Js)) + 1108 3 Js**two*Ks4*(d2zeta*d01Js + two*d1zeta*d11Js + 1109 4 zeta*d21Js) + Js4* 1110 5 (six*eta*d01Ks*d10Ks**two - 1111 6 two*Ks*(two*d1eta*d01Ks*d10Ks + 1112 7 two*eta*d10Ks*d11Ks + eta*d01Ks*d20Ks) + 1113 8 Ks**two*(d2eta*d01Ks + two*d1eta*d11Ks + 1114 9 eta*d21Ks))))/(Js4*Ks4)) 1115c 1116 d21Fxhse2 = (four*B*(d01lambda2*(six*d10lambda2**two - 1117 1 two*lambda2*d20lambda2) + 1118 2 lambda2*(-four*d10lambda2*d11lambda2 + lambda2*d21lambda2)))/ 1119 3 (nine*L4) 1120c 1121 d21Fxhse3 = (two*C*(r24*(two + chi)*Fs*d01lambda2* 1122 1 d10lambda2**two - six*lambda2* 1123 2 (two*(two + chi)*d1Fs*d01lambda2*d10lambda2 + 1124 3 Fs*(d10lambda2*(d01chi*d10lambda2 + two*(two + chi)* 1125 4 d11lambda2) + d01lambda2*(two*d10chi*d10lambda2 + 1126 5 (two + chi)*d20lambda2))) - L3*(d2Fs*d01chi + two*d1Fs* 1127 6 d11chi + Fs*d21chi) + two*L2*((two + chi)*d2Fs*d01lambda2 + 1128 7 two*d1Fs*(d01lambda2*d10chi + d01chi*d10lambda2 + 1129 8 (two + chi)*d11lambda2) + Fs*(two*d10lambda2*d11chi + 1130 9 two*d10chi*d11lambda2 + d01lambda2*d20chi + d01chi* 1131 A d20lambda2 + two*d21lambda2 + chi*d21lambda2))))/ 1132 B (nine*L5) 1133c 1134 d21Fxhse4 = (-two*lambda2**two*(lambda2*d1EGs - 1135 1 three*EGs*d10lambda2)* 1136 2 (two*d01chi*d10chi + (three + two*chi)*d11chi) + 1137 3 six*(three + two*chi)*lambda2*d10chi* 1138 4 (-four*EGs*d01lambda2*d10lambda2 + 1139 5 lambda2*(d1EGs*d01lambda2 + EGs*d11lambda2)) + 1140 6 three*EGs*lambda2**two*d01lambda2* 1141 7 (two*d10chi**two + (three + two*chi)*d20chi) - 1142 8 (three + two*chi)*lambda2*d01chi* 1143 9 (lambda2**two*d2EGs + r12*EGs*d10lambda2**two - 1144 A three*lambda2*(two*d1EGs*d10lambda2 + EGs*d20lambda2)) - 1145 B EGs*L3*(four*d10chi*d11chi + 1146 C two*d01chi*d20chi + (three + two*chi)*d21chi) + 1147 D (eight + nine*chi + three*chi**two)* 1148 E (r20*EGs*d01lambda2*d10lambda2**two - 1149 F four*lambda2*(two*d1EGs*d01lambda2*d10lambda2 + 1150 G EGs*(two*d10lambda2*d11lambda2 + 1151 H d01lambda2*d20lambda2)) + 1152 I lambda2**two*(d2EGs*d01lambda2 + two*d1EGs*d11lambda2 + 1153 J EGs*d21lambda2)))/(three*L6) 1154c 1155 d21Fxhse5 = (-two*D*(two*D**two*zeta*d01Ms*d10Ms**two - 1156 1 D*Ms*(two*D*d1zeta*d01Ms*d10Ms + 1157 2 zeta*(two*D*d10Ms*d11Ms + 1158 3 d01Ms*(six*d10Ms**two + D*d20Ms))) + 1159 4 Ms**four*(d2zeta*d01Ms + two*d1zeta*d11Ms + 1160 5 zeta*d21Ms) - two*Ms3* 1161 6 (D*d2zeta*d01Ms + 1162 7 two*d1zeta*(d01Ms*d10Ms + D*d11Ms) + 1163 8 zeta*(two*d10Ms*d11Ms + d01Ms*d20Ms + 1164 9 D*d21Ms)) + Ms**two* 1165 A (D**two*d2zeta*d01Ms + 1166 B two*D*d1zeta*(three*d01Ms*d10Ms + D*d11Ms) + 1167 C zeta*(three*d01Ms*(two*d10Ms**two + D*d20Ms) + 1168 D D*(six*d10Ms*d11Ms + D*d21Ms)))))/((D - Ms)**three*Ms3) 1169c 1170 d21Fxhse6 = (two*(A - D)*(Ns**two*(A - D + Ns)**two*d2eta*d01Ns + 1171 1 two*(A - D)*Ns*(A - D + Ns)*d1eta*d01Ns*d10Ns - 1172 2 four*Ns*(A - D + Ns)**two*d1eta*d01Ns*d10Ns + 1173 3 two*Ns**two*(A - D + Ns)**two*d1eta*d11Ns + 1174 4 eta*(d01Ns*(two*((A - D)**two + three*(A - D)*Ns + 1175 5 three*Ns**two)*d10Ns**two - 1176 6 Ns*((A - D)**two + three*(A - D)*Ns + two*Ns**two)*d20Ns) + 1177 7 Ns*(A - D + Ns)*(-two*(A - D + two*Ns)*d10Ns*d11Ns + 1178 8 Ns*(A - D + Ns)*d21Ns))))/(Ns3*(A - D + Ns)**three) 1179c 1180 d21Fxhse = d21Fxhse1 + d21Fxhse2 + d21Fxhse3 + d21Fxhse4 + 1181 1 d21Fxhse5 + d21Fxhse6 1182c 1183c Calculate the mixed third derivative of Fxhse (dnu^2 ds) 1184c 1185c 1186c Calculate mixed partial derivatives for the variables depending 1187c on both s and nu. 1188c 1189 d12chi = (three*nu*(-two*nu2 + three*lambda)*d1lambda)/ 1190 1 (two*(nu2 + lambda)**f72) 1191 d12lambda2 = ((-one + chi)*d1lambda*(-two*d01chi**two + 1192 1 (-one + chi)*d02chi) + 1193 2 lambda*(six*d01chi**two*d10chi - 1194 3 four*(-one + chi)*d01chi*d11chi + 1195 4 (-one + chi)*(-two*d02chi*d10chi + (-one + chi)*d12chi)))/ 1196 5 (-one + chi)**four 1197 d12Js = (eta**three*(two*nu2 - zeta)*d1zeta + 1198 1 nu2*zeta*((three*nu4 + five*nu2*zeta + two*zeta**two)*d1eta - 1199 2 three*nu3*(nu + dsqrt(nu2 + eta))*d1zeta) - 1200 3 eta*((nu2 + zeta)*(zeta**two + three*nu3*(nu + 1201 4 dsqrt(nu2 + zeta)) + 1202 5 nu*zeta*(two*nu + three*dsqrt(nu2 + zeta)))* 1203 6 d1eta - nu2*(three*nu4 - nu*(five*nu + 1204 7 six*dsqrt(nu2 + eta))*zeta + zeta**two)*d1zeta) + 1205 8 eta**two*(five*nu4*d1zeta + zeta**two*(d1eta + d1zeta) + 1206 9 nu*zeta*(nu*d1eta - three*(nu + dsqrt(nu2 + eta))*d1zeta)))/ 1207 A (two*(nu2 + eta)**f52*(nu2 + zeta)**f52) 1208 d12Ks = d12Js 1209 d12Ms = (lambda**three*(two*nu2 - zeta)*d1zeta + 1210 1 nu2*zeta*((three*nu4 + five*nu2*zeta + 1211 2 two*zeta**two)*d1lambda - three*nu3*(nu + dsqrt(nu2 + lambda))* 1212 3 d1zeta) - lambda*((nu2 + zeta)*(zeta**two + 1213 4 three*nu3*(nu + dsqrt(nu2 + zeta)) + 1214 5 nu*zeta*(two*nu + three*dsqrt(nu2 + zeta)))*d1lambda - 1215 6 nu2*(three*nu4 - nu*(five*nu + six*dsqrt(nu2 + lambda))*zeta + 1216 7 zeta**two)*d1zeta) + lambda**two*(five*nu4*d1zeta + 1217 8 zeta**two*(d1lambda + d1zeta) + 1218 9 nu*zeta*(nu*d1lambda - three*(nu + dsqrt(nu2 + lambda))* 1219 A d1zeta)))/(two*(nu2 + lambda)**f52*(nu2 + zeta)**f52) 1220 d12Ns = (eta**three*(two*nu2 - lambda)*d1lambda + nu2*lambda* 1221 1 ((three*nu4 + five*nu2*lambda + two*lambda**two)*d1eta - 1222 2 three*nu3*(nu + dsqrt(nu2 + eta))*d1lambda) - 1223 3 eta*((nu2 + lambda)*(lambda**two + 1224 4 three*nu3*(nu + dsqrt(nu2 + lambda)) + 1225 5 nu*lambda*(two*nu + three*dsqrt(nu2 + lambda)))*d1eta - 1226 6 nu2*(three*nu4 - nu*(five*nu + six*dsqrt(nu2 + eta))*lambda + 1227 7 lambda**two)*d1lambda) + eta**two*(five*nu4*d1lambda + 1228 8 lambda**two*(d1eta + d1lambda)+ nu*lambda*(nu*d1eta - 1229 9 three*(nu + dsqrt(nu2 + eta))*d1lambda)))/ 1230 A (two*(nu2 + eta)**f52*(nu2 + lambda)**f52) 1231c 1232 d12Fxhse1 = -((A*(six*Ks4*zeta*d01Js**two*d10Js - 1233 1 two*Js*Ks4*(d1zeta*d01Js**two + 1234 2 zeta*(d02Js*d10Js + two*d01Js*d11Js)) + 1235 3 Js**two*Ks4*(d1zeta*d02Js + zeta*d12Js) + 1236 4 Js4*(six*eta*d01Ks**two*d10Ks - 1237 5 two*Ks*(d1eta*d01Ks**two + eta*d02Ks*d10Ks + 1238 6 two*eta*d01Ks*d11Ks) + 1239 7 Ks**two*(d1eta*d02Ks + eta*d12Ks))))/(Js4*Ks4)) 1240c 1241 d12Fxhse2 = (four*B*(six*d01lambda2**two*d10lambda2 - 1242 1 four*lambda2*d01lambda2*d11lambda2 + 1243 2 lambda2*(-two*d02lambda2*d10lambda2 + lambda2*d12lambda2)))/ 1244 3 (nine*L4) 1245c 1246 d12Fxhse3 = (-five*C*nu*(two*(nu2 + lambda)*d1Fs - 1247 1 seven*Fs*d1lambda))/(three*(nu2 + lambda)**f92) 1248c 1249 d12Fxhse4 = (r20*(eight + nine*chi + three*chi**two)*EGs* 1250 1 d01lambda2**two*d10lambda2 - 1251 2 four*lambda2*((eight + nine*chi + three*chi**two)* 1252 3 d1EGs*d01lambda2**two + 1253 4 EGs*(three*(three + two*chi)*d01lambda2**two*d10chi + 1254 5 (eight + nine*chi + three*chi**two)*d02lambda2*d10lambda2 + 1255 6 two*d01lambda2*(three*(three + two*chi)*d01chi*d10lambda2 + 1256 7 (eight + nine*chi + three*chi**two)*d11lambda2))) - 1257 8 L3*(d1EGs*(two*d01chi**two + 1258 9 (three + two*chi)*d02chi) + 1259 A EGs*(two*d02chi*d10chi + four*d01chi*d11chi + 1260 B (three + two*chi)*d12chi)) + 1261 C lambda2**two*(d1EGs*(six*(three + two*chi)*d01chi*d01lambda2 + 1262 D (eight + nine*chi + three*chi**two)*d02lambda2) + 1263 E EGs*(three*(three + two*chi)*d02lambda2*d10chi + 1264 F six*d01chi**two*d10lambda2 + nine*d02chi*d10lambda2 + 1265 G six*chi*d02chi*d10lambda2 + 1266 H r18*d01lambda2*d11chi + 1267 I r12*chi*d01lambda2*d11chi + 1268 J six*d01chi*(two*d01lambda2*d10chi + 1269 K (three + two*chi)*d11lambda2) + eight*d12lambda2 + 1270 L nine*chi*d12lambda2 + three*chi**two*d12lambda2)))/ 1271 M (three*L6) 1272c 1273 d12Fxhse5 = (-two*D*(two*D**two*zeta*d01Ms**two*d10Ms - 1274 1 D*Ms*(D*d1zeta*d01Ms**two + zeta*(six*d01Ms**two* 1275 2 d10Ms + D*d02Ms*d10Ms + two*D*d01Ms*d11Ms)) + 1276 3 Ms4*(d1zeta*d02Ms + zeta*d12Ms) - 1277 4 two*Ms3*(d1zeta*(d01Ms**two + D*d02Ms) + 1278 5 zeta*(d02Ms*d10Ms + two*d01Ms*d11Ms + D*d12Ms)) + 1279 6 Ms2*(D*d1zeta*(three*d01Ms**two + D*d02Ms) + 1280 7 zeta*(six*d01Ms**two*d10Ms + three*D*d02Ms*d10Ms + 1281 8 six*D*d01Ms*d11Ms + D**two*d12Ms))))/ 1282 9 ((D - Ms)**three*Ms3) 1283c 1284 d12Fxhse6 = (two*(A - D)*(Ns*(A - D + Ns)*d1eta* 1285 1 ((-A + D - two*Ns)*d01Ns**two + Ns*(A - D + Ns)*d02Ns) + 1286 2 (A - D)*Ns*eta*(-two*d01Ns**two + (A - D + Ns)*d02Ns)* 1287 3 d10Ns + two*(A - D)*(A - D + Ns)*eta*d01Ns* 1288 4 (-two*d01Ns*d10Ns + Ns*d11Ns) + 1289 5 (A - D + Ns)**two*eta*(six*d01Ns**two*d10Ns - 1290 6 four*Ns*d01Ns*d11Ns + 1291 7 Ns*(-two*d02Ns*d10Ns + Ns*d12Ns))))/ 1292 8 (Ns3*(A - D + Ns)**three) 1293c 1294 d12Fxhse = d12Fxhse1 + d12Fxhse2 + d12Fxhse3 + d12Fxhse4 + 1295 1 d12Fxhse5 + d12Fxhse6 1296#endif 1297 1298 end 1299#ifndef SECOND_DERIV 1300#define SECOND_DERIV 1301c 1302c Compile source again for the 2nd derivative case 1303c 1304#include "hse08fx.F" 1305c 1306#endif 1307#ifndef THIRD_DERIV 1308#define THIRD_DERIV 1309c 1310c Compile source again for the 3rd derivative case 1311c 1312#include "hse08fx.F" 1313#endif 1314