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/ggauxc.f | 82 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 82 insertions(+) create mode 100644 src/labat/ggauxc.f (limited to 'src/labat/ggauxc.f') diff --git a/src/labat/ggauxc.f b/src/labat/ggauxc.f new file mode 100644 index 0000000..d5fcb61 --- /dev/null +++ b/src/labat/ggauxc.f @@ -0,0 +1,82 @@ + subroutine ggauxc(r,db,zet,i,imin,imax,donela,donegr,donez, + j uxcup,uxcdn) +c ggauxc calculates the exchange-correlation potential by taking +c the functional derivative of the energy with respect to the density +c input r : position coordinate (array) +c input db : spin up and down density. (array) +c input zet : relative spin polarization (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 donela : true if r^2*grad(density) is splined +c in/output donegr : true if abs(grad(density)) is splined +c in/output donez : true if zet is splined +c output uxcup : spin-up exchange-correlation potential +c output uxcdn : spin-down exchange-correlation potential + implicit logical (a-z) + double precision r,db,zet,uxcup,uxcdn,gradd,gradup,graddn + double precision fkup,fkdn,sk,g,sup,sdn,t,tempz,z2,gradz + double precision onethi,pi,d,rs,lapld,laplup,lapldn,deriv + double precision gagd,gagup,gagdn,uu,vv,ww,uup,vup,udn,vdn + double precision uclcup,uclcdn,ucgaup,ucgadn,uxup,uxdn + double precision ucup,ucdn,ec,ecrs,eczet + integer i,imin,imax,j + logical donela,donegr,donez + common/gga/gradd,gradup,graddn,fkup,fkdn,sk,g,sup,sdn,t + common/ggaz/tempz,z2 + dimension r(561),db(2,561),zet(561),tempz(2,561),z2(2,561) + dimension gradd(561),gradup(561),graddn(561),fkup(561) + dimension fkdn(561),sk(561),g(561),sup(561),sdn(561),t(561) + onethi = 1.d0/3.d0 + pi = 4.d0*datan(1.d0) + call laplac(r,gradup,graddn,i,imin,imax,donela,lapld,laplup, + j lapldn) + call grabgr(r,gradup,graddn,i,imin,imax,donegr,gagd,gagup,gagdn) + +c calculate the exchange potential + if (db(1,i).gt.1.d-100) then + uup = gradup(i)*gagup/(db(1,i)**2*(2.d0*fkup(i))**3) + vup = laplup/(db(1,i)*(2.d0*fkup(i))**2) + call exchpt(2.d0*db(1,i),sup(i),uup,vup,uxup) + else + uxup = 0.d0 + endif + if (db(2,i).gt.1.d-100) then + udn = graddn(i)*gagdn/(db(2,i)**2*(2.d0*fkdn(i))**3) + vdn = lapldn/(db(2,i)*(2.d0*fkdn(i))**2) + call exchpt(2.d0*db(2,i),sdn(i),udn,vdn,uxdn) + else + uxdn = 0.d0 + endif + +c exchange potential done, now calculate the correlation potential + + if (.not.donez) then + do 10 j = imin,imax + tempz(1,j) = zet(j) + 10 continue + call spline(r,imin+1,imax,tempz,z2) + donez = .true. + endif + gradz = deriv(r,tempz,z2,1,i,imin,imax) + d = db(1,i) + db(2,i) + if (d.gt.1.d-18) then + rs = (3.d0/(4.d0*pi*d))**onethi + uu = gradd(i)*gagd/(d**2*(2.d0*sk(i)*g(i))**3) + vv = lapld/(d*(2.d0*sk(i)*g(i))**2) + ww = gradd(i)*gradz/(d*(2.d0*sk(i)*g(i))**2) + call ldauc(rs,zet(i),ec,ecrs,eczet,uclcup,uclcdn) + call ggauc(rs,zet(i),t(i),uu,vv,ww,ec,ecrs,eczet,ucgaup, + j ucgadn) + ucup = uclcup + ucgaup + ucdn = uclcdn + ucgadn + else + ucdn = 0.d0 + ucup = 0.d0 + endif + uxcup = uxup + ucup + uxcdn = uxdn + ucdn + + return + end + -- cgit v1.2.3