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/ggaexc.f | 76 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 src/labat/ggaexc.f (limited to 'src/labat/ggaexc.f') diff --git a/src/labat/ggaexc.f b/src/labat/ggaexc.f new file mode 100644 index 0000000..135b42e --- /dev/null +++ b/src/labat/ggaexc.f @@ -0,0 +1,76 @@ + subroutine ggaexc(r,db,i,imin,imax,doned,exc) +c ggaexc calculates the exchange-correlation energy +c input r : position coordinate (array) +c input db : spin up and down density. (array) +c input i : counter, i.e. present position = r(i) +c input imin : min value of i +c input imax : max value of i +c in/output doned : true if density is splined +c output exc : exchange-correlation energy + implicit logical (a-z) + logical doned + double precision r,db,gradd,gradup,graddn,fk,fkup,fkdn,sk,g + double precision sup,sdn,t,exc,ec + double precision onethi,twothi,pi,conkf,d,zet,rs,exup,exdn + double precision ex,eclda,ecgga + integer i,imin,imax + dimension r(561),db(2,561) + dimension gradd(561),gradup(561),graddn(561),fk(561),fkup(561) + dimension fkdn(561),sk(561),g(561),sup(561),sdn(561),t(561) + common/gga/gradd,gradup,graddn,fkup,fkdn,sk,g,sup,sdn,t + onethi = 1.d0/3.d0 + twothi = 2.d0/3.d0 + pi = 4.d0*datan(1.d0) + conkf = (3.d0*pi**2)**onethi + d = db(1,i) + db(2,i) + call grad(r,db,i,imin,imax,doned,gradd(i),gradup(i),graddn(i)) + +c calculate the exchange energy, first spin up, then spin down + + if (db(1,i).gt.1.d-100) then + fkup(i) = conkf*(2.d0*db(1,i))**onethi + sup(i) = dabs(gradup(i))/(2.d0*fkup(i)*db(1,i)) + call exchen(2.d0*db(1,i),sup(i),exup) + else + fkup(i) = 0.d0 + sup(i) = 0.d0 + exup = 0.d0 + endif + if (db(2,i).gt.1.d-100) then + fkdn(i) = conkf*(2.d0*db(2,i))**onethi + sdn(i) = dabs(graddn(i))/(2.d0*fkdn(i)*db(2,i)) + call exchen(2.d0*db(2,i),sdn(i),exdn) + else + fkdn(i) = 0.d0 + sdn(i) = 0.d0 + exdn = 0.d0 + endif + if (d.gt.1d-100) then + ex = (exup*db(1,i) + exdn*db(2,i))/d + else + ex = 0.d0 + endif + +c exchange energy done, now calculate the correlation energy + + fk(i) = conkf*d**onethi + sk(i) = dsqrt(4.d0*fk(i)/pi) + if (d.gt.1.d-100) then + zet = (db(1,i) - db(2,i))/d + else + zet = 0.d0 + endif + g(i) = ((1.d0+zet)**twothi+(1.d0-zet)**twothi)/2.d0 + if (d.gt.1.d-18) then + rs = (3.d0/(4.d0*pi*d))**onethi + t(i) = dabs(gradd(i))/(d*2*sk(i)*g(i)) + call ldaec(rs,zet,eclda) + call ggaec(rs,zet,t(i),eclda,ecgga) + ec = eclda + ecgga + else + ec = 0.d0 + endif + exc = ex + ec + + return + end -- cgit v1.2.3