diff --git a/src/PBEFunctional.C b/src/PBEFunctional.C
index b482c991..ac1b141b 100644
--- a/src/PBEFunctional.C
+++ b/src/PBEFunctional.C
@@ -23,8 +23,11 @@
 #include <vector>
 using namespace std;
 
-PBEFunctional::PBEFunctional(const vector<vector<double> > &rhoe,
+PBEFunctional::PBEFunctional(const vector<vector<double> > &rhoe, bool sol,
   double x_coeff, double c_coeff)
+: sol(sol),
+  um(sol ? 10./81 : 0.2195149727645171),
+  bet(sol ? 0.046 : 0.06672455060314922)
 {
   x_coeff_ = x_coeff;
   c_coeff_ = c_coeff;
@@ -171,7 +174,6 @@ void PBEFunctional::excpbe(double rho, double grad,
   const double third  = 1.0 / 3.0;
   const double third4 = 4.0 / 3.0;
   const double ax = -0.7385587663820224058; /* -0.75*pow(3.0/pi,third) */
-  const double um = 0.2195149727645171;
   const double uk = 0.804;
   const double ul = um / uk;
   const double pi32third = 3.09366772628014; /* (3*pi^2 ) ^(1/3) */
@@ -179,7 +181,6 @@ void PBEFunctional::excpbe(double rho, double grad,
   const double seven_sixth  =  7.0 / 6.0;
   const double four_over_pi = 1.27323954473516;
   const double gamma = 0.03109069086965489; /* gamma = (1-ln2)/pi^2 */
-  const double bet = 0.06672455060314922; /* see [a] (4) */
   const double delt = bet / gamma;
 
   double rtrs,fk,twoks,rs,t,h,
@@ -290,7 +291,6 @@ void PBEFunctional::excpbe_sp(double rho_up, double rho_dn,
   const double third4 =  4.0 / 3.0;
   const double sixthm = -1.0 / 6.0;
   const double ax = -0.7385587663820224058; /* -0.75*pow(3.0/pi,third) */
-  const double um = 0.2195149727645171;
   const double uk = 0.804;
   const double ul = um / uk;
   const double pi32third = 3.09366772628014; /* (3*pi^2 ) ^(1/3) */
@@ -300,7 +300,6 @@ void PBEFunctional::excpbe_sp(double rho_up, double rho_dn,
   const double gam = 0.5198420997897463; /* gam = 2^(4/3) - 2 */
   const double fzz = 8.0 / ( 9.0 * gam );
   const double gamma = 0.03109069086965489; /* gamma = (1-ln2)/pi^2 */
-  const double bet = 0.06672455060314922; /* see [a] (4) */
   const double delt = bet / gamma;
   const double eta = 1.e-12; // small number to avoid blowup as |zeta|->1
 
diff --git a/src/PBEFunctional.h b/src/PBEFunctional.h
index 8116706b..a467ba74 100644
--- a/src/PBEFunctional.h
+++ b/src/PBEFunctional.h
@@ -24,6 +24,8 @@
 
 class PBEFunctional : public XCFunctional
 {
+  const bool sol; //!< PBEsol, if true
+  const double um, bet; //!< only constants that change between PBE and PBEsol
   double x_coeff_, c_coeff_;
   std::vector<double> _exc, _exc_up, _exc_dn;
   std::vector<double> _vxc1, _vxc1_up, _vxc1_dn,
@@ -46,12 +48,12 @@ class PBEFunctional : public XCFunctional
   public:
 
   // constructor with variable coefficients for exchange and correlation
-  // with default values 1.0
-  PBEFunctional(const std::vector<std::vector<double> > &rhoe,
+  // with default values 1.0, and sol = true/false switches between PBE and PBEsol
+  PBEFunctional(const std::vector<std::vector<double> > &rhoe, bool sol,
                 double x_coeff=1.0, double c_coeff=1.0);
 
   bool isGGA() const { return true; };
-  std::string name() const { return "PBE"; };
+  std::string name() const { return (sol ? "PBEsol" : "PBE"); };
   void setxc(void);
 };
 #endif
diff --git a/src/XCOperator.C b/src/XCOperator.C
index f43dfe83..2973a413 100644
--- a/src/XCOperator.C
+++ b/src/XCOperator.C
@@ -40,6 +40,7 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
   if ( ( functional_name ==  "LDA" ) ||
        ( functional_name ==  "VWN" ) ||
        ( functional_name ==  "PBE" ) ||
+       ( functional_name == "PBEsol" ) ||
        ( functional_name == "BLYP" ) )
   {
     // create only an xc potential
diff --git a/src/XCPotential.C b/src/XCPotential.C
index 045aa9f8..931f446e 100644
--- a/src/XCPotential.C
+++ b/src/XCPotential.C
@@ -48,7 +48,11 @@ XCPotential::XCPotential(const ChargeDensity& cd, const string functional_name,
   }
   else if ( functional_name == "PBE" )
   {
-    xcf_ = new PBEFunctional(rhototal_r_);
+    xcf_ = new PBEFunctional(rhototal_r_, false);
+  }
+  else if ( functional_name == "PBEsol" )
+  {
+    xcf_ = new PBEFunctional(rhototal_r_, true);
   }
   else if ( functional_name == "BLYP" )
   {
@@ -58,7 +62,7 @@ XCPotential::XCPotential(const ChargeDensity& cd, const string functional_name,
   {
     const double x_coeff = 1.0 - ctrl.alpha_PBE0;
     const double c_coeff = 1.0;
-    xcf_ = new PBEFunctional(rhototal_r_,x_coeff,c_coeff);
+    xcf_ = new PBEFunctional(rhototal_r_,false,x_coeff,c_coeff);
   }
   else if ( functional_name == "HSE" )
   {
diff --git a/src/Xc.h b/src/Xc.h
index de79ae63..edd134e7 100644
--- a/src/Xc.h
+++ b/src/Xc.h
@@ -47,6 +47,7 @@ class Xc : public Var
     if ( !( v == "LDA" ||
             v == "VWN" ||
             v == "PBE" ||
+            v == "PBEsol" ||
             v == "BLYP" ||
             v == "HF" ||
             v == "PBE0" ||
@@ -56,7 +57,7 @@ class Xc : public Var
             v == "BHandHLYP" ) )
     {
       if ( ui->onpe0() )
-        cout << " xc must be LDA, VWN, PBE, BLYP, "
+        cout << " xc must be LDA, VWN, PBE, PBEsol, BLYP, "
              << "HF, PBE0, HSE, RSH, B3LYP or BHandHLYP" << endl;
       return 1;
     }
diff --git a/src/testXCFunctional.C b/src/testXCFunctional.C
index b390ba62..da0edd75 100644
--- a/src/testXCFunctional.C
+++ b/src/testXCFunctional.C
@@ -51,7 +51,7 @@ int main(int argc, char **argv)
 
   XCFunctional *xcf_list[2];
   xcf_list[0] = new LDAFunctional(rh);
-  xcf_list[1] = new PBEFunctional(rh);
+  xcf_list[1] = new PBEFunctional(rh, false);
 
   for ( int ixcf = 0; ixcf < 2; ixcf++ )
   {