/* DFT functional repository: xc_blyp fortran77 source
c DFT repository Quantum Chemistry Group
c      subroutine xc_rks_blyp(r,g,f,dfdra,dfdgaa,dfdgab)
c subroutine xc_uks_blyp(ra,rb,gaa,gbb,gab,f,dfdra,dfdrb,
c     +                        dfdgaa,dfdgbb,dfdgab)
c subroutine x_rks_b3(r,g,f,dfdra,dfdgaa)
c subroutine x_uks_b3(ra,rb,gaa,gbb,f,dfdra,dfdrb,dfdgaa,dfdgbb)
c subroutine c_rks_lyp(r,g,f,dfdra,dfdgaa,dfdgab)
c  subroutine c_uks_lyp(ra,rb,gaa,gbb,gab,f,dfdra,dfdrb,
c       &                     dfdgaa,dfdgbb,dfdgab)
c subroutine c_rks_vwnrpa(r,f,dfdra)
c subroutine c_uks_vwnrpa(ra,rb,f,dfdra,dfdrb)
c-------------------------------------------------------------------------------------
c call c_rks_vwnrpa(r,f1,dfdra1)
c call x_rks_b3(r,g,f1,dfdra1,dfdgaa1)
c call c_rks_lyp(r,g,f1,dfdra1,dfdgaa1,dfdgab1)
c call x_uks_b3(ra,rb,gaa,gbb,f,dfdra,dfdrb,dfdgaa,dfdgbb)
c call c_uks_vwnrpa(ra,rb,f1,dfdra1,dfdrb1)
c call c_uks_lyp(ra,rb,gaa,gbb,gab,f1,dfdra1,dfdrb1,
c     +               dfdgaa1,dfdgbb1,dfdgab1)
c-----------------------------------------------------------------------
c----------------------------------------------------------------------- */
#include <math.h>
#include <iostream>
using namespace std;
void x_rks_slater(double r,double &f,double &dfdra);
void xc_rks_blyp(double r,double g,double &f,double &dfdra,double &dfdgaa,double &dfdgab);
void xc_uks_blyp(double ra,double rb,double gaa,double gbb,double gab,double &f,double &dfdra,double &dfdrb,double &dfdgaa,double &dfdgbb,double &dfdgab);
void x_rks_b88(double r,double g,double &f,double &dfdra,double &dfdgaa);
void x_uks_b88(double ra,double rb,double gaa,double gbb,double &f,double &dfdra,double &dfdrb,double &dfdgaa,double &dfdgbb);
void c_rks_lyp(double r,double g,double &f,double &dfdra,double &dfdgaa,double &dfdgab);
void c_uks_lyp(double ra,double rb,double gaa,double gbb,double gab,double &f,double &dfdra,double &dfdrb,double &dfdgaa,double &dfdgbb,double &dfdgab);
void c_rks_vwnrpa(double r,double &f,double &dfdra);
void c_uks_vwnrpa(double ra,double rb,double &f,double &dfdra,double &dfdrb);
//void xc_rks_general(double r,double g,double lyp_wt,double vwnrpa_wt,double &f,double &dfdra,double &dfdgaa,double &dfdgab);
void xc_rks_general(double r,double g,double LDA_wt,double b88_wt,double lyp_wt,double vwnrpa_wt,double pz81_wt,double pbe_wt,double &f,double &dfdra,double &dfdgaa,double &dfdgab);

//void xc_rks_blyp(double r,double g,double &f,double &dfdra,double &dfdgaa,double &dfdgab);
void c_rks_pz81(double r,double &f,double &dfdra);
void c_uks_pz81(double ra,double rb,double &f,double &dfdra,double &dfdrb);
inline double fzeta(double x);
inline double fzzprime(double x);
void x_rks_pbe(double r,double g,double &f,double &dfdra,double &dfdgaa);
void x_uks_pbe(double ra,double rb,double gaa,double gbb,double &f,double &dfdra,double &dfdrb,double &dfdgaa,double &dfdgbb);
//c-----------------------------------------------------------------------

void xc_rks_general(double r,double g,double LDA_wt,double b88_wt,double lyp_wt,double vwnrpa_wt,double pz81_wt,double pbe_wt,double &f,double &dfdra,double &dfdgaa,double &dfdgab){
//      implicit none
/*
c     This subroutine evaluates the blyp exchange-correlation
c     functional [1,2,3] for the closed shell case and its derivative
c     with respect to the alpha density, and the density gradient
c     dot products.
c
c     [1] P.J. Stephens, F.J. Devlin, C.F. Chabalowski, M.J. Frisch,
c         "Ab initio calculation of vibrational absorption and circular
c          dichroism spectra using density functional force fields",
c         J.Phys.Chem. Vol 98 (1994) 11623-11627.
c
c     [2] R.H. Hertwig, W. Koch,
c         "On the parameterization of the local correlation functional:
c          What is Becke-3-LYP?",
c         Chem.Phys.Lett. Vol 268 (1997) 345-351.
c
c     [3] "Gaussian 98 User's Reference
c          Density Functional Methods (DFT) Keywords",
c         M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria,
c         M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery,
c         R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam,
c         A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi,
c         V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli,
c         C. Adamo, S. Clifford, J. Ochterski, G.A. Petersson,
c         P.Y. Ayala, Q. Cui, K. Morokuma, D.K. Malick, A.D. Rabuck,
c         K. Raghavachari, J.B. Foresman, J. Cioslowski, J.V. Ortiz,
c         B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi,
c         R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham,
c         C.Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe,
c         P.M.W. Gill, B.G. Johnson, W. Chen, M.W. Wong, J.L. Andres,
c         M. Head-Gordon, E.S. Replogle, J.A. Pople,
c         Gaussian, Inc., Pittsburgh PA, 1998.
c
c     Parameters:
c
c     r      the total electron density
c     g      the dot product of the total density gradient with itself
c     f      On return the functional value
c     dfdra  On return the derivative of f with respect to the alpha
c            electron density
c     dfdgaa On return the derivative of f with respect to the
c            dot product of the alpha density gradient with itself
c     dfdgab On return the derivative of f with respect to the dot
c            product of the alpha density gradient with the beta density
c            gradient
*/
//      double r, g;
//      double f, dfdra, dfdgaa, dfdgab;
//c
      double f1, dfdra1, dfdgaa1, dfdgab1;
//c
//      double vwnrpa_wt, lyp_wt;
//      const double vwnrpa_wt= 0.0;//default= 0.19;   changed 26 July 2003
//      const double lyp_wt   = 1.0;//default= 0.81;      changed 26 July 2003
//      const double LDA_wt   = 0.0;//0.0;      changed 28 July 2003
//      const double b88_wt   = 1.0;//default= 1.0;      changed 28 July 2003      
//c
      f      = 0.0;
      dfdra  = 0.0;
      dfdgaa = 0.0;
      dfdgab = 0.0;
      if (LDA_wt > 0){
      x_rks_slater(r,f1,dfdra1);
      f      = f      + LDA_wt*f1;
      dfdra  = dfdra  + LDA_wt*dfdra1;
      }
      
      if (vwnrpa_wt > 0 ){
      c_rks_vwnrpa(r,f1,dfdra1);
//      cout<<"slater  r  =   "<< r <<  "   f1  =   "<< f1 << "  dfdra1  =   "<< dfdra1 <<  endl;

      f      = f      + vwnrpa_wt*f1;
      dfdra  = dfdra  + vwnrpa_wt*dfdra1;
      }
//      cout<<"vwnrpa  r  =   "<< r <<  "   f1  =   "<< f1 << "  dfdra1  =   "<< dfdra1 <<  endl;

      if (b88_wt > 0 ){
      x_rks_b88(r,g,f1,dfdra1,dfdgaa1);
      f      = f      + b88_wt*f1;
      dfdra  = dfdra  + b88_wt*dfdra1;
      dfdgaa = dfdgaa + b88_wt*dfdgaa1;
//      cout<<"b88  r  =   "<< r <<  "   g  =   "<< g<<  "   f1  =   "<< f1 << "  dfdra1  =   "<< dfdra1<< "  dfdgaa1  =   "<< dfdgaa1 <<  endl;
      }

      if (lyp_wt > 0 ){
      c_rks_lyp(r,g,f1,dfdra1,dfdgaa1,dfdgab1);
      f      = f      + lyp_wt*f1;
      dfdra  = dfdra  + lyp_wt*dfdra1;
      dfdgaa = dfdgaa + lyp_wt*dfdgaa1;
      dfdgab = dfdgab + lyp_wt*dfdgab1;
//      cout<<"lyp  r  =   "<< r <<  "   g  =   "<< g<<  "   f1  =   "<< f1 << "  dfdra1  =   "<< dfdra1<< "  dfdgaa1  =   "<< dfdgaa1 <<  endl;
      }

      if (pz81_wt > 0 ){      
      c_rks_pz81(r,f1,dfdra1);
      f      = f      + pz81_wt*f1;
      dfdra  = dfdra  + pz81_wt*dfdra1;
      }

      if (pbe_wt > 0 ){
      x_rks_pbe(r,g,f1,dfdra1,dfdgaa1);
      f      = f      + pbe_wt*f1;
      dfdra  = dfdra  + pbe_wt*dfdra1;
      dfdgaa = dfdgaa + pbe_wt*dfdgaa1;
//      cout<<"b88  r  =   "<< r <<  "   g  =   "<< g<<  "   f1  =   "<< f1 << "  dfdra1  =   "<< dfdra1<< "  dfdgaa1  =   "<< dfdgaa1 <<  endl;
      }
      
// Added 11 July 2003 to stop infinities
/* */
if ((r==0.0)){
      f      = 0.0;
      dfdra  = 0.0;
      dfdgaa = 0.0;
      dfdgab = 0.0; }
//      cout<<"Total  f  =   "<< f <<  "   dfdra  =   "<< dfdra << "  dfdgaa  =   "<< dfdgaa<< "  dfdgab  =   "<< dfdgab <<  endl;
//c
}
//c-----------------------------------------------------------------------


void x_rks_slater(double r,double &f,double &dfdra){
//c
//c     A = -3/4*(6/pi)**(1/3)
//      double A
      const double A    = -0.930525736349100013861335177933319;
//c
//      double third, third4;
      const double third  =  0.333333333333333333333333333333333;
      const double third4 =  1.333333333333333333333333333333333;

      double ra, gaa;
      double alpha_rho13, alpha_rho43;
//c
      double  alpha_G;

//c
      ra          = 0.5*r;
      alpha_rho13 = pow(ra,third);
      alpha_rho43 = ra*alpha_rho13;
//c
      if (alpha_rho43 > 0.0){

         alpha_G  = A;

//c        FOURTHIRDS*CD13a*(Ga - Xa*Pa)*fac
// switched to B88 26 June 2003    dfdra = third4*alpha_rho13*(alpha_G - aX*alpha_X*alpha_Gp);
         dfdra = third4*alpha_rho13*alpha_G;
//c

      }else{
         dfdra   = 0.0;

         alpha_G = 0.0;
      }
//c
      f = 2.0*alpha_rho43*alpha_G;
}//      end
//c-----------------------------------------------------------------------

void xc_rks_blyp(double r,double g,double &f,double &dfdra,double &dfdgaa,double &dfdgab){
//      implicit none
/*
c     This subroutine evaluates the blyp exchange-correlation
c     functional [1,2,3] for the closed shell case and its derivative
c     with respect to the alpha density, and the density gradient
c     dot products.
c
c     [1] P.J. Stephens, F.J. Devlin, C.F. Chabalowski, M.J. Frisch,
c         "Ab initio calculation of vibrational absorption and circular
c          dichroism spectra using density functional force fields",
c         J.Phys.Chem. Vol 98 (1994) 11623-11627.
c
c     [2] R.H. Hertwig, W. Koch,
c         "On the parameterization of the local correlation functional:
c          What is Becke-3-LYP?",
c         Chem.Phys.Lett. Vol 268 (1997) 345-351.
c
c     [3] "Gaussian 98 User's Reference
c          Density Functional Methods (DFT) Keywords",
c         M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria,
c         M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery,
c         R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam,
c         A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi,
c         V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli,
c         C. Adamo, S. Clifford, J. Ochterski, G.A. Petersson,
c         P.Y. Ayala, Q. Cui, K. Morokuma, D.K. Malick, A.D. Rabuck,
c         K. Raghavachari, J.B. Foresman, J. Cioslowski, J.V. Ortiz,
c         B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi,
c         R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham,
c         C.Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe,
c         P.M.W. Gill, B.G. Johnson, W. Chen, M.W. Wong, J.L. Andres,
c         M. Head-Gordon, E.S. Replogle, J.A. Pople,
c         Gaussian, Inc., Pittsburgh PA, 1998.
c
c     Parameters:
c
c     r      the total electron density
c     g      the dot product of the total density gradient with itself
c     f      On return the functional value
c     dfdra  On return the derivative of f with respect to the alpha
c            electron density
c     dfdgaa On return the derivative of f with respect to the
c            dot product of the alpha density gradient with itself
c     dfdgab On return the derivative of f with respect to the dot
c            product of the alpha density gradient with the beta density
c            gradient
*/
//      double r, g;
//      double f, dfdra, dfdgaa, dfdgab;
//c
      double f1, dfdra1, dfdgaa1, dfdgab1;
//c
//      double vwnrpa_wt, lyp_wt;
      const double vwnrpa_wt=0.19;
      const double lyp_wt   =0.81;
//c
      f      = 0.0;
      dfdra  = 0.0;
      dfdgaa = 0.0;
      dfdgab = 0.0;
      c_rks_vwnrpa(r,f1,dfdra1);
      f      = f      + vwnrpa_wt*f1;
      dfdra  = dfdra  + vwnrpa_wt*dfdra1;
      x_rks_b88(r,g,f1,dfdra1,dfdgaa1);
      f      = f      + f1;
      dfdra  = dfdra  + dfdra1;
      dfdgaa = dfdgaa + dfdgaa1;
      c_rks_lyp(r,g,f1,dfdra1,dfdgaa1,dfdgab1);
      f      = f      + lyp_wt*f1;
      dfdra  = dfdra  + lyp_wt*dfdra1;
      dfdgaa = dfdgaa + lyp_wt*dfdgaa1;
      dfdgab = dfdgab + lyp_wt*dfdgab1;
// Added 11 July 2003 to stop infinities
/* */
if ((r==0.0)){
      f      = 0.0;
      dfdra  = 0.0;
      dfdgaa = 0.0;
      dfdgab = 0.0; }      
//c
}
//c-----------------------------------------------------------------------
void xc_uks_blyp(double ra,double rb,double gaa,double gbb,double gab,double &f,double &dfdra,double &dfdrb,double &dfdgaa,double &dfdgbb,double &dfdgab){

//      xc_uks_blyp(ra,rb,gaa,gbb,gab,f,dfdra,dfdrb,dfdgaa,dfdgbb,dfdgab)
/*      implicit none
c
c     This subroutine evaluates the B3LYP exchange-correlation
c     functional [1,2,3] for the closed shell case and its derivative
c     with respect to the alpha density, and the density gradient
c     dot products.
c
c     [1] P.J. Stephens, F.J. Devlin, C.F. Chabalowski, M.J. Frisch,
c         "Ab initio calculation of vibrational absorption and circular
c          dichroism spectra using density functional force fields",
c         J.Phys.Chem. Vol 98 (1994) 11623-11627.
c
c     [2] R.H. Hertwig, W. Koch,
c         "On the parameterization of the local correlation functional:
c          What is Becke-3-LYP?",
c         Chem.Phys.Lett. Vol 268 (1997) 345-351.
c
c     [3] "Gaussian 98 User's Reference
c          Density Functional Methods (DFT) Keywords",
c         M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria,
c         M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery,
c         R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam,
c         A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi,
c         V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli,
c         C. Adamo, S. Clifford, J. Ochterski, G.A. Petersson,
c         P.Y. Ayala, Q. Cui, K. Morokuma, D.K. Malick, A.D. Rabuck,
c         K. Raghavachari, J.B. Foresman, J. Cioslowski, J.V. Ortiz,
c         B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi,
c         R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham,
c         C.Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe,
c         P.M.W. Gill, B.G. Johnson, W. Chen, M.W. Wong, J.L. Andres,
c         M. Head-Gordon, E.S. Replogle, J.A. Pople,
c         Gaussian Inc., Pittsburgh PA, 1998.
c
c     Parameters:
c
c     ra     the alpha-electron density
c     rb     the beta-electron density
c     gaa    the dot product of the alpha density gradient with itself
c     gbb    the dot product of the beta density gradient with itself
c     gab    the dot product of the alpha density gradient with
c            the beta density
c     f      On return the functional value
c     dfdra  On return the derivative of f with respect to ra
c     dfdrb  On return the derivative of f with respect to rb
c     dfdgaa On return the derivative of f with respect to gaa
c     dfdgbb On return the derivative of f with respect to gbb
c     dfdgab On return the derivative of f with respect to gab
*/
//      double ra, rb, gaa, gbb, gab;
//      double f, dfdra, dfdrb, dfdgaa, dfdgbb, dfdgab;

      double f1, dfdra1, dfdrb1, dfdgaa1, dfdgbb1, dfdgab1;

//      double vwnrpa_wt, lyp_wt;
      const double vwnrpa_wt=0.19;
      const double lyp_wt   =0.81;

      dfdgab = 0.0;
      x_uks_b88(ra,rb,gaa,gbb,f,dfdra,dfdrb,dfdgaa,dfdgbb);
      c_uks_vwnrpa(ra,rb,f1,dfdra1,dfdrb1);
      f      = f      + vwnrpa_wt*f1;
      dfdra  = dfdra  + vwnrpa_wt*dfdra1;
      dfdrb  = dfdrb  + vwnrpa_wt*dfdrb1;
      c_uks_lyp(ra,rb,gaa,gbb,gab,f1,dfdra1,dfdrb1,dfdgaa1,dfdgbb1,dfdgab1);
      f      = f      + lyp_wt*f1;
      dfdra  = dfdra  + lyp_wt*dfdra1;
      dfdrb  = dfdrb  + lyp_wt*dfdrb1;
      dfdgaa = dfdgaa + lyp_wt*dfdgaa1;
      dfdgab = dfdgab + lyp_wt*dfdgab1;
      dfdgbb = dfdgbb + lyp_wt*dfdgbb1;

if ((ra<=1.0e-142) && (gaa<=1.0e-142)){
      f      = 0.0;
      dfdra  = 0.0;
      dfdrb  = 0.0;
      dfdgaa = 0.0;
      dfdgbb = 0.0;
      dfdgab = 0.0; }

      
}
//c-----------------------------------------------------------------------
//c-----------------------------------------------------------------------
      void x_rks_b88(double r,double g,double &f,double &dfdra,double &dfdgaa){
//      subroutine x_rks_b3(r,g,f,dfdra,dfdgaa)
/*      implicit none
c
c     This subroutine evaluates the exchange part of the B3xxx
c     functionals excluding the exact exchange term for closed shell
c     cases. In the context of Becke's paper [1] this means that this
c     subroutine evaluates the part
c
c         (1-a0) E_X^LSDA + aX dE_X^B88
c
c     of equation (2). To complete the exchange part one has to add
c
c         a0 E_X^exact
c
c     In [1] a0 and aX were given as
c
c         a0 = 0.20
c         aX = 0.72
c
c     In addition to the functional the ingredients of the corresponding
c     potential are calculated as well.
c
c
c     The original code was provided by Dr. Phillip Young.
c
c     [1] A.D. Becke,
c         "Density-functional thermochemistry. III. The role of
c          exact exchange",
c         J.Chem.Phys. Vol. 98 (1993) 5648-5652.
c
c     Parameters:
c
c     r      the total electron density
c     g      the dot product of the total density gradient with itself
c     f      On return the functional value
c     dfdra  On return the derivative of f with respect to the alpha-
c            electron density
c     dfdgaa On return the derivative of f with respect to the
c            dot product of the alpha-density gradient with itself
*/
//      double a0, aX;
// switched to B88 26 June 2003      const double a0 = 0.20;
// switched to B88 26 June 2003      const double aX = 0.72;
//c
//c     A = -3/4*(6/pi)**(1/3)
//      double A
      const double A    = -0.930525736349100013861335177933319;
//c
//      double third, third4;
      const double third  =  0.333333333333333333333333333333333;
      const double third4 =  1.333333333333333333333333333333333;
/*
c...  BeckeP : The parameter beta in Becke's paper
c...  BeckeP6: 6 Times the parameter beta in Becke's paper
*/
//      double BeckeP, BeckeP6;
      const double BeckeP  = 0.0042;
      const double BeckeP6 = 0.0252;
//c
//      double r, g;
//      double f, dfdra, dfdgaa;
      double ra, gaa;
      double alpha_rho13, alpha_rho43;
//c
      double alpha_X, alpha_X2, alpha_Xt, alpha_Xhi, alpha_de, alpha_G;
      double alpha_Gp, gamma_a;
//c
      gaa         = 0.25*g;
      ra          = 0.5*r;
      alpha_rho13 = pow(ra,third);
      alpha_rho43 = ra*alpha_rho13;
//c
      if (alpha_rho43 > 0.0){
         gamma_a     = sqrt(gaa);
//c...     Xa = the reduced density gradient
         alpha_X     = gamma_a/alpha_rho43;
         alpha_X2    = alpha_X * alpha_X;
//c ...    Sa = sqrt(Xa*Xa + ONE)
         alpha_Xt    = sqrt(1.0 +alpha_X2);
//c ...    Ha = log(Xa + Sa)
         alpha_Xhi   = log(alpha_X+alpha_Xt);
//         cout << " ra = " <<  ra << " gaa = " <<  gaa << " alpha_X = " <<  alpha_X << " Arcsinh(x) = " <<  alpha_Xhi <<  endl;
         alpha_Xt    = alpha_X/alpha_Xt;
//c
//c ...    Ta = 1/(ONE + SIX*(BETA*Xa)*Ha)
         alpha_de = 1.0/(1.0 +BeckeP6*alpha_X*alpha_Xhi);
//c ...    Ga = P2LDA + Ca
// switched to B88 26 June 2003   alpha_G  = (1.0-a0)*A-aX*(BeckeP*alpha_X2)*alpha_de;
         alpha_G  = A-(BeckeP*alpha_X2)*alpha_de;         
//c
//c ...    Pa = ( (SIX*(BETA*Xa)**2)*(Xa/Sa - Ha) - TWO*(BETA*Xa) )*Ta**2
//c
         alpha_Gp =(6.0*pow((BeckeP*alpha_X),2)*(alpha_Xt-alpha_Xhi)-
                (2.0*BeckeP*alpha_X) ) * (alpha_de*alpha_de);
//c        FOURTHIRDS*CD13a*(Ga - Xa*Pa)*fac
// switched to B88 26 June 2003    dfdra = third4*alpha_rho13*(alpha_G - aX*alpha_X*alpha_Gp);
         dfdra = third4*alpha_rho13*(alpha_G - alpha_X*alpha_Gp);
//c
         if(gamma_a > 0.0)
            dfdgaa = (0.5/gamma_a)*alpha_Gp; // dfdgaa = (0.5/gamma_a)*aX*alpha_Gp;
         else
            dfdgaa = 0.0;
         //endif
      }else{
         dfdra   = 0.0;
         dfdgaa  = 0.0;
         alpha_G = 0.0;
      }
//c
      f = 2.0*alpha_rho43*alpha_G;

// Added 11 July 2003 to stop infinities
if ((r<=1.0e-142) && (g<=1.0e-142)){
      f      = 0.0;
      dfdra  = 0.0;
      dfdgaa = 0.0;
//      dfdgab = 0.0;
      }      
      
}//      end
//c-----------------------------------------------------------------------
      void x_uks_b88(double ra,double rb,double gaa,double gbb,double &f,double &dfdra,double &dfdrb,double &dfdgaa,double &dfdgbb) {
/*      implicit none
c
c     This subroutine evaluates the exchange part of the B3xxx
c     functionals excluding the exact exchange term. In the context of
c     Becke's paper [1] this means that this subroutine evaluates the
c     part
c
c         (1-a0) E_X^LSDA + aX dE_X^B88
c
c     of equation (2). To complete the exchange part one has to add
c
c         a0 E_X^exact
c
c     In [1] a0 and aX were given as
c
c         a0 = 0.20
c         aX = 0.72
c
c     In addition to the functional the ingredients of the corresponding
c     potential are calculated as well.
c
c
c     The original code was provided by Dr. Phillip Young.
c
c     [1] A.D. Becke,
c         "Density-functional thermochemistry. III. The role of
c          exact exchange",
c         J.Chem.Phys. Vol. 98 (1993) 5648-5652.
c
c     Parameters:
c
c     ra     the alpha-electron density
c     rb     the beta-electron density
c     gaa    the dot product of the alpha density gradient with itself
c     gbb    the dot product of the beta density gradient with itself
c     f      On return the functional value
c     dfdra  On return the derivative of f with respect to ra
c     dfdrb  On return the derivative of f with respect to rb
c     dfdgaa On return the derivative of f with respect to gaa
c     dfdgbb On return the derivative of f with respect to gbb
*/
//      double a0, aX
      const double a0 = 0.20;
      const double aX = 0.72;
//c
//c     A = -3/4*(6/pi)**(1/3)
//      double A
      const double A    = -0.930525736349100013861335177933319;
//c
//      double third, third4
      const double third  =  0.333333333333333333333333333333333;
      const double third4 =  1.333333333333333333333333333333333;
/*c
c...  BeckeP : The parameter beta in Becke's paper
c...  BeckeP6: 6 Times the parameter beta in Becke's paper
c*/
//      double BeckeP, BeckeP6
      const double BeckeP  = 0.0042;
      const double BeckeP6 = 0.0252;
//c
//      double ra, rb, gaa, gbb;
//      double f, dfdra, dfdrb, dfdgaa, dfdgbb;
      double alpha_rho13, alpha_rho43, beta_rho13, beta_rho43;
//c
      double alpha_X, alpha_X2, alpha_Xt, alpha_Xhi, alpha_de, alpha_G;
      double alpha_Gp, gamma_a;
      double beta_X, beta_X2, beta_Xt, beta_Xhi, beta_de, beta_G;
      double beta_Gp, gamma_b;
//c
      alpha_rho13 = pow(ra,third);
      beta_rho13  = pow(rb,third);
      alpha_rho43 = ra*alpha_rho13;
      beta_rho43  = rb*beta_rho13;
//c
      if (alpha_rho43>0.0){
         gamma_a     = sqrt(gaa);
//c...     Xa = the reduced density gradient
         alpha_X     = gamma_a/alpha_rho43;
         alpha_X2    = alpha_X * alpha_X;
//c ...    Sa = sqrt(Xa*Xa + ONE)
         alpha_Xt    = sqrt(1.0+alpha_X2);
//c ...    Ha = log(Xa + Sa)
         alpha_Xhi   = log(alpha_X+alpha_Xt);
         alpha_Xt    = alpha_X/alpha_Xt;
//c
//c ...    Ta = 1/(ONE + SIX*(BETA*Xa)*Ha)
         alpha_de = 1.0/(1.0+BeckeP6*alpha_X*alpha_Xhi);
//c ...    Ga = P2LDA + Ca
//// 27 June 2003 b3 to B88         alpha_G  = (1.0-a0)*A-aX*(BeckeP*alpha_X2)*alpha_de ;
alpha_G  = A-(BeckeP*alpha_X2)*alpha_de;
//c
//c ...    Pa = ( (SIX*(BETA*Xa)**2)*(Xa/Sa - Ha) - TWO*(BETA*Xa) )*Ta**2
//c
         alpha_Gp =(6.0*pow((BeckeP*alpha_X),2)*(alpha_Xt-alpha_Xhi)-
                  (2.0*BeckeP*alpha_X) ) * (alpha_de*alpha_de);
//c        FOURTHIRDS*CD13a*(Ga - Xa*Pa)*fac
// 27 June 2003       dfdra = third4*alpha_rho13*(alpha_G - aX*alpha_X*alpha_Gp);
         dfdra = third4*alpha_rho13*(alpha_G - alpha_X*alpha_Gp);
//c
         if(gamma_a > 0.0) 
            dfdgaa = (0.5/gamma_a)*alpha_Gp ;  // 27 June 2003 dfdgaa = (0.5/gamma_a)*aX*alpha_Gp ;
         else
            dfdgaa = 0.0 ;
         //endif
      }else{
         dfdra   = 0.0;
         dfdgaa  = 0.0 ;
         alpha_G = 0.0;
      }//endif
//c
      if (beta_rho43 > 0.0){
         gamma_b     = sqrt(gbb);
         beta_X      = gamma_b/beta_rho43;
         beta_X2     = beta_X * beta_X ;
         beta_Xt     = sqrt(beta_X*beta_X+1.0);
         beta_Xhi    = log(beta_X+beta_Xt);
         beta_Xt     = beta_X/beta_Xt ;
/*
C        energy
c
c        Tb = 1/(ONE + SIX*(BETA*Xb)*Hb)
*/
         beta_de  = 1.0/(1.0+BeckeP6*beta_X*beta_Xhi);
//c        Gb = P2LDA + Cb
//         beta_G   = (1.0-a0)*A-aX*(BeckeP*beta_X2)*beta_de ;
        beta_G   = A-(BeckeP*beta_X2)*beta_de ;
/*
C        potential and gradient potential
c
c        Pb = ( (SIX*(BETA*Xb)**2)*(Xb/Sb - Hb) - TWO*(BETA*Xb) )*Tb**2
*/
         beta_Gp  =(6.0*pow((BeckeP*beta_X),2)*(beta_Xt -beta_Xhi )-
                 (2.0*BeckeP*beta_X ) ) * (beta_de *beta_de );
//         dfdrb = third4*beta_rho13 *(beta_G   - aX*beta_X *beta_Gp );
dfdrb = third4*beta_rho13 *(beta_G   - beta_X *beta_Gp );
//c
//c        gaa = Pa/GNa*fac
         if(gamma_b > 0.0) 
            dfdgbb = (0.5/gamma_b)*beta_Gp;  //  dfdgbb = (0.5/gamma_b)*aX*beta_Gp;
         else
            dfdgbb = 0.0 ;
        // endif
      }else {
         dfdrb  = 0.0 ;
         dfdgbb = 0.0;
         beta_G = 0.0;
     }// endif
//c
      f = alpha_rho43*alpha_G+beta_rho43*beta_G;
      

      
   }//  end
//c-----------------------------------------------------------------------
//c-----------------------------------------------------------------------
      void c_rks_lyp(double r,double g,double &f,double &dfdra,double &dfdgaa,double &dfdgab){
/*      implicit none
c
c     This subroutine evaluates LYP correlation functional [1] and
c     the ingredients of the corresponding potential for the closed
c     shell case. The implementation is based on the expression by
c     Miehlich et al. [2]. This expression is simpler than the original
c     one in the sense that the second derivatives of the density have
c     been removed through partial integration.
c
c     The original code was provided by Dr. Phillip Young.
c
c     [1] C. Lee, W. Yang, and R.G. Parr,
c         "Development of the Colle-Salvetti correlation-energy
c          formula into a functional of the electron density"
c         Phys.Rev. Vol. B37 (1988) 785-789.
c
c     [2] B. Miehlich, A. Savin, H. Stoll, and H. Preuss,
c         "Results obtained with the correlation energy density
c          functionals of Becke and Lee, Yang and Parr"
c         Chem.Phys.Lett. Vol. 157 (1989) 200-206.
c
c     Parameters:
c
c     r      the total electron density
c     g      the dot product of the total density gradient with itself
c     f      On return the functional value
c     dfdra  On return the derivative of f with respect to the alpha
c            electron density
c     dfdgaa On return the derivative of f with respect to the
c            dot product of the alpha density gradient with itself
c     dfdgab On return the derivative of f with respect to the dot
c            product of the alpha density gradient with the beta density
c            gradient
*/
//      double r, g;
//      double f, dfdra, dfdgaa, dfdgab;
//c
      double rhoa13,rho13,rho43;
      double rhoa83;
      double r_a,rab,raa,rab9,ra,gaa,gab ;
      double faa,fab,dgaa,dgab;
      double d2agaa,d2bgaa;
      double d2agab;
      double dfaaa,dfaab;
      double dfaba;
      double del,delp,ome,ro,xtc,xtd,rx,lt1,lt2a,lt2 ;
      double dlt1a,dlt2a;
      double zeta, cd;
//c
//      double t1,ta,tb,tc,td,t2,t3,t43,t83,t19,n13
      const double n13 = 0.3333333333333333;
      const double t1  = 0.3646239897876487e02; //  ! p1
      const double ta  = 0.04918;               // ! a
      const double tb  = 0.13200;              //  ! b
      const double tc  = 0.25330;              //  ! c
      const double td  = 0.34900;              //  ! d
      const double t2  = ta*tb/td ;            //    ! p2
      const double t3  = 4.0*ta/td;            //  ! p3
      const double t43 = 4.0*n13;
      const double t83 = 8.0*n13;
      const double t19 = 1.0/9.0;

      gaa    = 0.25*g;
      gab    = gaa;   
      ra     = 0.5*r;
      cd     = r ;
      rho13  = pow(cd,n13);
      rhoa13 = pow(ra,n13);

      rho43  = cd*rho13;
      rhoa83 = rhoa13*rhoa13*ra*ra;

      r_a  = 0.5;

      raa = 0.25;
      rab = 0.25;

      zeta = 0.0 ;
//c
      xtc = tc/rho13 ;
      xtd = td/rho13;
//c
      rx = xtd/(1. + xtd);
//c
      del = xtc + rx;
      ome = t2*rx*exp(-xtc)/rho43 ;
//c
      rab9 = t19*rab;
//c
      faa = rab9*(1. - 3.*del - (del - 11.)*r_a);
      fab = rab9*(47. - 7.*del);
//c
      dgaa = -ome*(faa - raa);
      dgab = -ome*(fab - t43);
//c
      lt1  = -t3*rx*rab*rho43;
      lt2a = -(t1*ome*rab)*rhoa83;
      lt2  = 2.0*lt2a;
//c
      f = lt1 + lt2 + 2.0*dgaa*gaa + dgab*gab;
//c
      delp = n13*(rx*rx - del)/cd;
      ro   = n13*(del - 5.)/cd;
//c
      dfaaa = -rab9*( (3. + r_a)*delp + (del - 11.)*r_a/cd );
      dfaab = -rab9*( (3. + r_a)*delp - (del - 11.)*r_a/cd );
      dfaba = - 7.*rab9*delp;
//c
      dlt1a = (lt1/cd)*(n13*rx + 1.0);
      dlt2a = ro*lt2a + (ro + (t83/ra))*lt2a;
//c
      d2agaa = ro*dgaa - ome * (dfaaa + raa*(2./cd));
      d2bgaa = ro*dgaa - ome * (dfaab - rab*(2./cd)) ;
      d2agab = ro*dgab - ome* dfaba ;
//c
      dfdra = dlt1a + dlt2a +
             d2agaa*gaa + d2bgaa*gaa + d2agab*gab;
//c Changed 11 July 2003
      dfdgaa = dgaa;
//      if (r==0.0) dfdgaa = 0.0;
      dfdgab = dgab;
//c
// Added 11 July 2003 to stop infinities
if ((r<=1.0e-142) && (g<=1.0e-142)){
      f      = 0.0;
      dfdra  = 0.0;
      dfdgaa = 0.0;
      dfdgab = 0.0; }

      }//end
//c-----------------------------------------------------------------------
      void c_uks_lyp(double ra,double rb,double gaa,double gbb,double gab,double &f,double &dfdra,double &dfdrb,double &dfdgaa,double &dfdgbb,double &dfdgab){
/*      implicit none
c
c     This subroutine evaluates LYP correlation functional [1] and
c     the ingredient of the corresponding potential. The implementation
c     is based on the expression by Miehlich et al. [2]. This expression
c     is simpler than the original one in the sense that the second
c     derivatives of the density have been removed through partial
c     integration.
c
c     The original code was provided by Dr. Phillip Young.
c
c     [1] C. Lee, W. Yang, and R.G. Parr,
c         "Development of the Colle-Salvetti correlation-energy
c          formula into a functional of the electron density"
c         Phys.Rev. Vol. B37 (1988) 785-789.
c
c     [2] B. Miehlich, A. Savin, H. Stoll, and H. Preuss,
c         "Results obtained with the correlation energy density
c          functionals of Becke and Lee, Yang and Parr"
c         Chem.Phys.Lett. Vol. 157 (1989) 200-206.
c
c     Parameters:
c
c     ra     the alpha-electron density
c     rb     the beta-electron density
c     gaa    the dot product of the alpha density gradient with itself
c     gbb    the dot product of the beta density gradient with itself
c     gab    the dot product of the alpha density gradient with
c            the beta density
c     f      On return the functional value
c     dfdra  On return the derivative of f with respect to ra
c     dfdrb  On return the derivative of f with respect to rb
c     dfdgaa On return the derivative of f with respect to gaa
c     dfdgbb On return the derivative of f with respect to gbb
c     dfdgab On return the derivative of f with respect to gab
*/
//      double ra, rb, gaa, gbb, gab;
//      double f, dfdra, dfdrb, dfdgaa, dfdgbb, dfdgab;
//c
      double rhoa13,rhob13,rho13,rho43;
      double rhoa83,rhob83;
      double r_a,r_b,rab,raa,rbb,rab9;
      double faa,fbb,fab,dgaa,dgbb,dgab;
      double d2agaa,d2bgaa;
      double d2agbb,d2bgbb;
      double d2agab,d2bgab;
      double dfaaa,dfaab;
      double dfbba,dfbbb;
      double dfaba,dfabb;
      double del,delp,ome,ro,xtc,xtd,rx,lt1,lt2a,lt2b,lt2;
      double dlt1a,dlt1b,dlt2a,dlt2b;
      double zeta, cd ;
//c
//      double t1,ta,tb,tc,td,t2,t3,t43,t83,t19,n13;
      const double n13 = 0.3333333333333333;
      const double t1  = 0.3646239897876487e02;     //   ! p1
      const double ta  = 0.04918;       //                ! a
      const double tb  = 0.13200;     //                ! b
      const double tc  = 0.25330;     //                ! c
      const double td  = 0.34900;     //                ! d
      const double t2  = ta*tb/td ;     //                ! p2
      const double t3  = 4.0*ta/td;     //              ! p3
      const double t43 = 4.0*n13;
      const double t83 = 8.0*n13;
      const double t19 = 1.0/9.0;
//c
      cd     = ra + rb;
      rho13  = pow(cd,n13);
      rhoa13 = pow(ra,n13);
      rhob13 = pow(rb,n13);
//c
      rho43  = cd*rho13;
      rhoa83 = rhoa13*rhoa13*ra*ra;
      rhob83 = rhob13*rhob13*rb*rb;
//c
      r_a  = ra/cd;
      r_b  = rb/cd;
//c
      rab = r_a*r_b;
      raa = r_a*r_a;
      rbb = r_b*r_b;
//c
      zeta = r_a - r_b;
//c
      xtc = tc/rho13;
      xtd = td/rho13;
//c
      rx = xtd/(1. + xtd);
//c
      del = xtc + rx;
      ome = t2*rx*exp(-xtc)/rho43 ;
//c
      rab9 = t19*rab;
//c
      faa = rab9*(1. - 3.*del - (del - 11.)*r_a);
      fbb = rab9*(1. - 3.*del - (del - 11.)*r_b);
      fab = rab9*(47. - 7.*del);
//c
      dgaa = -ome*(faa - rbb);
      dgbb = -ome*(fbb - raa);
      dgab = -ome*(fab - t43);
//c
      lt1  = -t3*rx*rab*rho43;
      lt2a = -(t1*ome*rab)*rhoa83;
      lt2b = -(t1*ome*rab)*rhob83 ;
      lt2  = lt2a + lt2b;
//c
      f = lt1 + lt2 + dgaa*gaa + dgbb*gbb + dgab*gab;
//c
      delp = n13*(rx*rx - del)/cd;
      ro   = n13*(del - 5.)/cd;
//c
      if (ra <= 0.0 || rb <= 0.0){
         dfaaa = 0.0;
         dfaab = 0.0;
         dfbba = 0.0;
         dfbbb = 0.0;
         dfaba = 0.0;
         dfabb = 0.0;
         dlt1a = 0.0;
         dlt1b = 0.0;
         dlt2a = 0.0;
         dlt2b = 0.0;
      }else{
         dfaaa = -(zeta/ra)*faa - rab9*( (3. + r_a)*delp + (del - 11.)*r_b/cd );
         dfaab =  (zeta/rb)*faa - rab9*( (3. + r_a)*delp - (del - 11.)*r_a/cd );
         dfbba = -(zeta/ra)*fbb - rab9*( (3. + r_b)*delp - (del - 11.)*r_b/cd );
         dfbbb =  (zeta/rb)*fbb - rab9*( (3. + r_b)*delp + (del - 11.)*r_a/cd );
         dfaba = -(zeta/ra)*fab - 7.*rab9*delp;
         dfabb =  (zeta/rb)*fab - 7.*rab9*delp;
//c
         dlt1a = (lt1/cd)*(n13*rx + r_b/r_a);
         dlt1b = (lt1/cd)*(n13*rx + r_a/r_b);
         dlt2a = (ro - (zeta/ra))*lt2b +  (ro + (t83 - zeta)/ra)*lt2a;
         dlt2b = (ro + (zeta/rb))*lt2a +  (ro + (t83 + zeta)/rb)*lt2b;
}//      endif
//c
      d2agaa = ro*dgaa - ome * (dfaaa + rbb*(2./cd));
      d2bgaa = ro*dgaa - ome * (dfaab - rab*(2./cd));
      d2agbb = ro*dgbb - ome * (dfbba - rab*(2./cd));
      d2bgbb = ro*dgbb - ome * (dfbbb + raa*(2./cd));
      d2agab = ro*dgab - ome* dfaba;
      d2bgab = ro*dgab - ome* dfabb;
//c
      dfdra = dlt1a + dlt2a +  d2agaa*gaa + d2agbb*gbb + d2agab*gab;
      dfdrb = dlt1b + dlt2b +  d2bgaa*gaa + d2bgbb*gbb + d2bgab*gab;
//c
      dfdgaa = dgaa;
      dfdgbb = dgbb;
      dfdgab = dgab;
//c
     }// end
//c-----------------------------------------------------------------------
//c-----------------------------------------------------------------------
      void c_rks_vwnrpa(double r,double &f,double &dfdra){
//      implicit none
/*c
c     This subroutine evaluates the Vosko, Wilk and Nusair correlation
c     functional for the closed shell case using [4.5] with the RPA
c     parametrisation given below equation [4.4].
c
c     The original code was obtained from Dr. Phillip Young.
c
c     [1] S.H. Vosko, L. Wilk, and M. Nusair
c         "Accurate spin-dependent electron liquid correlation energies
c          for local spin density calculations: a critical analysis",
c         Can.J.Phys, Vol. 58 (1980) 1200-1211.
c
c     Parameters:
c
c     r      the total electron density
c     f      On return the functional value
c     dfdra  On return the derivative of f with respect to the alpha-
c            electron density
*/
//      double r;
//      double f, dfdra;
//c
      double t4,t5,t6,t7;
      double a2,b2,c2,d2;
      double p1,p2,p3,p4;
      double srho,srho13;
      double iv,iv2,inv,i1,i2,i3;
      double pp1,pp2;
//c
//      double n13, n43, n16, one, two, four;
      const double n13 = 0.333333333333333333333333333333;
      const double n43 = 1.333333333333333333333333333333;
      const double n16 = 0.166666666666666666666666666666;
      const double one = 1.0;
      const double two = 2.0;
      const double four= 4.0;
/*
C VWN (RPA) interpolation parameters
C
C paramagnetic
*/
      a2 =  0.0621814;
      b2 =  13.0720;    //       ! ch
      c2 =  42.7198;    //       ! ch
      d2 = -0.409286;   //       ! ch
//C
//C t4 = (1/(4/3)*pi)**(1/3)
      t4=0.620350490899399531;
//C
//C t5 = 0.5/(2**(1/3)-1)
      t5 = 1.92366105093153617;
//C
//C t6 = 2/(3*(2**(1/3)-1))
      t6 = 2.56488140124204822;
//C
//C t7 = 2.25*(2**(1/3)-1)
      t7 = 0.584822362263464735;
/*
C
C Paramagnetic interpolation constants
*/
      p1 = 0.448998886415767975e-01;
      p2 = 0.310907e-01;
      p3 = 0.443137364463981956e-02;
      p4 = 20.5219723675762680;
//C
//C closed shell case
      srho        = r;
      srho13      = pow(srho,n13);
      iv2         = t4/srho13;
      iv          = sqrt(iv2);
//C
//C paramagnetic
      inv = 1.0/(iv2+b2*iv+c2);
      i1  = log(iv2*inv);
      i2  = log((iv-d2)*(iv-d2)*inv);
//c corrected b1->b2 ps Apr98
      i3  = atan(p1/(2.0*iv+b2));
      pp1 = p2*i1 + p3*i2 + p4*i3;
      pp2 = a2*(1.0/iv-iv*inv*(1.0+b2/(iv-d2)));
//C
      f     = pp1*srho;
      dfdra = pp1 - n16*iv*pp2;
//c
// Added 15 July 2003 to stop infinities
if (r==0.0){
      f      = 0.0;
      dfdra  = 0.0;
 }


 }//    end
//-----------------------------------------------------------------------
      void c_uks_vwnrpa(double ra,double rb,double &f,double &dfdra,double &dfdrb){
/*      implicit none

c     This subroutine evaluates the Vosko, Wilk and Nusair correlation
c     functional using [4.5] with the RPA parametrisation given below
c     equation [4.4].
c
c     The original code was obtained from Dr. Phillip Young.
c
c     [1] S.H. Vosko, L. Wilk, and M. Nusair
c         "Accurate spin-dependent electron liquid correlation energies
c          for local spin density calculations: a critical analysis",
c         Can.J.Phys, Vol. 58 (1980) 1200-1211.
c
c     Parameters:
c
c     ra     the alpha-electron density
c     rb     the beta-electron density
c     f      On return the functional value
c     dfdra  On return the derivative of f with respect to ra
c     dfdrb  On return the derivative of f with respect to rb
*/
//      double ra, rb;
//      double f, dfdra, dfdrb;
//c
      double t1,t2,t4,t5,t6,t7;
      double a1,b1,c1,d1;
      double a2,b2,c2,d2;
      double a3,b3,c3,d3;
      double S1,S2,S3,S4;
      double p1,p2,p3,p4;
      double f1,f2,f3,f4;
      double inter1,inter2;
      double alpha_rho13,beta_rho13;
      double srho,srho13;
      double iv,iv2,inv,i1,i2,i3;
      double vwn1,vwn2;
      double ss1,ss2,pp1,pp2,ff1,ff2;
      double zeta,zeta3,zeta4;
      double tau,dtau,v;
//c
//      double n13, n43, n16, one, two, four, tn13, tn43
      const double n13 = 0.333333333333333333333333333333;
      const double n43 = 1.333333333333333333333333333333;
      const double n16 = 0.166666666666666666666666666666;
      const double tn13= 1.25992104989487316476721060727823;
      const double tn43= 2.51984209978974632953442121455646;
      const double one = 1.0;
      const double two = 2.0;
      const double four= 4.0;
/*
C VWN (RPA) interpolation parameters
C
C spin stiffness
*/
      a1 = -0.337737278807791058e-01;
      b1 =  1.13107;
      c1 = 13.0045;
      d1 = -0.00475840;
// paramagnetic
      a2 =  0.0621814;
      b2 =  13.0720;         //  ! ch
      c2 =  42.7198;         //  ! ch
      d2 = -0.409286;        //  ! ch
// ferromagnetic
// try cadpac/nwchem value (.5*a2)
      a3 = 0.0310907;
      b3 =  20.1231;         //  ! ch
      c3 = 101.578;          //  ! ch
      d3 = -0.743294;        //  ! ch
//
// t4 = (1/(4/3)*pi)**(1/3)
      t4=0.620350490899399531;
//
// t5 = 0.5/(2**(1/3)-1)
      t5 = 1.92366105093153617;
//
// t6 = 2/(3*(2**(1/3)-1))
      t6 = 2.56488140124204822;
//
// t7 = 2.25*(2**(1/3)-1)
      t7 = 0.584822362263464735;
/*
C Spin stiffness interpolation constants
*/
      S1 = 7.12310891781811772;
      S2 = a1 *0.5;
      S3 = -0.139834647015288626e-04 *0.5;
      S4 = -0.107301836977671539e-01 *0.5;
/*
C Paramagnetic interpolation constants
C
c     p1 = sqrt(4.0d0*c2 - b2*b2)
c     p2 = 0.5d0 * a2
c     p3 = -0.5d0 * a2*b2*d2/(d2*d2 + b2*d2 + c2)
c     p4 = 0.5d0 * a2*(two*b2/p1)*((c2-d2*d2)/(d2*d2 + b2*d2 + c2))
*/
      p1 = 0.448998886415767975e-01;
      p2 = 0.310907e-01;
      p3 = 0.443137364463981956e-02;
      p4 = 20.5219723675762680;
/*
C Ferromagnetic interpolation constants
C
c     f1 = sqrt(four*c3 - b3*b3)
c     f2 = 0.5d0 * a3
c     f3 = -0.5d0 * a3*b3*d3/(d3*d3 + b3*d3 + c3)
c     f4 = 0.5d0 * a3*(two*b3/f1)*((c3-d3*d3)/(d3*d3 + b3*d3 + c3))
*/
      f1 = 1.17168527770897146;
      f2 = 0.155453495681285858e-01;
      f3 = 0.266730993317173832e-02;
      f4 = 0.618818012598993383;
/*
C Interpolation intervals
*/
      inter1 =  1.0 -1.0e-10;
      inter2 = -1.0 +1.0e-10;
//
// open shell case
      alpha_rho13 = pow(ra,n13);
      beta_rho13  = pow(rb,n13);
      srho        = ra+rb;
      srho13      = pow(srho,n13);
      iv2         = t4/srho13;
      iv          = sqrt(iv2);
//
// spin-stiffness
      inv = 1.0/(iv2+b1*iv+c1);
      i1  = log(iv2*inv);
      i2  = log((iv-d1)*(iv-d1)*inv);
      i3  = atan(S1/(2.0*iv+b1));
      ss1 = S2*i1 + S3*i2 + S4*i3;
      ss2 = a1*(1.0/iv-iv*inv*(1.0+b1/(iv-d1)));
//
// paramagnetic
      inv = 1.0/(iv2+b2*iv+c2);
      i1  = log(iv2*inv);
      i2  = log((iv-d2)*(iv-d2)*inv);
// corrected b1->b2 ps Apr98
      i3  = atan(p1/(2.0*iv+b2));
      pp1 = p2*i1 + p3*i2 + p4*i3 ;
      pp2 = a2*(1.0/iv-iv*inv*(1.0+b2/(iv-d2))) ;
//
// ferromagnetic
      inv = 1.0/(iv2+b3*iv+c3) ;
      i1  = log(iv2*inv) ;
      i2  = log((iv-d3)*(iv-d3)*inv);
      i3  = atan(f1/(2.0*iv+b3));
      ff1 = f2*i1 + f3*i2 + f4*i3 ;
      ff2 = a3*(1.0/iv-iv*inv*(1.0+b3/(iv-d3)));
/*
C polarisation function
c     zeta  = (alpha_rho13-beta_rho13)/srho13
*/
      zeta  = (ra-rb)/srho;
      zeta3 = zeta*zeta*zeta;
      zeta4 = zeta3*zeta;
      if(zeta > inter1){
         vwn1 = (tn43-two)*t5 ;
         vwn2 = (tn13)*t6 ;
      }else if (zeta <  inter2){
         vwn1 = (tn43-two)*t5;
         vwn2 = -(tn13)*t6;
     }else{
         vwn1 = (pow((one+zeta),n43)+pow((one-zeta),n43)-two)*t5 ;
         vwn2 = (pow((one+zeta),n13)-pow((one-zeta),n13))*t6;
      }
      ss1  = ss1*t7;
      ss2  = ss2*t7;
//     tau  = ff1-pp1-ss1
//     dtau = ff2-pp2-ss2
      tau  = ff1-pp1;
      dtau = ff2-pp2;

//     v = pp1+vwn1*(ss1+tau*zeta4)
      v = pp1+vwn1*tau ;
      f = v*srho;

//     t1 = v - n16*iv*(pp2+vwn1*(ss2+dtau*zeta4))
//c     t2 = vwn2*(ss1+tau*zeta4)+vwn1*four*tau*zeta3
      t1 = v - n16*iv*(pp2+vwn1*dtau);
      t2 = vwn2*tau;
      dfdra = t1+t2*(one-zeta);
      dfdrb = t1-t2*(one+zeta);

}//    end
//c-----------------------------------------------------------------------

//c-----------------------------------------------------------------------
      void c_rks_pz81(double r,double &f,double &dfdra){
/*      implicit none
c
c     This subroutine evaluates the functional proposed by Perdew and
c     Zunger [1] for the closed shell case using the Ceperley-Alder
c     parameterization [2]. The corresponding potential is evaluated also.
c
c     [1] J.P. Perdew, and A. Zunger
c         "Self-interaction correction to density-functional
c          approximations for many-electron systems"
c         Physical Review B, Vol. 23 (1981) 5048-5079.
c         See appendix C on page 5075.
c
c     [2] D.M. Ceperley, and B.J. Alder,
c         Physical Review Letters Vol. 45 (1980) 566.
c
c
c     The original code was written by Dr. Alex de Vries, 1998. (Niet waar!!!)
c
c     Parameters:
c
c     r     the total electron density
c     f     On return the functional value
c     dfdra On return the derivative of f with respect to alpha-electron
c           density
c */
//      double r;
//      double f, dfdra;
//c
      
      double rhoval, rs, alnrs, eU, vU, drs, adU, ePol;
      double fz, zeta, eP, vP, adP, AAadd;
//c
//      double AU, BU, CU, DU, GU, B1U, B2U, AP, BP, CP, DP, GP;
//      double B1P, B2P, rsfact, ONE3, SEV6, FOUR3, TWO3;
      const double AU    =  0.03110;
      const double BU    = -0.0480;
      const double CU    =  0.00200;
      const double DU    = -0.01160;
      const double GU    = -0.1423;
      const double B1U   =  1.0529;
      const double B2U   =  0.3334;
      const double AP    =  0.01555;
      const double BP    = -0.0269;
      const double CP    =  0.0007;
      const double DP    = -0.0048;
      const double GP    = -0.0843;
      const double B1P   =  1.3981;
      const double B2P   =  0.2611;
//c     rsfact = (0.75/pi)**(1/3)
      const double rsfact=  0.620350490899400016668006812047778;
      const double ONE3  =  1.0/3.0;
      const double SEV6  =  7.0/6.0;
      const double FOUR3 =  4.0/3.0;
      const double TWO3  =  2.0/3.0;
/*c
c     real*8 x, fzeta, fzzprime
c     fzeta(x) = ((1.0d0+x)**FOUR3 +
c    &amp;            (1.0d0-x)**FOUR3 - 2.0d0) / (2.0d0**FOUR3-2.0d0)
c     fzzprime(x) = FOUR3*((1.0d0+x)**ONE3 -
c    &amp;                     (1.0d0-x)**ONE3) / (2.0d0**FOUR3-2.0d0)
c  */
    if (r> 1.0e-50){
      rhoval = r;
      zeta = 0.0;
//c     fz = fzeta(zeta)
      fz = 0.0;
      rs = pow(rsfact/rhoval,ONE3);
      if (rs < 1.0){
         alnrs = log(rs);
         eU = AU*alnrs+BU+CU*rs*alnrs+DU*rs;
         eP = AP*alnrs+BP+CP*rs*alnrs+DP*rs;
         vU = AU*alnrs + (BU-ONE3*AU) + TWO3*CU*rs*alnrs + ONE3*(DU+DU-CU)*rs;
         vP = AP*alnrs + (BP-ONE3*AP) + TWO3*CP*rs*alnrs + ONE3*(DP+DP-CP)*rs;
      }else{
         drs = sqrt(rs);
         adU = 1.0/(1.0 +B1U*drs+B2U*rs);
         adP = 1.0/(1.0 +B1P*drs+B2P*rs);
         eU  = GU*adU;
         eP  = GP*adP;
         vU  = eU*(1.0 +SEV6*B1U*drs+FOUR3*B2U*rs)*adU;
         vP  = eP*(1.0 +SEV6*B1P*drs+FOUR3*B2P*rs)*adP;
}//      endif
//c     AAadd = vU + fz*(vP-vU)
      AAadd = vU;
//c     ePOL = (eP-eU)*fzzprime(zeta)
      ePol = 0.0;
//c     f     = (eU+fz*(eP-eU))*rhoval
      f     = eU*rhoval;
//c     dfdra = (AAadd + (1.d0-zeta)*ePOL)
      dfdra = AAadd;
//c
  }else{
      f      = 0.0;
      dfdra  = 0.0;
 }


}//      end
//c-----------------------------------------------------------------------
      void c_uks_pz81(double ra,double rb,double &f,double &dfdra,double &dfdrb){
/*      implicit none
c
c     This subroutine evaluates the functional proposed by Perdew and
c     Zunger [1] using the Ceperley-Alder parameterization [2]. The
c     corresponding potential is evaluated also.
c
c     [1] J.P. Perdew, and A. Zunger
c         "Self-interaction correction to density-functional
c          approximations for many-electron systems"
c         Physical Review B, Vol. 23 (1981) 5048-5079.
c         See appendix C on page 5075.
c
c     [2] D.M. Ceperley, and B.J. Alder,
c         Physical Review Letters Vol. 45 (1980) 566.
c
c
c     The original code was written by Dr. Alex de Vries, 1998. (Niet waar!!!)
c
c     Parameters:
c
c     ra    the alpha electron density
c     rb    the beta  electron density
c     f     On return the functional value
c     dfdra On return the derivative of f with respect to ra
c     dfdrb On return the derivative of f with respect to rb
c  */
//      double ra, rb;
//      double f, dfdra, dfdrb;
//c
      double rhoval, rs, alnrs, eU, vU, drs, adU, ePol;
      double fz, zeta, eP, vP, adP, AAadd;
//c
//      double AU, BU, CU, DU, GU, B1U, B2U, AP, BP, CP, DP, GP;
//      double B1P, B2P, rsfact, ONE3, SEV6, FOUR3, TWO3;
      const double AU    =  0.0311;
      const double BU    = -0.048;
      const double CU    =  0.0020;
      const double DU    = -0.0116;
      const double GU    = -0.1423;
      const double B1U   =  1.0529;
      const double B2U   =  0.3334;
      const double AP    =  0.01555;
      const double BP    = -0.0269;
      const double CP    =  0.0007;
      const double DP    = -0.0048;
      const double GP    = -0.0843;
      const double B1P   =  1.3981;
      const double B2P   =  0.2611;
//c     rsfact = (0.75/pi)**(1/3)
      const double rsfact=  0.620350490899400016668006812047778;
      const double ONE3  =  1.0/3.0;
      const double SEV6  =  7.0/6.0;
      const double FOUR3 =  4.0/3.0;
      const double TWO3  =  2.0/3.0;
//c

//c
      rhoval = ra + rb;
      zeta = (ra-rb)/rhoval;
      fz = fzeta(zeta);
      rs = pow(rsfact/rhoval,ONE3);
      if (rs <1.0 ){
         alnrs = log(rs);
         eU = AU*alnrs+BU+CU*rs*alnrs+DU*rs;
         eP = AP*alnrs+BP+CP*rs*alnrs+DP*rs;
         vU = AU*alnrs + (BU-ONE3*AU) + TWO3*CU*rs*alnrs + ONE3*(DU+DU-CU)*rs;
         vP = AP*alnrs + (BP-ONE3*AP) + TWO3*CP*rs*alnrs +  ONE3*(DP+DP-CP)*rs;
      }else{
         drs = sqrt(rs);
         adU = 1.0/(1.0+B1U*drs+B2U*rs);
         adP = 1.0/(1.0+B1P*drs+B2P*rs);
         eU  = GU*adU ;
         eP  = GP*adP;
         vU  = eU*(1.0+SEV6*B1U*drs+FOUR3*B2U*rs)*adU;
         vP  = eP*(1.0+SEV6*B1P*drs+FOUR3*B2P*rs)*adP;
}//      endif
      AAadd = vU + fz*(vP-vU);
      ePol = (eP-eU)*fzzprime(zeta);
      f     = (eU+fz*(eP-eU))*rhoval ;
      dfdra = (AAadd + (1.0-zeta)*ePol);
      dfdrb = (AAadd - (1.0+zeta)*ePol);
//c
}//      end


//      double x, fzeta, fzzprime;
      inline double fzeta(double x){
      const double FOUR3 =  4.0/3.0;  
      return (pow((1.0+x),FOUR3) + pow((1.0-x),FOUR3) - 2.0) / (pow(2.0,FOUR3)-2.0);
      }

      inline double fzzprime(double x){
      const double ONE3  =  1.0/3.0;
      const double FOUR3 =  4.0/3.0;        
      return FOUR3*(pow((1.0+x),ONE3) - pow((1.0-x),ONE3)) / (pow(2.0,FOUR3)-2.0);
     }


//c-----------------------------------------------------------------------
      void x_rks_pbe(double r,double g,double &f,double &dfdra,double &dfdgaa){
/*      implicit none
c
c     This subroutine evaluates the PBE exchange functional [1] for
c     the closed shell case. The original PBE expression requires
c     second derivatives of the density to evaluate the potential. This
c     requirement was lifted by Matt Challacombe.
c
c     The original subroutine was provided by
c
c        Matt Challacombe
c        Los Alamos National Laboratory
c        The University of California
c
c     and is part of the MondoSCF suite of linear scaling electronic
c     structure codes. It was modified to conform with the design of
c     repository.
c
c     [1] John P. Perdew, Kieron Burke, Matthias Ernzerhof
c         "Generalized gradient approximation made simple"
c         Physical review letters Vol. 77 (1996) 3865-3868.
c
c     Parameters:
c
c     r      the total electron density
c     g      the dot product of the total density gradient with itself
c     f      On return the functional value
c     dfdra  On return the derivative of f with respect to the alpha
c            electron density
c     dfdgaa On return the derivative of f with respect to the
c            dot product of the alpha density gradient with itself
c     */
   //   double r, g;
   //   double f, dfdra, dfdgaa;
//c
      double A, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11;
//c
      if (1.0e-50 < r){
         A=g;
         r1=pow(r,(1.0/3.0));
         r2=r*r1;
         r3=pow(r,2);
         r4=pow(r1,2);
         r5=r3*r4;
         r6=7.13182658760049e-3*A;
         r7=r5+r6;
         r8=pow(r7,2);
         r9=1.0/r8;
         r10=pow(A,2);
         r11=pow(r3,2);
         f=-1.332360014553170*r2+(5.93801248171146e-1*r2)/(1.0+ (7.13182658760049e-3*A)/r5);
         dfdra=(-9.03570152478593e-5*r1*r10-8.39954475162682e-3*A*r*r3-9.84745021842697e-1*r*r11*r4)*r9;
         dfdgaa=-4.23488752945734e-3*r11*r9;
         dfdgaa= 0.50*dfdgaa;
     }else{
         f      = 0.00;
         dfdra  = 0.00;
         dfdgaa = 0.00;
}//      ENDIF
}//      END
//c-----------------------------------------------------------------------
      void x_uks_pbe(double ra,double rb,double gaa,double gbb,double &f,double &dfdra,double &dfdrb,double &dfdgaa,double &dfdgbb){
/*      implicit none
c
c     This subroutine evaluates the PBE exchange functional [1] for
c     the open shell case. The original PBE expression requires
c     second derivatives of the density to evaluate the potential. This
c     requirement was lifted by Matt Challacombe. This subroutine
c     uses x_rks_pbe to evaluate the functional using the spin-scaling
c     relationship of the exact exchange energy (e.g. see [1] Eq. (11)),
c     in analogy to the original PBE code [2].
c
c     The original subroutine was provided by
c
c        Matt Challacombe
c        Los Alamos National Laboratory
c        The University of California
c
c     and is part of the MondoSCF suite of linear scaling electronic
c     structure codes. It was modified to conform with the design of
c     repository.
c
c     [1] John P. Perdew, Kieron Burke, Matthias Ernzerhof
c         "Generalized gradient approximation made simple"
c         Physical review letters Vol. 77 (1996) 3865-3868.
c
c     [2] John P. Perdew, Kieron Burke, Matthias Ernzerhof
c         "PBE alpha2.1"
c         http://crab.rutgers.edu/~kieron/pubs/PBE.asc
c         September 20, 1996.
c
c     Parameters:
c
c     ra     the alpha electron density
c     rb     the beta  electron density
c     gaa    the dot product of the alpha density gradient with itself
c     gbb    the dot product of the beta  density gradient with itself
c     f      On return the functional value
c     dfdra  On return the derivative of f with respect to the alpha
c            electron density
c     dfdrb  On return the derivative of f with respect to the beta
c            electron density
c     dfdgaa On return the derivative of f with respect to the
c            dot product of the alpha density gradient with itself
c     dfdgbb On return the derivative of f with respect to the
c            dot product of the beta  density gradient with itself
c   */
//      double ra, rb, gaa, gbb;
//      double f, dfdra, dfdrb, dfdgaa, dfdgbb;
//c
      double ra2, rb2, ga2, gb2, fa, fb;
//c
      if (1.0e-20 < ra){
         ra2 = 2.00*ra;
         ga2 = 4.00*gaa;
         x_rks_pbe(ra2,ga2,fa,dfdra,dfdgaa);
      }else{
         fa     = 0.00;
         dfdra  = 0.00;
         dfdgaa = 0.00;
     }// endif
      if (1.0e-20 < rb){
         rb2 = 2.00*rb;
         gb2 = 4.00*gbb;
         x_rks_pbe(rb2,gb2,fb,dfdrb,dfdgbb);
      }else{
         fb     = 0.00;
         dfdrb  = 0.00;
         dfdgbb = 0.00;
      }
      f = 0.50*(fa+fb);
}//      END
//c-----------------------------------------------------------------------