1 files changed, 32 insertions, 0 deletions

diff --git a/src/labat/ldauc.f b/src/labat/ldauc.f
new file mode 100644
index 0000000..4285c9a
--- /dev/null
+++ b/src/labat/ldauc.f
@@ -0,0 +1,32 @@
+      subroutine ldauc(rs,zet,ec,ecrs,eczet,vcup,vcdn)
+c  uniform-gas correlation of perdew and wang 1991
+c  calculates the local correlation potential within the lda approx.
+c  input rs           : seitz radius
+c  input zet          : relative spin polarization 
+c  output ec          : correlation energy
+c  output ecrs        : derivative of ec with respect to rs
+c  output eczet       : derivative of ec with respect to zet
+c  output vcup, vcdn  : up- and down-spin potentials
+      implicit double precision (a-h,o-z)
+      data gam,fzz/0.5198421d0,1.709921d0/
+      data thrd,thrd4/0.333333333333d0,1.333333333333d0/
+      f = ((1.d0+zet)**thrd4+(1.d0-zet)**thrd4-2.d0)/gam
+      call gcor(0.0310907d0,0.21370d0,7.5957d0,3.5876d0,1.6382d0,
+     1     0.49294d0,1.00d0,rs,eu,eurs)
+      call gcor(0.01554535d0,0.20548d0,14.1189d0,6.1977d0,3.3662d0,
+     1     0.62517d0,1.00d0,rs,ep,eprs)
+      call gcor(0.0168869d0,0.11125d0,10.357d0,3.6231d0,0.88026d0,
+     1     0.49671d0,1.00d0,rs,alfm,alfrsm)
+c  alfm is minus the spin stiffness alfc
+      z4 = zet**4
+      ec = eu*(1.d0-f*z4)+ep*f*z4-alfm*f*(1.d0-z4)/fzz
+c  energy done. now the potential:
+      ecrs = eurs*(1.d0-f*z4)+eprs*f*z4-alfrsm*f*(1.d0-z4)/fzz
+      fz = thrd4*((1.d0+zet)**thrd-(1.d0-zet)**thrd)/gam
+      eczet = 4.d0*(zet**3)*f*(ep-eu+alfm/fzz)+fz*(z4*ep-z4*eu
+     1     -(1.d0-z4)*alfm/fzz)
+      comm = ec -rs*ecrs/3.d0-zet*eczet
+      vcup = comm + eczet
+      vcdn = comm - eczet
+      return
+      end