First RSH implementation

RSH is implemented as a separate functional (so far).
New variables alpha_RSH, beta_RSH and mu_RSH are added. There is
still duplicated code with HSEFunctional and PBEFunctional.
The InteractionPotential class is removed and ExchangeOperator is
modified to implement the interaction potential depending on the
parameters alpha_RSH, beta_RSH and mu_RSH.
Functionality needs testing for HF and all hybrid functionals.

src/AlphaRSH.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// Copyright (c) 2014 The Regents of the University of California
+//
+// This file is part of Qbox
+//
+// Qbox is distributed under the terms of the GNU General Public License
+// as published by the Free Software Foundation, either version 2 of
+// the License, or (at your option) any later version.
+// See the file COPYING in the root directory of this distribution
+// or <http://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// AlphaRSH.h
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALPHARSH_H
+#define ALPHARSH_H
+
+#include<iostream>
+#include<iomanip>
+#include<sstream>
+#include<stdlib.h>
+
+#include "Sample.h"
+
+class AlphaRSH : public Var
+{
+  Sample *s;
+
+  public:
+
+  const char *name ( void ) const { return "alpha_RSH"; };
+
+  int set ( int argc, char **argv )
+  {
+    if ( argc != 2 )
+    {
+      if ( ui->onpe0() )
+      cout << " alpha_RSH takes only one value" << endl;
+      return 1;
+    }
+
+    double v = atof(argv[1]);
+    if ( v < 0.0 )
+    {
+      if ( ui->onpe0() )
+        cout << " alpha_RSH must be non-negative" << endl;
+      return 1;
+    }
+
+    s->ctrl.alpha_RSH = v;
+
+    return 0;
+  }
+
+  string print (void) const
+  {
+     ostringstream st;
+     st.setf(ios::left,ios::adjustfield);
+     st << setw(10) << name() << " = ";
+     st.setf(ios::right,ios::adjustfield);
+     st << s->ctrl.alpha_RSH;
+     return st.str();
+  }
+
+  AlphaRSH(Sample *sample) : s(sample)
+  {
+    s->ctrl.alpha_RSH = 0;
+  }
+};
+#endif

src/BetaRSH.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// Copyright (c) 2014 The Regents of the University of California
+//
+// This file is part of Qbox
+//
+// Qbox is distributed under the terms of the GNU General Public License
+// as published by the Free Software Foundation, either version 2 of
+// the License, or (at your option) any later version.
+// See the file COPYING in the root directory of this distribution
+// or <http://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// BetaRSH.h
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#ifndef BETARSH_H
+#define BETARSH_H
+
+#include<iostream>
+#include<iomanip>
+#include<sstream>
+#include<stdlib.h>
+
+#include "Sample.h"
+
+class BetaRSH : public Var
+{
+  Sample *s;
+
+  public:
+
+  const char *name ( void ) const { return "beta_RSH"; };
+
+  int set ( int argc, char **argv )
+  {
+    if ( argc != 2 )
+    {
+      if ( ui->onpe0() )
+      cout << " beta_RSH takes only one value" << endl;
+      return 1;
+    }
+
+    double v = atof(argv[1]);
+    if ( v < 0.0 )
+    {
+      if ( ui->onpe0() )
+        cout << " beta_RSH must be non-negative" << endl;
+      return 1;
+    }
+
+    s->ctrl.beta_RSH = v;
+
+    return 0;
+  }
+
+  string print (void) const
+  {
+     ostringstream st;
+     st.setf(ios::left,ios::adjustfield);
+     st << setw(10) << name() << " = ";
+     st.setf(ios::right,ios::adjustfield);
+     st << s->ctrl.beta_RSH;
+     return st.str();
+  }
+
+  BetaRSH(Sample *sample) : s(sample)
+  {
+    s->ctrl.beta_RSH = 0.25;
+  }
+};
+#endif

src/Control.h

@@ -54,6 +54,9 @@ struct Control
 
   std::string xc;
   double alpha_PBE0;
+  double alpha_RSH;
+  double beta_RSH;
+  double mu_RSH;
   std::string spin;
   int delta_spin;
 

src/ExchangeOperator.C

@@ -35,11 +35,10 @@ using namespace std;
 #define Tag_States 5
 
 ////////////////////////////////////////////////////////////////////////////////
-ExchangeOperator::ExchangeOperator( Sample& s, double HFCoeff,
-  const InteractionPotential& interaction_potential ) :
-  s_(s), wf0_(s.wf), dwf0_(s.wf), wfc_(s.wf), HFCoeff_(HFCoeff),
-  interaction_potential_(interaction_potential),
-  coulomb_(interaction_potential.coulomb()),
+ExchangeOperator::ExchangeOperator( Sample& s, double alpha_sx,
+  double beta_sx, double mu_sx ) :
+  s_(s), wf0_(s.wf), dwf0_(s.wf), wfc_(s.wf),
+  coulomb_(alpha_sx==beta_sx),
   gcontext_(s.wf.sd(0,0)->context())
 {
   eex_ = 0.0; // exchange energy
@@ -541,9 +540,9 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
             else
             {
               int_pot1_[ig] = ( qpG21_[ig] > 0.0 ) ?
-                  interaction_potential_(qpG21_[ig]) : 0;
+                  vint(qpG21_[ig]) : 0;
               int_pot2_[ig] = ( qpG22_[ig] > 0.0 ) ?
-                  interaction_potential_(qpG22_[ig]) : 0;
+                  vint(qpG22_[ig]) : 0;
             }
 
             // if iKpi=0 (first k point)
@@ -658,7 +657,6 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
                     // acumulate weighted contributions
                     // Psi_j,kj(r) * TF( rhog[G]/|q+G|^2 ) and symmetric
                     // in dpsi_i. We take now into account the mixing coeff
-                    weight *= HFCoeff_;
                     for ( int ir = 0; ir < np012loc_; ir++ )
                     {
                       dstatei_[i][ir] += ( statej_[j][ir] * rhor1_[ir] +
@@ -713,8 +711,6 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
                     //
                     // We take also into account the mixing coefficient.
                     //
-                    weighti *= HFCoeff_;
-                    weightj *= HFCoeff_;
                     for ( int ir = 0; ir < np012loc_; ir++ )
                     {
                       dstatei_[i][ir] += ( statej_[j][ir] * rhor1_[ir] +
@@ -803,8 +799,7 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
       }
 
       // divergence corrections
-      const double factor = (coulomb_) ? 1.0 :
-        interaction_potential_.divergence_scaling(rcut_);
+      const double factor = (coulomb_) ? beta_sx_ : vint_div_scal(rcut_);
       const double integ = 4.0 * M_PI * sqrt(M_PI) / ( 2.0 * rcut_ ) * factor;
       const double vbz = pow(2.0*M_PI,3.0) / omega;
 
@@ -852,7 +847,7 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
           complex<double> *pf1=dci.valptr(i*dci.mloc());
           complex<double> *pf2=&force_kpi_[i*dci.mloc()];
           for ( int j = 0; j < dci.mloc(); j++ )
-            pf1[j] += pf2[j] - ps[j] * div_corr * HFCoeff_;
+            pf1[j] += pf2[j] - ps[j] * div_corr;
         }
 
       } // for i
@@ -863,7 +858,6 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
   // reduce the total energy
   gcontext_.dsum(1, 1, &exchange_sum, 1);
   extot = exchange_sum;
-  extot *= HFCoeff_;
 
   tm.stop();
 #ifdef DEBUG
@@ -1324,7 +1318,7 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
       // factor 2.0: real basis
       const double tg2i = g2i[ig];
       // V(G) = 1/G^2 for Coulomb potential
-      const double int_pot = (coulomb_) ? tg2i : interaction_potential_(g2[ig]);
+      const double int_pot = (coulomb_) ? beta_sx_*tg2i : vint(g2[ig]);
       const double expG2 = exp( - rc2 * g2[ig] );
       double t = 2.0 * expG2 * int_pot;
       SumExpG2 += t;
@@ -1335,9 +1329,8 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
         const double tgy = g_y[ig];
         const double tgz = g_z[ig];
         // factor 2.0: derivative of G^2
-        const double fac = 2.0 * ( (coulomb_) ? t * ( rc2 + tg2i ) :
-          ( t * rc2 - 2.0 * expG2 *
-            interaction_potential_.derivative(g2[ig]) ) );
+        const double fac = 2.0 * ( (coulomb_) ? beta_sx_*t * ( rc2 + tg2i ) :
+          ( t * rc2 - 2.0 * expG2 * dvint(g2[ig]) ) );
         sigma_sumexp[0] += fac * tgx * tgx;
         sigma_sumexp[1] += fac * tgy * tgy;
         sigma_sumexp[2] += fac * tgz * tgz;
@@ -1596,8 +1589,7 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
             // and |rho2(G)|^2*V(|G+q2|) to the exchange energy.
             // factor 2.0: real basis
             const double tg2i = g2i[ig];
-            const double int_pot = ( coulomb_ ) ? tg2i :
-              interaction_potential_(g2[ig]);
+            const double int_pot = ( coulomb_ ) ? beta_sx_*tg2i : vint(g2[ig]);
             const double factor2 = 2.0;
             const double t1 = factor2 * norm(rhog1_[ig]) * int_pot;
             const double t2 = factor2 * norm(rhog2_[ig]) * int_pot;
@@ -1610,9 +1602,8 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
               const double tgy = g_y[ig];
               const double tgz = g_z[ig];
               // factor 2.0: derivative of 1/G^2
-              const double fac1 = 2.0 * ( coulomb_ ? (t1 * tg2i) :
-                -factor2 * norm(rhog1_[ig]) *
-                interaction_potential_.derivative(g2[ig]) );
+              const double fac1 = 2.0 * ( coulomb_ ? beta_sx_*(t1 * tg2i) :
+                -factor2 * norm(rhog1_[ig]) * dvint(g2[ig]) );
               sigma_sum_1[0] += fac1 * tgx * tgx;
               sigma_sum_1[1] += fac1 * tgy * tgy;
               sigma_sum_1[2] += fac1 * tgz * tgz;
@@ -1620,9 +1611,8 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
               sigma_sum_1[4] += fac1 * tgy * tgz;
               sigma_sum_1[5] += fac1 * tgz * tgx;
 
-              const double fac2 = 2.0 * ( coulomb_ ? (t2 * tg2i) :
-                -factor2 * norm(rhog2_[ig]) *
-                interaction_potential_.derivative(g2[ig]) );
+              const double fac2 = 2.0 * ( coulomb_ ? beta_sx_*(t2 * tg2i) :
+                -factor2 * norm(rhog2_[ig]) * dvint(g2[ig]) );
               sigma_sum_2[0] += fac2 * tgx * tgx;
               sigma_sum_2[1] += fac2 * tgy * tgy;
               sigma_sum_2[2] += fac2 * tgz * tgz;
@@ -1661,7 +1651,7 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
 
             if (dwf)
             {
-              const double weight = exfac * occ_kj_[j1] * HFCoeff_;
+              const double weight = exfac * occ_kj_[j1];
               double *dp = (double *) &dstatei_[i1][0];
               double *pj = (double *) &statej_[j1][0];
               double *pr = (double *) &rhor1_[0];
@@ -1683,8 +1673,8 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
 
             if (dwf)
             {
-              double weighti = exfac * occ_ki_[i1] * HFCoeff_;
-              double weightj = exfac * occ_kj_[j1] * HFCoeff_;
+              double weighti = exfac * occ_ki_[i1];
+              double weightj = exfac * occ_kj_[j1];
               double *dpi = (double *) &dstatei_[i1][0];
               double *dpj = (double *) &dstatej_[j1][0];
               double *pi = (double *) &statei_[i1][0];
@@ -1715,7 +1705,7 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
 
             if (dwf)
             {
-              const double weight = exfac * occ_kj_[j2] * HFCoeff_;
+              const double weight = exfac * occ_kj_[j2];
               double *dp = (double *) &dstatei_[i2][0];
               double *pj = (double *) &statej_[j2][0];
               double *pr = (double *) &rhor1_[0];
@@ -1739,8 +1729,8 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
 
             if (dwf)
             {
-              double weighti = exfac * occ_ki_[i2] * HFCoeff_;
-              double weightj = exfac * occ_kj_[j2] * HFCoeff_;
+              double weighti = exfac * occ_ki_[i2];
+              double weightj = exfac * occ_kj_[j2];
               double *dpi = (double *) &dstatei_[i2][0];
               double *dpj = (double *) &dstatej_[j2][0];
               double *pi = (double *) &statei_[i2][0];
@@ -1796,8 +1786,7 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
             // and |rho2(G)|^2*V(|G|) to the exchange energy.
             // factor 2.0: real basis
             const double tg2i = g2i[ig];
-            const double int_pot = ( coulomb_ ) ? tg2i :
-              interaction_potential_(g2[ig]);
+            const double int_pot = ( coulomb_ ) ? beta_sx_*tg2i : vint(g2[ig]);
             const double factor2 = 2.0;
             const double t1 = factor2 * norm(rhog1_[ig]) * int_pot;
             ex_sum_1 += t1;
@@ -1808,9 +1797,8 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
               const double tgy = g_y[ig];
               const double tgz = g_z[ig];
               // factor 2.0: derivative of 1/G^2
-              const double fac = 2.0 * ( coulomb_ ? (t1 * tg2i) :
-                -factor2 * norm(rhog1_[ig]) *
-                interaction_potential_.derivative(g2[ig]) );
+              const double fac = 2.0 * ( coulomb_ ? beta_sx_*(t1 * tg2i) :
+                -factor2 * norm(rhog1_[ig]) * dvint(g2[ig]) );
               sigma_sum_1[0] += fac * tgx * tgx;
               sigma_sum_1[1] += fac * tgy * tgy;
               sigma_sum_1[2] += fac * tgz * tgz;
@@ -1845,7 +1833,7 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
 
             if (dwf)
             {
-              const double weight = exfac * occ_kj_[j1] * HFCoeff_;
+              const double weight = exfac * occ_kj_[j1];
               double *dp = (double *) &dstatei_[i1][0];
               double *pj = (double *) &statej_[j1][0];
               double *pr = (double *) &rhor1_[0];
@@ -1867,8 +1855,8 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
 
             if (dwf)
             {
-              double weighti = exfac * occ_ki_[i1] * HFCoeff_;
-              double weightj = exfac * occ_kj_[j1] * HFCoeff_;
+              double weighti = exfac * occ_ki_[i1];
+              double weightj = exfac * occ_kj_[j1];
               double *dpi = (double *) &dstatei_[i1][0];
               double *dpj = (double *) &dstatej_[j1][0];
               double *pi = (double *) &statei_[i1][0];
@@ -2123,8 +2111,7 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
 
       // analytical part
       // scaling factor relative to Coulomb potential
-      const double factor = (coulomb_) ? 1.0 :
-        interaction_potential_.divergence_scaling(rcut_);
+      const double factor = (coulomb_) ? 1.0 : vint_div_scal(rcut_);
       const double integ = 4.0 * M_PI * sqrt(M_PI) / ( 2.0 * rcut_ ) * factor;
       const double vbz = pow(2.0*M_PI,3.0) / omega;
 
@@ -2147,7 +2134,7 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
         complex<double> *ps=c.valptr(i*c.mloc());
         complex<double> *pf=dc.valptr(i*dc.mloc());
         for ( int j = 0; j < dc.mloc(); j++ )
-          pf[j] -= ps[j] * div_corr * HFCoeff_;
+          pf[j] -= ps[j] * div_corr;
       }
 
     } // for i
@@ -2165,14 +2152,9 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
   gcontext_.dsum(1, 1, &exchange_sum, 1);
   extot = exchange_sum;
 
-  extot *= HFCoeff_;
-
   // accumulate stress tensor contributions
   gcontext_.dsum(6,1,&sigma_exhf_[0],6);
 
-  // scale stress tensor with HF coefficient
-  sigma_exhf_ *= HFCoeff_;
-
   tm.stop();
 
 #ifdef DEBUG
@@ -2442,3 +2424,91 @@ void ExchangeOperator::CompleteSendingOccupations(int iRotationStep)
     wait_send_occupations_=0;
   }
 }
+
+////////////////////////////////////////////////////////////////////////////////
+// interaction potential
+// vint(r) = alpha_sx_ * erf(mu_sx_*r)/r + beta_sx * erfc(mu_sx_*r)/r
+//         = beta_sx_ / r + (alpha_sx_ - beta_sx_) * erf(mu_sx_*r)/r
+// Fourier transform:
+// vint(g2) = ( beta_sx + (alpha_sx-beta_sx) * exp( -g2 / 4*mu^2 ) ) / g2
+double ExchangeOperator::vint(double g2)
+{
+  const double r1_4w2 = 0.25 / ( mu_sx_ * mu_sx_ );
+  const double x = g2 * r1_4w2;
+
+  if ( g2 == 0 )
+  {
+    // if alpha_sx_ == 0, finite limit for g2 = 0
+    if ( alpha_sx_ == 0.0 )
+      return beta_sx_ * r1_4w2;
+    else
+      return 0.0;
+  }
+  else if ( g2 < 1e-6 )
+  {
+    if ( alpha_sx_ == 0.0 )
+    {
+      // Taylor expansion near origin
+      return r1_4w2 * beta_sx_ * ( 1.0 - 0.5 * x );
+    }
+    else
+      return 0.0;
+  }
+  else
+  {
+    return ( beta_sx_ + ( alpha_sx_ - beta_sx_ ) * exp(-x) ) / g2;
+  }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Derivative of the interaction potential vint(g2) w.r.t g2
+// dvint(g2) = d (vint(g2)) / d g2
+// dvint(g2) = vint(g2) * ( -1/g2 - 1/(4*w^2) )
+//
+//             exp( -g2 / 4 w^2 )   V(g2)
+// dvint(g2) = ------------------ - -----
+//                  4 g2 w^2         g2
+double ExchangeOperator::dvint(double g2)
+{
+  const double r1_4w2 = 0.25 / ( mu_sx_ * mu_sx_ );
+  const double x = g2 * r1_4w2;
+  const double third = 1.0 / 3.0;
+
+  if ( g2 == 0 )
+  {
+    if ( alpha_sx_ == 0.0 )
+    {
+      // if alpha_sx_ == 0, finite limit for g2 = 0
+      return -0.5 * beta_sx_ * r1_4w2 * r1_4w2;
+    }
+    else
+      return 0.0;
+  }
+  else if ( g2 < 1e-6 )
+  {
+    // Taylor expansion near origin
+    return beta_sx_ * ( -0.5 + x * third ) * r1_4w2 * r1_4w2;
+  }
+  else
+  {
+    // exact derivative
+    return - ( beta_sx_ / g2 +
+      ( alpha_sx_ - beta_sx_ ) * exp(-x) * ( r1_4w2 + 1.0/g2 ) ) / g2;
+  }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// scaling of the divergence correction relative to the Coulomb potential
+//
+// integral( exp(-rcut^2 G^2) * V(G^2) )            x
+// ------------------------------------- = 1 - -------------
+//  integral( exp(-rcut^2 G^2) / G^2 )         sqrt(x^2 + 1)
+//
+// with x = 2 * rcut * mu_sx_
+double ExchangeOperator::vint_div_scal(double rc)
+{
+  //!! next lines only correct if alpha_sx_==0
+  assert(alpha_sx_==0.0);
+  const double x = 2.0 * rc * mu_sx_;
+  return 1 - x / sqrt(x * x + 1);
+}

src/ExchangeOperator.h

@@ -19,7 +19,6 @@
 #include "Sample.h"
 #include "SlaterDet.h"
 #include "FourierTransform.h"
-#include "InteractionPotential.h"
 
 #ifndef EXCHANGEOPERATOR_H
 #define EXCHANGEOPERATOR_H
@@ -34,8 +33,6 @@ class ExchangeOperator
   double eex_;
   // constant of support function for exchange integration
   double rcut_;
-  // mixing coefficient for exchange energy and dwf accumulation
-  double HFCoeff_;
 
   // HF stress tensor
   std::valarray<double> sigma_exhf_;
@@ -172,25 +169,26 @@ class ExchangeOperator
   vector<DoubleMatrix*> uc_;
   vector<long int> localization_;
 
-  // Fourier transform of interaction potential
-  const InteractionPotential interaction_potential_;
+  // screened interaction potential paramters
+  double alpha_sx_, beta_sx_, mu_sx_;
+  // interaction potential. g2 is the wave vector squared.
+  double vint(double g2);
+  // derivative of the interaction potential w.r.t. g2
+  double dvint(double g2);
+  double vint_div_scal(double rc);
 
-  // coulomb potential
+  // coulomb potential flag (true if alpha_sx==beta_sx)
   bool coulomb_;
 
   public:
 
   // constructor
-  ExchangeOperator(Sample& s_, double HFCoeff,
-    const InteractionPotential& interaction_potential = InteractionPotential());
+  // screened interaction potential: alpha*erf(mu*r)/r + beta*erfc(mu*r)/r
+  ExchangeOperator(Sample& s_, double alpha_sx, double beta_sx, double mu_sx);
 
   // destructor
   ~ExchangeOperator();
 
-  // parameters
-  void setmixCoeff(double value) { HFCoeff_ = value; };
-  double HFCoeff() { return HFCoeff_; };
-
   // exchange energy and forces computation
   double eex() { return eex_; };
   double update_energy(bool compute_stress);

src/HSEFunctional.C

@@ -42,7 +42,7 @@ const double A = 1.0161144, B = -0.37170836, C = -0.077215461, D = 0.57786348,
 
 // constructor
 HSEFunctional::HSEFunctional(const vector<vector<double> > &rhoe) :
-  x_coeff_(0.75), c_coeff_(1.0)
+  x_coeff_(0.75), c_coeff_(1.0), omega(0.11)
 {
   // nonmagnetic or magnetic
   _nspin = rhoe.size();
@@ -968,7 +968,7 @@ void HSEFunctional::setxc(void)
     }
   }
 }
-
+#if 0
 // evaluate fourier transform of nonlocal potential for given
 // input g2 = G^2
 // fourier transform of erfc ( w r ) / r
@@ -1041,3 +1041,4 @@ double HSEFunctional::divergence_scaling(const double& rcut)
   const double x = 2.0 * rcut * omega;
   return 1 - x / sqrt(x * x + 1);
 }
+#endif

src/HSEFunctional.h

@@ -27,7 +27,6 @@
 #define HSEFUNCTIONAL_H
 
 #include "XCFunctional.h"
-#include "InteractionPotential.h"
 #include <vector>
 
 class HSEFunctional : public XCFunctional
@@ -35,7 +34,7 @@ class HSEFunctional : public XCFunctional
   const double x_coeff_, c_coeff_;
 
   // screening parameter of the HSE functional
-  static const double omega = 0.11;
+  const double omega; // == 0.11
 
   // vectors common to all GGA exchange functionals
   std::vector<double> _exc, _exc_up, _exc_dn;
@@ -61,22 +60,5 @@ class HSEFunctional : public XCFunctional
   }
 
   void setxc(void);
-
-  // evaluate fourier transform of nonlocal potential for given G vector
-  // input g2 = G^2
-  static double interaction_potential(const double& g2);
-  // derivative of interaction potential
-  static double derivative_interaction_potential(const double& g2);
-  // scaling of the divergence correction relative to the Coulomb potential
-  static double divergence_scaling(const double& rcut);
-
-  // construct interaction potential class
-  static const InteractionPotential make_interaction_potential()
-  {
-    return InteractionPotential(&interaction_potential,
-      &derivative_interaction_potential,&divergence_scaling);
-  }
-
 };
-
 #endif

src/InteractionPotential.h

-////////////////////////////////////////////////////////////////////////////////
-//
-// Copyright (c) 2008 The Regents of the University of California
-//
-// This file is part of Qbox
-//
-// Qbox is distributed under the terms of the GNU General Public License
-// as published by the Free Software Foundation, either version 2 of
-// the License, or (at your option) any later version.
-// See the file COPYING in the root directory of this distribution
-// or <http://www.gnu.org/licenses/>.
-//
-////////////////////////////////////////////////////////////////////////////////
-//
-// InteractionPotential.h
-//
-////////////////////////////////////////////////////////////////////////////////
-//
-// Implements InteractionPotential class that evaluates the potential for given
-// norm of G vectors and the derivative w.r.t. this argument
-//
-////////////////////////////////////////////////////////////////////////////////
-
-#ifndef INTERACTIONPOTENTIAL_H
-#define INTERACTIONPOTENTIAL_H
-
-#include <cassert>
-
-class InteractionPotential
-{
-  public:
-
-  // default constructor = Coulomb potential
-  InteractionPotential() : coulomb_(true) {}
-
-  // constructor - define function and derivative
-  InteractionPotential(double(*V)(const double&), double(*dV)(const double&),
-    double(*div_scale)(const double&)) :
-    V_(V), dV_(dV), div_scale_(div_scale), coulomb_(false) {}
-
-  // is the interaction potential a coulomb potential?
-  inline bool coulomb() const
-  {
-    return coulomb_;
-  }
-
-  // evaluate the interaction potential for given norm of G vector
-  inline double operator()(const double G2) const
-  {
-    // the current implementation expects that the coulomb potential
-    // is treated externaly
-    assert(not coulomb_);
-    return V_(G2);
-  }
-
-  // evaluate the derivative of the interaction potential w.r.t. G^2
-  inline double derivative(const double G2) const
-  {
-    // the current implementation expects that the coulomb potential
-    // is treated externaly
-    assert(not coulomb_);
-    return dV_(G2);
-  }
-
-  inline double divergence_scaling(const double rcut) const
-  {
-    // the current implementation expects that the coulomb potential
-    // is treated externaly
-    assert(not coulomb_);
-    return div_scale_(rcut);
-  }
-
-  private:
-
-  const bool coulomb_;
-  double (*V_)(const double&);
-  double (*dV_)(const double&);
-  double (*div_scale_)(const double&);
-
-};
-
-#endif

src/Makefile

@@ -31,6 +31,7 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
 	XCPotential.o LDAFunctional.o VWNFunctional.o \
         PBEFunctional.o BLYPFunctional.o B3LYPFunctional.o  \
 	ExponentialIntegral.o HSEFunctional.o \
+        RSHFunctional.o \
 	NonLocalPotential.o SampleReader.o StructuredDocumentHandler.o \
 	SampleHandler.o AtomSetHandler.o WavefunctionHandler.o \
         XMLGFPreprocessor.o Base64Transcoder.o \
@@ -162,6 +163,9 @@ CXXFLAGS += -DTARGET='"$(TARGET)"'
 AlphaPBE0.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
 AlphaPBE0.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
 AlphaPBE0.o: Wavefunction.h Control.h
+AlphaRSH.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h UnitCell.h
+AlphaRSH.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
+AlphaRSH.o: Control.h
 AndersonMixer.o: AndersonMixer.h blas.h
 AngleCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
 AngleCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
@@ -225,6 +229,9 @@ Basis.o: Basis.h D3vector.h UnitCell.h
 Basis.o: D3vector.h UnitCell.h
 BasisMapping.o: Basis.h D3vector.h UnitCell.h Context.h blacs.h
 BasisMapping.o: BasisMapping.h
+BetaRSH.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h UnitCell.h
+BetaRSH.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
+BetaRSH.o: Control.h
 Bisection.o: Bisection.h Context.h blacs.h SlaterDet.h Basis.h D3vector.h
 Bisection.o: UnitCell.h Matrix.h Timer.h jade.h FourierTransform.h
 Bisection.o: Context.h blacs.h SlaterDet.h Basis.h D3vector.h UnitCell.h
@@ -291,7 +298,7 @@ CellStepper.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
 CellStepper.o: Wavefunction.h Control.h
 ChargeDensity.o: ChargeDensity.h Timer.h Context.h blacs.h Basis.h D3vector.h
 ChargeDensity.o: UnitCell.h Wavefunction.h FourierTransform.h SlaterDet.h
-ChargeDensity.o: Matrix.h
+ChargeDensity.o: Matrix.h blas.h
 ChargeDensity.o: Timer.h Context.h blacs.h
 ChargeMixCoeff.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
 ChargeMixCoeff.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
@@ -383,6 +390,7 @@ ExchangeOperator.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
 ExchangeOperator.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h
 ExchangeOperator.o: ExtForceSet.h Wavefunction.h Control.h SlaterDet.h
 ExchangeOperator.o: Basis.h Matrix.h Timer.h FourierTransform.h
+ExponentialIntegral.o: ExponentialIntegral.h
 ExtForce.o: ExtForce.h D3vector.h
 ExtForce.o: D3vector.h
 ExtForceCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
@@ -406,6 +414,8 @@ FourierTransform.o: Timer.h
 GlobalExtForce.o: GlobalExtForce.h ExtForce.h D3vector.h AtomSet.h Context.h
 GlobalExtForce.o: blacs.h Atom.h UnitCell.h D3tensor.h blas.h Species.h
 GlobalExtForce.o: ExtForce.h D3vector.h
+HSEFunctional.o: HSEFunctional.h XCFunctional.h ExponentialIntegral.h
+HSEFunctional.o: XCFunctional.h
 HelpCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
 HelpCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
 HelpCmd.o: ExtForceSet.h Wavefunction.h Control.h
@@ -466,6 +476,9 @@ Matrix.o: Context.h blacs.h
 MoveCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
 MoveCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
 MoveCmd.o: ExtForceSet.h Wavefunction.h Control.h
+MuRSH.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h UnitCell.h
+MuRSH.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
+MuRSH.o: Control.h
 Nempty.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h UnitCell.h
 Nempty.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
 Nempty.o: Control.h
@@ -534,6 +547,8 @@ PrintCmd.o: ExtForceSet.h Wavefunction.h Control.h
 QuitCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
 QuitCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
 QuitCmd.o: ExtForceSet.h Wavefunction.h Control.h
+RSHFunctional.o: RSHFunctional.h XCFunctional.h ExponentialIntegral.h
+RSHFunctional.o: XCFunctional.h
 RandomizeRCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
 RandomizeRCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
 RandomizeRCmd.o: ExtForceSet.h Wavefunction.h Control.h
@@ -702,14 +717,16 @@ XCOperator.o: XCOperator.h Sample.h AtomSet.h Context.h blacs.h Atom.h
 XCOperator.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
 XCOperator.o: ExtForceSet.h Wavefunction.h Control.h ChargeDensity.h Timer.h
 XCOperator.o: XCPotential.h ExchangeOperator.h SlaterDet.h Basis.h Matrix.h
-XCOperator.o: FourierTransform.h
+XCOperator.o: FourierTransform.h HSEFunctional.h XCFunctional.h
+XCOperator.o: RSHFunctional.h
 XCOperator.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
 XCOperator.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
 XCOperator.o: Wavefunction.h Control.h
 XCPotential.o: XCPotential.h Control.h D3vector.h ChargeDensity.h Timer.h
 XCPotential.o: Context.h blacs.h LDAFunctional.h XCFunctional.h
 XCPotential.o: VWNFunctional.h PBEFunctional.h BLYPFunctional.h
-XCPotential.o: B3LYPFunctional.h Basis.h UnitCell.h FourierTransform.h blas.h
+XCPotential.o: HSEFunctional.h RSHFunctional.h B3LYPFunctional.h Basis.h
+XCPotential.o: UnitCell.h FourierTransform.h blas.h
 XCPotential.o: Control.h D3vector.h ChargeDensity.h Timer.h Context.h blacs.h
 XMLGFPreprocessor.o: Timer.h Context.h blacs.h Base64Transcoder.h Matrix.h
 XMLGFPreprocessor.o: XMLGFPreprocessor.h
@@ -731,13 +748,13 @@ qb.o: LoadCmd.h MoveCmd.h PartialChargeCmd.h PlotCmd.h PrintCmd.h QuitCmd.h
 qb.o: RandomizeRCmd.h RandomizeVCmd.h RandomizeWfCmd.h ResetVcmCmd.h
 qb.o: RescaleVCmd.h RseedCmd.h RunCmd.h SaveCmd.h SetCmd.h SetVelocityCmd.h
 qb.o: SpeciesCmd.h StatusCmd.h StrainCmd.h TorsionCmd.h BisectionCmd.h
-qb.o: Bisection.h SlaterDet.h Basis.h Matrix.h AlphaPBE0.h AtomsDyn.h BlHF.h
-qb.o: BtHF.h Cell.h CellDyn.h CellLock.h CellMass.h ChargeMixCoeff.h
-qb.o: ChargeMixNdim.h ChargeMixRcut.h Debug.h Dspin.h Ecut.h Ecutprec.h
-qb.o: Ecuts.h Efield.h Polarization.h Emass.h ExtStress.h FermiTemp.h
-qb.o: IterCmd.h IterCmdPeriod.h Dt.h Nempty.h NetCharge.h Nrowmax.h Nspin.h
-qb.o: RefCell.h ScfTol.h Stress.h Thermostat.h ThTemp.h ThTime.h ThWidth.h
-qb.o: WfDiag.h WfDyn.h Xc.h
+qb.o: Bisection.h SlaterDet.h Basis.h Matrix.h AlphaPBE0.h AlphaRSH.h
+qb.o: AtomsDyn.h BetaRSH.h BlHF.h BtHF.h Cell.h CellDyn.h CellLock.h
+qb.o: CellMass.h ChargeMixCoeff.h ChargeMixNdim.h ChargeMixRcut.h Debug.h
+qb.o: Dspin.h Ecut.h Ecutprec.h Ecuts.h Efield.h Polarization.h Emass.h
+qb.o: ExtStress.h FermiTemp.h IterCmd.h IterCmdPeriod.h Dt.h MuRSH.h Nempty.h
+qb.o: NetCharge.h Nrowmax.h Nspin.h RefCell.h ScfTol.h Stress.h Thermostat.h
+qb.o: ThTemp.h ThTime.h ThWidth.h WfDiag.h WfDyn.h Xc.h
 qbox_xmlns.o: qbox_xmlns.h
 release.o: release.h
 sinft.o: sinft.h

src/MuRSH.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// Copyright (c) 2014 The Regents of the University of California
+//
+// This file is part of Qbox
+//
+// Qbox is distributed under the terms of the GNU General Public License
+// as published by the Free Software Foundation, either version 2 of
+// the License, or (at your option) any later version.
+// See the file COPYING in the root directory of this distribution
+// or <http://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// MuRSH.h
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#ifndef MURSH_H
+#define MURSH_H
+
+#include<iostream>
+#include<iomanip>
+#include<sstream>
+#include<stdlib.h>
+
+#include "Sample.h"
+
+class MuRSH : public Var
+{
+  Sample *s;
+
+  public:
+
+  const char *name ( void ) const { return "mu_RSH"; };
+
+  int set ( int argc, char **argv )
+  {
+    if ( argc != 2 )
+    {
+      if ( ui->onpe0() )
+      cout << " mu_RSH takes only one value" << endl;
+      return 1;
+    }
+
+    double v = atof(argv[1]);
+    if ( v < 0.0 )
+    {
+      if ( ui->onpe0() )
+        cout << " mu_RSH must be non-negative" << endl;
+      return 1;
+    }
+
+    s->ctrl.mu_RSH = v;
+
+    return 0;
+  }
+
+  string print (void) const
+  {
+     ostringstream st;
+     st.setf(ios::left,ios::adjustfield);
+     st << setw(10) << name() << " = ";
+     st.setf(ios::right,ios::adjustfield);
+     st << s->ctrl.mu_RSH;
+     return st.str();
+  }
+
+  MuRSH(Sample *sample) : s(sample)
+  {
+    s->ctrl.mu_RSH = 0.11;
+  }
+};
+#endif

src/RSHFunctional.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// Copyright (c) 2008 The Regents of the University of California
+//
+// This file is part of Qbox
+//
+// Qbox is distributed under the terms of the GNU General Public License
+// as published by the Free Software Foundation, either version 2 of
+// the License, or (at your option) any later version.
+// See the file COPYING in the root directory of this distribution
+// or <http://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// RSHFunctional.h
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// Range-separated hybrid functional (RSH)
+// J.H.Skone et al. Phys. Rev. B93, 235106 (2016)
+// RSH is defined by alpha_RSH, beta_RSH, mu_RSH
+// sigma = alpha_RSH * rho(r,r') * erf(r-r')/(r-r') +
+//         beta_RSH * rho(r,r') * erfc(r-r')/(r-r') +
+//         (1 - alpha_RSH) * Vx_LR(r,mu_RSH) +
+//         (1 - beta_RSH) * Vx_SR(r,mu_RSH)
+// The HSE functional is obtained using alpha_RSH=0, beta_RSH=0.25, mu_RSH=0.11
+// Heyd et al.,   J. Chem. Phys. 118, 8207 (2003)
+// Heyd et al.,   J. Chem. Phys. 120, 7274 (2004)
+// Krukau et al., J. Chem. Phys. 125, 224106 (2006)
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#ifndef RSHFUNCTIONAL_H
+#define RSHFUNCTIONAL_H
+
+#include "XCFunctional.h"
+#include <vector>
+
+class RSHFunctional : public XCFunctional
+{
+  const double x_coeff_, c_coeff_;
+  const double alpha_RSH_, beta_RSH_, mu_RSH_;
+  const double omega; // == mu_RSH
+
+  // vectors common to all GGA exchange functionals
+  std::vector<double> _exc, _exc_up, _exc_dn;
+  std::vector<double> _vxc1, _vxc1_up, _vxc1_dn, _vxc2, _vxc2_upup, _vxc2_updn,
+    _vxc2_dnup, _vxc2_dndn;
+  std::vector<double> _grad_rho[3], _grad_rho_up[3], _grad_rho_dn[3];
+
+  void RSH_exchange(const double rho, const double grad,
+  const double a_ex, const double w, double *ex, double *vx1, double *vx2);
+
+  void gcor2(double a, double a1, double b1, double b2,
+  double b3, double b4, double rtrs, double *gg, double *ggrs);
+
+  void PBE_correlation(const double rho, const double grad,
+  double *ec, double *vc1, double *vc2);
+
+  void PBE_correlation_sp(const double rho_up, const double rho_dn,
+  const double grad_up, const double grad_dn, const double grad, double *ec,
+  double *vc1_up, double *vc1_dn, double *vc2);
+
+  void approximateIntegral(const double omega_kF, const double Hs2,
+  const double D_term, const double dHs2_ds, double *appInt,
+  double *dAppInt_ds, double *dAppInt_dkF);
+
+  void RSH_enhance(const double s_inp, const double kF,
+  const double w, double *fx, double *dfx_ds, double* dfx_dkf);
+
+  public:
+
+  // constructor
+  RSHFunctional(const std::vector<std::vector<double> > &rhoe,
+                double alpha_RSH, double beta_RSH, double mu_RSH);
+
+  bool isGGA() const
+  {
+    return true;
+  }
+
+  // return the name of the functional
+  std::string name() const
+  {
+    return "RSH";
+  }
+
+  void setxc(void);
+};
+
+#endif

src/XCOperator.C

@@ -20,6 +20,7 @@
 #include "XCPotential.h"
 #include "ExchangeOperator.h"
 #include "HSEFunctional.h"
+#include "RSHFunctional.h"
 using namespace std;
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -46,16 +47,14 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
     hasPotential_ = true;
     hasGGA_ = xcp_->isGGA();
     hasHF_ = false;
-    HFmixCoeff_ = 0.0;
   }
   else if ( functional_name == "HF" )
   {
     // create exchange operator with mixing coeff=1
-    xop_ = new ExchangeOperator(s, 1.0);
+    xop_ = new ExchangeOperator(s, 1.0, 1.0, 0.0);
     hasPotential_ = false;
     hasGGA_ = false;
     hasHF_ = true;
-    HFmixCoeff_ = 1.0;
   }
   else if ( functional_name == "PBE0" )
   {
@@ -63,11 +62,10 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
     xcp_ = new XCPotential(cd, functional_name, s.ctrl);
 
     // create the exchange operator with mixing coeff=0.25
-    xop_ = new ExchangeOperator(s, s.ctrl.alpha_PBE0);
+    xop_ = new ExchangeOperator(s, s.ctrl.alpha_PBE0, s.ctrl.alpha_PBE0, 0.0);
     hasPotential_ = true;
     hasGGA_ = xcp_->isGGA();
     hasHF_ = true;
-    HFmixCoeff_ = s.ctrl.alpha_PBE0;;
   }
   else if ( functional_name == "HSE" )
   {
@@ -75,12 +73,22 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
     xcp_ = new XCPotential(cd, functional_name, s.ctrl);
 
     // create the exchange operator with mixing coeff=0.25
-    xop_ = new ExchangeOperator(s, 0.25,
-      HSEFunctional::make_interaction_potential() );
+    xop_ = new ExchangeOperator(s, 0.0, 0.25, 0.11);
+    hasPotential_ = true;
+    hasGGA_ = xcp_->isGGA();
+    hasHF_ = true;
+  }
+  else if ( functional_name == "RSH" )
+  {
+    // create an exchange potential
+    xcp_ = new XCPotential(cd, functional_name, s.ctrl);
+
+    // create the exchange operator with mixing coeff=beta_RSH
+    xop_ = new ExchangeOperator(s, s.ctrl.alpha_RSH, s.ctrl.beta_RSH,
+      s.ctrl.mu_RSH);
     hasPotential_ = true;
     hasGGA_ = xcp_->isGGA();
     hasHF_ = true;
-    HFmixCoeff_ = 0.25;
   }
   else if ( functional_name == "B3LYP" )
   {
@@ -88,11 +96,10 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
     xcp_ = new XCPotential(cd, functional_name, s.ctrl);
 
     // create the exchange operator with mixing coeff=0.20
-    xop_ = new ExchangeOperator(s, 0.20);
+    xop_ = new ExchangeOperator(s, 0.20, 0.20, 0.0);
     hasPotential_ = true;
     hasGGA_ = xcp_->isGGA();
     hasHF_ = true;
-    HFmixCoeff_ = 0.20;
   }
   else
   {

src/XCOperator.h

@@ -33,7 +33,6 @@ class XCOperator
   ExchangeOperator* xop_;
 
   const ChargeDensity& cd_;
-  double HFmixCoeff_ ;
   double exc_; // XC energy: includes local and HF terms
   double dxc_;
 

src/XCPotential.C

@@ -22,6 +22,7 @@
 #include "PBEFunctional.h"
 #include "BLYPFunctional.h"
 #include "HSEFunctional.h"
+#include "RSHFunctional.h"
 #include "B3LYPFunctional.h"
 #include "Basis.h"
 #include "FourierTransform.h"
@@ -59,6 +60,10 @@ XCPotential::XCPotential(const ChargeDensity& cd, const string functional_name,
   {
     xcf_ = new HSEFunctional(cd_.rhor);
   }
+  else if ( functional_name == "RSH" )
+  {
+    xcf_ = new RSHFunctional(cd_.rhor,ctrl.alpha_RSH,ctrl.beta_RSH,ctrl.mu_RSH);
+  }
   else if ( functional_name == "B3LYP" )
   {
     xcf_ = new B3LYPFunctional(cd_.rhor);

src/qb.C

@@ -77,7 +77,9 @@ using namespace std;
 #include "BisectionCmd.h"
 
 #include "AlphaPBE0.h"
+#include "AlphaRSH.h"
 #include "AtomsDyn.h"
+#include "BetaRSH.h"
 #include "BlHF.h"
 #include "BtHF.h"
 #include "Cell.h"
@@ -100,6 +102,7 @@ using namespace std;
 #include "IterCmd.h"
 #include "IterCmdPeriod.h"
 #include "Dt.h"
+#include "MuRSH.h"
 #include "Nempty.h"
 #include "NetCharge.h"
 #include "Nrowmax.h"
@@ -292,7 +295,9 @@ int main(int argc, char **argv, char **envp)
   ui.addCmd(new TorsionCmd(s));
 
   ui.addVar(new AlphaPBE0(s));
+  ui.addVar(new AlphaRSH(s));
   ui.addVar(new AtomsDyn(s));
+  ui.addVar(new BetaRSH(s));
   ui.addVar(new BlHF(s));
   ui.addVar(new BtHF(s));
   ui.addVar(new Cell(s));
@@ -314,6 +319,7 @@ int main(int argc, char **argv, char **envp)
   ui.addVar(new FermiTemp(s));
   ui.addVar(new IterCmd(s));
   ui.addVar(new IterCmdPeriod(s));
+  ui.addVar(new MuRSH(s));
   ui.addVar(new Nempty(s));
   ui.addVar(new NetCharge(s));
   ui.addVar(new Nrowmax(s));