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/ggaec.f | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 src/labat/ggaec.f (limited to 'src/labat/ggaec.f') diff --git a/src/labat/ggaec.f b/src/labat/ggaec.f new file mode 100644 index 0000000..3177859 --- /dev/null +++ b/src/labat/ggaec.f @@ -0,0 +1,43 @@ + subroutine ggaec(rs,zet,t,ec,h) +c called by subroutine ggaexc +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 : correlation energy per electron (ec) +c output h : nonlocal part of correlation energy per electron + 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 thrd2/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 + 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 + return + end -- cgit v1.2.3