From 5df79c53745fde5d6c3340a2979b1429cd5892c1 Mon Sep 17 00:00:00 2001 From: Henrik Rydberg Date: Sat, 8 Oct 2011 20:30:28 +0200 Subject: Initial import of htcd system 1.0 Signed-off-by: Henrik Rydberg --- src/labat/ggauc.f | 96 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) create mode 100644 src/labat/ggauc.f (limited to 'src/labat/ggauc.f') diff --git a/src/labat/ggauc.f b/src/labat/ggauc.f new file mode 100644 index 0000000..f3adcf5 --- /dev/null +++ b/src/labat/ggauc.f @@ -0,0 +1,96 @@ + subroutine ggauc(rs,zet,t,uu,vv,ww,ec,ecrs,eczet,dvcup,dvcdn) +c gga91 correlation +c input rs : seitz radius +c input zet : relative spin polarization +c input t : abs(grad d)/(d*2.*ks*g) +c input uu : (grad d)*grad(abs(grad d))/(d**2 * (2*ks*g)**3) +c input vv : (laplacian d)/(d * (2*ks*g)**2) +c input ww : (grad d)*(grad zet)/(d * (2*ks*g)**2 +c input ec : correlation energy +c input ecrs : derivative of ec with respect to rs +c input eczet : derivative of ec with respect to zet +c output dvcup,dvcdn : nonlocal parts of correlation potentials + implicit double precision (a-h,o-z) + data xnu,cc0,cx,alf/15.75592d0,0.004235d0,-0.001667212d0,0.09d0/ + data c1,c2,c3,c4/0.002568d0,0.023266d0,7.389d-6,8.723d0/ + data c5,c6,a4/0.472d0,7.389d-2,100.d0/ + data thrdm,thrd2/-0.333333333333d0,0.666666666667d0/ + pi = 4.d0*datan(1.d0) + fk = 1.91915829d0/rs + sk = dsqrt(4.d0*fk/pi) + g = ((1.d0+zet)**thrd2+(1.d0-zet)**thrd2)/2.d0 + bet = xnu*cc0 + delt = 2.d0*alf/bet + g3 = g**3 + g4 = g3*g + pon = -delt*ec/(g3*bet) + b = delt/(dexp(pon)-1.d0) + b2 = b*b + t2 = t*t + t4 = t2*t2 + t6 = t4*t2 + rs2 = rs*rs + rs3 = rs2*rs + q4 = 1.d0+b*t2 + q5 = 1.d0+b*t2+b2*t4 + q6 = c1+c2*rs+c3*rs2 + q7 = 1.d0+c4*rs+c5*rs2+c6*rs3 + cc = -cx + q6/q7 + r0 = (sk/fk)**2 + r1 = a4*r0*g4 + coeff = cc-cc0-3.d0*cx/7.d0 + r2 = xnu*coeff*g3 + r3 = dexp(-r1*t2) + h0 = g3*(bet/delt)*dlog(1.d0+delt*q4*t2/q5) + h1 = r3*r2*t2 + h = h0+h1 +c energy done. now the potential: + ccrs = (c2+2.*c3*rs)/q7 - q6*(c4+2.*c5*rs+3.*c6*rs2)/q7**2 + rsthrd = rs/3.d0 + r4 = rsthrd*ccrs/coeff +c========================================================= +c fix made by Tomas Holmquist + if ((zet.eq.1.d0).or.(zet.eq.-1.d0)) then + gz = 0.d0 + else + gz = ((1.d0+zet)**thrdm - (1.d0-zet)**thrdm)/3.d0 + endif +c========================================================= + fac = delt/b+1.d0 + bg = -3.d0*b2*ec*fac/(bet*g4) + bec = b2*fac/(bet*g3) + q8 = q5*q5+delt*q4*q5*t2 + q9 = 1.d0+2.d0*b*t2 + h0b = -bet*g3*b*t6*(2.d0+b*t2)/q8 + h0rs = -rsthrd*h0b*bec*ecrs + fact0 = 2.d0*delt-6.d0*b + fact1 = q5*q9+q4*q9*q9 + h0bt = 2.d0*bet*g3*t4*((q4*q5*fact0-delt*fact1)/q8)/q8 + h0rst = rsthrd*t2*h0bt*bec*ecrs + h0z = 3.d0*gz*h0/g + h0b*(bg*gz+bec*eczet) + h0t = 2.*bet*g3*q9/q8 + h0zt = 3.d0*gz*h0t/g+h0bt*(bg*gz+bec*eczet) + fact2 = q4*q5+b*t2*(q4*q9+q5) + fact3 = 2.d0*b*q5*q9+delt*fact2 + h0tt = 4.d0*bet*g3*t*(2.d0*b/q8-(q9*fact3/q8)/q8) + h1rs = r3*r2*t2*(-r4+r1*t2/3.d0) + fact4 = 2.d0-r1*t2 + h1rst = r3*r2*t2*(2.d0*r4*(1.d0-r1*t2)-thrd2*r1*t2*fact4) + h1z = gz*r3*r2*t2*(3.d0-4.d0*r1*t2)/g + h1t = 2.d0*r3*r2*(1.d0-r1*t2) + h1zt = 2.d0*gz*r3*r2*(3.d0-11.d0*r1*t2+4.d0*r1*r1*t4)/g + h1tt = 4.d0*r3*r2*r1*t*(-2.d0+r1*t2) + hrs = h0rs+h1rs + hrst = h0rst+h1rst + ht = h0t+h1t + htt = h0tt+h1tt + hz = h0z+h1z + hzt = h0zt+h1zt + comm = h+hrs+hrst+t2*ht/6.d0+7.d0*t2*t*htt/6.d0 + pref = hz-gz*t2*ht/g + fact5 = gz*(2.d0*ht+t*htt)/g + comm = comm-pref*zet-uu*htt-vv*ht-ww*(hzt-fact5) + dvcup = comm + pref + dvcdn = comm - pref + return + end -- cgit v1.2.3