1#if defined(FUJITSU_VPP) 2!ocl scalar 3#endif 4#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 5 Subroutine xc_att_xc(rho,ipol,Ex,Amat,Cmat) 6#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 7 Subroutine xc_att_xc_d2(rho,ipol,Ex,Amat,Cmat,Amat2,Cmat2,Cmat3) 8#else 9 Subroutine xc_att_xc_d3(rho,ipol,Ex,Amat,Cmat,Amat2,Cmat2,Cmat3, 10 1 Amat3,Cmat4,Cmat5,Cmat6) 11#endif 12c 13 implicit none 14c 15#include "stdio.fh" 16#include "case.fh" 17c 18 double precision rho, Ex, Amat, Cmat 19 integer ipol 20 21#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 22c Second Derivatives of the Exchange Energy Functional 23 double precision Amat2, Cmat2, Cmat3 24#endif 25#ifdef THIRD_DERIV 26c Third Derivatives of the Exchange Energy Functional 27 double precision Amat3, Cmat4, Cmat5, Cmat6 28#endif 29c 30c 31c References: 32c 33c 34c*************************************************************************** 35c 36 double precision a, b, c, btmp,bfactor 37c 38 double precision a_first,a2_first,btmp_first, btmp1 39c 40 double precision sqrt_pi,t1,t2,t3,t4,t5,t6,t7 41 double precision alpha,beta, DERF 42 double precision f10, f01, b_first 43 double precision a2, a3, a4, a5, a6, a7, a8, a9, a10, a11 44 double precision ta, ta2, ta3, ta4, ta5, ta6, ta7, ta8, ta9, 45 1 ta10 46 double precision f43, f23 47 double precision expf, erff 48 49 Parameter (sqrt_pi = 1.77245385090552d0) 50 Parameter (t7 = 2.666666666666666667d0) 51 Parameter (f43 = 4.0d0/3.0d0) 52 Parameter (f23 = 2.0d0/3.0d0) 53c 54#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 55c Second Derivatives of the Exchange Energy Functional 56 double precision a_second, a2_second, f20 57 double precision b_second, btmp_second, t8 58 double precision a3_second 59 double precision f11, f02 60#endif 61#ifdef THIRD_DERIV 62 double precision a_third, a2_third, a3_third, a4_third 63 double precision f30, f21, f12, f03, f02a 64 double precision b_third, btmp_third 65 double precision t9 66#endif 67 68 69c calculate the a_sigma parameter 70 71c write(luout,*) 'alpha',alpha 72c write(luout,*) 'beta',beta 73c write(luout,*) 'mu',mu 74c 75 if (ipol.eq.1) then 76 Ex = Ex/2d0 77 rho = rho/2d0 78 endif 79 if(ex.le.0) then 80 a = cam_omega*sqrt(-2d0*Ex)/(6d0*sqrt_pi*rho) 81 else 82 write(6,*) ' negative Ex ',Ex 83 call errquit(' xc_att_xc: negative Ex ',0,0) 84 endif 85 alpha = cam_alpha 86 beta = cam_beta 87c 88 f10 = Amat/(2d0*Ex) -1d0/rho 89 a_first = f10*a 90 f01 = Cmat/(2d0*Ex) 91 a2_first = f01*a 92#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 93 f20 = Amat2/(2d0*Ex) - Amat*Amat/(2d0*Ex*Ex) 94 & + 1d0/(rho*rho) 95 f11 = Cmat2 - Amat*Cmat/Ex 96 f11 = f11/(2.0d0*Ex) 97c 98 f02 = Cmat3 - Cmat*Cmat/(2.0d0*Ex) 99 f02 = f02/(2.0d0*Ex) 100c 101 a_second = a*(f10*f10 + f20) 102 103c a2_second = a*(f10*f01 + Cmat2/(2d0*Ex) 104c & - Amat*Cmat/(2d0*Ex*Ex)) 105 a2_second = a*(f10*f01 + f11) 106 107c a3_second = a*(Cmat3/(2d0*Ex) - Cmat*Cmat/(4d0*Ex*Ex)) 108 a3_second = a*f02 109#endif 110#ifdef THIRD_DERIV 111c Amat3 = drdrdr 112c Cmat4 = drdrdg 113c Cmat5 = drdgdg 114c Cmat6 = dgdgdg 115c f02a = Cmat3/(2.0d0*Ex) - Cmat*Cmat/(2.0d0*Ex*Ex) 116 f02a = Cmat3 - Cmat*Cmat/Ex 117 f02a = f02a/(2.0d0*Ex) 118c 119 f30 = Amat3/(2.0d0*Ex) 120 1 - 3.0d0*Amat2*Amat/(2.0d0*Ex*Ex) 121 2 + Amat*Amat*Amat/(Ex*Ex*Ex) 122 3 - 2.0d0/(rho*rho*rho) 123c 124 f21 = Cmat4/(2.0d0*Ex) 125 1 - Cmat2*Amat/(Ex*Ex) 126 2 - Amat2*Cmat/(2.0d0*Ex*Ex) 127 3 + Amat*Amat*Cmat/(Ex*Ex*Ex) 128c 129 f12 = Cmat5/(2.0d0*Ex) 130 1 - Cmat2*Cmat/(Ex*Ex) 131 2 - Amat*Cmat3/(2.0d0*Ex*Ex) 132 3 + Amat*Cmat*Cmat/(Ex*Ex*Ex) 133c 134 f03 = Cmat6/(2.0d0*Ex) 135 1 - Cmat3*Cmat/(Ex*Ex) 136 2 + Cmat*Cmat*Cmat/(2.0d0*Ex*Ex*Ex) 137c 138 a_third = a*( f10*f10*f10 + 3.0d0*f10*f20 + f30 ) 139c 140 a2_third = a*( f10*f10*f01 + f20*f01 + 2.0d0*f10*f11 + f21 ) 141c 142 a3_third = a*( f10*f01*f01 + 2.0d0*f11*f01 + f10*f02a + f12 ) 143c 144 a4_third = a*( f01*f02 + f03 ) 145#endif 146 a2 = a*a 147#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 148 a4 = a2*a2 149 a6 = a4*a2 150 a8 = a6*a2 151 a10 = a8*a2 152#endif 153c 154#ifdef THIRD_DERIV 155 a3 = a2*a 156 a5 = a4*a 157 a7 = a6*a 158 a9 = a8*a 159 a11 = a10*a 160#endif 161 ta = 2d0*a 162 ta2 = ta*ta 163 ta3 = ta2*ta 164 ta4 = ta3*ta 165 ta5 = ta4*ta 166 ta6 = ta5*ta 167 ta7 = ta6*ta 168 ta8 = ta7*ta 169 ta9 = ta8*ta 170 ta10 = ta9*ta 171c 172 if (a .lt. 1.0d-10) then 173 expf = 0.0d0 174 erff = 1.0d0 175 else 176 expf = exp(-1d0/(4d0*a2)) 177 erff = DERF(1d0/(2d0*a)) 178 endif 179c 180 if (a .lt. 0.14d0) then 181c write(luout,*) 'a is small' 182c OLD CODE 183c a = 2d0*a 184c btmp = 1d0-(4d0/3d0)*sqrt_pi*a + 2d0*a*a 185c & - (2d0/3d0)*a*a*a*a 186c btmp = 1d0-btmp 187c 188c btmp_first = (4d0/3d0)*(-sqrt_pi + 3d0*a + 189c & (2d0*exp(-1/(a*a)) - 2d0)*a*a*a) 190c btmp_first = 2d0*btmp_first 191c a = a /2d0 192c OLD CODE 193 btmp = 1.0d0 - f43*sqrt_pi*ta 194 1 + 2.0d0*ta2 - f23*ta4 195 btmp = 1.0d0 - btmp 196 197 btmp_first = f43*( -sqrt_pi + 3.0d0*ta + 198 & (2.0d0*expf - 2.0d0)*ta3 ) 199 btmp_first = 2.0d0*btmp_first 200 else if (a .lt. 4.25d0) then 201c write(luout,*) 'a is medium' 202c stop 203 b = expf - 1d0 204c c = 2d0*a*a*b + 0.5d0 205 c = 2d0*a2*b + 0.5d0 206c btmp = (8d0/3d0)*a*(sqrt_pi*DERF(1/(2d0*a)) + 2d0*a*(b-c)) 207 btmp = (8d0/3d0)*a*(sqrt_pi*erff + 2d0*a*(b-c)) 208c OLD CODE 209c t1 = 1/a 210c t2 = a*a 211c t3 = 1/t2 212cc t4 = exp(-0.25d0*t3) 213c t4 = expf 214c t5 = t4 -1d0 215c t6 = t4 -2d0*t2*t5 - 1.5d0 216c btmp_first = -t7*a * 217c & (2*a*(t4/(2*a**3) - 4d0*a*t5 - t1*t4) + 2d0*t6 -t3*t4) - 218c & t7*(2*a*t6 + sqrt_pi*DERF(0.5d0*t1)) 219c btmp_first = -t7*a * 220c & (2*a*(t4/(2*a**3) - 4d0*a*t5 - t1*t4) + 2d0*t6 -t3*t4) - 221c & t7*(2*a*t6 + sqrt_pi*erff) 222c OLD CODE 223c Daniel (4-12-13) This is a simplified form of what was written above. 224c It is more efficient and likely more stable. 225 btmp_first = -2.0d0*t7*a * 226 & ( -8.0d0*a2*b + expf - 3.0d0 ) - t7*sqrt_pi*erff 227 else 228c write(luout,*) 'a is large' 229c stop 230c OLD CODE 231c a = 2d0*a 232c btmp = 1d0 - 1d0/(9d0*a*a) + 1d0/(60d0*a**4d0) - 233c & 1d0/(420d0*a**6d0) + 1d0/(3240d0*a**8d0) - 234c & 1d0/(27720d0*a**10d0) 235c 236c btmp_first = -1d0/(4.5d0*a**3) + 1d0/(15d0*a**5d0) - 237c & 1d0/(70d0*a**7d0) + 1d0/(405d0*a**9d0) 238c btmp_first = btmp_first*2d0 239c a = a /2d0 240c OLD CODE 241c 242 btmp = 1.0d0 - 1.0d0/(9.0d0*ta2) + 1.0d0/(60.0d0*ta4) 243 1 - 1.0d0/(420.0d0*ta6) + 1.0d0/(3240.0d0*ta8) 244 2 - 1.0d0/(27720.0d0*ta10) 245 246 btmp_first = -1.0d0/(4.5d0*ta3) + 1.0d0/(15.0d0*ta5) 247 1 - 1.0d0/(70.0d0*ta7) + 1.0d0/(405.0d0*ta9) 248 btmp_first = btmp_first*2.0d0 249 end if 250#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 251 if (abs(a) .lt. 1d-40) then 252c btmp_second = 16d0 253 btmp_second = 16.0d0 254c Daniel (4-5-13): Should this be 4.25d0? (This probably doesn't 255c matter since both exponentials and error functions vary slowly 256c as the argument gets large). 257 else if (a .ge. 5d0) then 258c OLD CODE 259c btmp_second = 1d0/(6d0*a**4d0) - 1d0/(48d0*a**6d0) + 260c & 1d0/(640d0*a**8d0) - 1d0/(11520d0*a**10d0) 261c OLD CODE 262c 263 btmp_second = 1.0d0/(6.0d0*a4) - 1.0d0/(48.0d0*a6) 264 1 + 1.0d0/(640.0d0*a8) - 1.0d0/(11520.0d0*a10) 265 266 else 267c OLD CODE 268c t1 = a*a 269c t2 = 1d0/t1 270cc t3 = exp(-0.25d0*t2) 271c t3 = expf 272c t4 = 1d0/(a*a*a) 273c t5 = t3 - 1d0 274c t6 = -t2*t3 275c t8 = -t3/a + 0.5d0*t4*t3 - 4d0*a*t5 276cc 277c btmp_second = -(8d0*a*(2d0*a*(t3/(4*a**6d0) - 278c & 2d0*t3/(a**4d0) +t6 - 4d0*t5) -t3/(2*a**5d0) + 279c & 4d0*t8 + 2d0*t4*t3)/3d0 + 16d0*(2d0*a*t8 + 280c & 2d0*(t3 - 2d0*t1*t5-1.5d0) + t6)/3d0) 281c OLD CODE 282c Daniel (4-12-13): This is a simplified form of what was written above. 283c It is likely more stable and efficient. 284 btmp_second = 16.0d0 - 128.0d0*a2 285 & + (16.0d0 + 128.0d0*a2)*expf 286 end if 287#endif 288#ifdef THIRD_DERIV 289 if (abs(a) .lt. 1.0d-40) then 290 btmp_third = 0.0d0 291 else if (a .ge. 5.0d0) then 292 btmp_third = -2.0d0/(3.0d0*a5) 293 1 + 1.0d0/(8.0d0*a7) 294 2 - 1.0d0/(80.0d0*a9) 295 3 + 1.0d0/(1152.0d0*a11) 296 else 297c OLD CODE 298c t1 = 1d0/(a*a) 299c t2 = expf 300c t3 = t1*t1*t1 301c t4 = t1*t1 302c t5 = t3*a 303c t6 = t2 - 1d0 304c t8 = -t1*t2 - 2d0*t4*t2 + t3*t2/4d0 - 4d0*t6 305c t9 = t4*a 306c 307c btmp_third = -8d0*( a*( 2d0*a*( t2/(8*a**9d0) 308c 1 - 5d0*t2/(2d0*a**7d0) 309c 2 + 15d0*t5*t2/2d0 ) 310c 3 - t2/(4d0*a**8d0) + 6d0*t8 311c 4 - 6d0*t4*t2 + 7d0*t3*t2/2d0 )/3d0 312c 5 + ( 4d0*( -t2/a + t9*t2/2d0 - 4d0*a*t6 ) 313c 6 + 2d0*a*t8 + 2d0*t9*t2 314c 7 - t5*t2/2d0 ) ) 315c OLD CODE 316c 317c Daniel (4-12-13): This is a simplified form of what was written above. 318c It is likely more stable and efficient. 319 btmp_third = 8.0d0*( -32.0d0*a4 320 1 + ( 1.0d0 + 8.0d0*a2 321 2 + 32.0d0*a4 )*expf )/a3 322 endif 323#endif 324 bfactor = 1d0 - alpha - beta*btmp 325 b_first = beta*btmp_first 326#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 327 b_second = beta*btmp_second 328#endif 329#ifdef THIRD_DERIV 330 b_third = beta*btmp_third 331c 332 Amat3 = bfactor*Amat3 333 1 + 3.0d0*Amat2*b_first*a_first 334 2 + 3.0d0*Amat*( b_second*a_first*a_first 335 3 + b_first*a_second ) 336 4 + Ex*( b_third*a_first*a_first*a_first 337 5 + 3.0d0*b_second*a_first*a_second 338 6 + b_first*a_third ) 339c 340 Cmat4 = bfactor*Cmat4 341 1 + 2.0d0*Cmat2*b_first*a_first 342 2 + Amat2*b_first*a2_first 343 3 + 2.0d0*Amat*( b_second*a_first*a2_first 344 4 + b_first*a2_second ) 345 5 + Cmat*( b_second*a_first*a_first 346 6 + b_first*a_second ) 347 7 + Ex*( b_third*a_first*a_first*a2_first 348 8 + b_second*( a2_first*a_second 349 9 + 2.0d0*a_first*a2_second ) 350 A + b_first*a2_third ) 351c 352 Cmat5 = bfactor*Cmat5 353 1 + 2.0d0*Cmat2*b_first*a2_first 354 2 + Amat*( b_second*a2_first*a2_first 355 3 + b_first*a3_second ) 356 4 + Cmat3*b_first*a_first 357 5 + 2.0d0*Cmat*( b_second*a_first*a2_first 358 6 + b_first*a2_second ) 359 7 + Ex*( b_third*a_first*a2_first*a2_first 360 8 + b_second*( a_first*a3_second 361 9 + 2.0d0*a2_first*a2_second ) 362 A + b_first*a3_third ) 363c 364 Cmat6 = bfactor*Cmat6 365 1 + 3.0d0*Cmat3*b_first*a2_first 366 2 + 3.0d0*Cmat*( b_second*a2_first*a2_first 367 3 + b_first*a3_second ) 368 4 + Ex*( b_third*a2_first*a2_first*a2_first 369 5 + 3.0d0*b_second*a2_first*a3_second 370 6 + b_first*a4_third ) 371#endif 372c 373#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 374c b_second = beta*btmp_second 375 Amat2 = bfactor*Amat2 + 2d0*Amat*b_first*a_first 376 & + Ex*b_second*a_first*a_first 377 & + Ex*b_first*a_second 378 379 Cmat2 = bfactor*Cmat2 + Amat*b_first*a2_first 380 & + Cmat*b_first*a_first 381 & + Ex*b_second*a_first*a2_first 382 & + Ex*b_first*a2_second 383 384 Cmat3 = bfactor*Cmat3 + 2d0*Cmat*b_first*a2_first 385 & + Ex*b_second*a2_first*a2_first 386 & + Ex*b_first*a3_second 387#endif 388 Amat = bfactor*Amat + Ex*b_first*a_first 389 Cmat = bfactor*Cmat + Ex*b_first*a2_first 390 Ex = Ex*bfactor 391 392 if (ipol.eq.1) then 393 Ex = 2d0*Ex 394 rho = 2d0*rho 395 endif 396c 397 return 398 end 399#ifndef SECOND_DERIV 400#define SECOND_DERIV 401c 402c Compile source again for the 2nd derivative case 403c 404#include "xc_att_xc.F" 405#endif 406#ifndef THIRD_DERIV 407#define THIRD_DERIV 408c 409c Compile source again for the 3rd derivative case 410c 411#include "xc_att_xc.F" 412#endif 413c $Id$ 414