Merge branch 'rsh-dev' into develop

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.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,23 @@ class ExchangeOperator
   vector<DoubleMatrix*> uc_;
   vector<long int> localization_;
 
-  // Fourier transform of interaction potential
-  const InteractionPotential interaction_potential_;
-
-  // coulomb potential
-  bool coulomb_;
+  // 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);
 
   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
@@ -378,12 +385,11 @@ ExchangeOperator.o: VectorLess.h ExchangeOperator.h Sample.h AtomSet.h
 ExchangeOperator.o: Context.h blacs.h Atom.h D3vector.h UnitCell.h D3tensor.h
 ExchangeOperator.o: blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
 ExchangeOperator.o: Control.h SlaterDet.h Basis.h Matrix.h Timer.h
-ExchangeOperator.o: FourierTransform.h InteractionPotential.h Bisection.h
+ExchangeOperator.o: FourierTransform.h Bisection.h
 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
-ExchangeOperator.o: InteractionPotential.h
 ExponentialIntegral.o: ExponentialIntegral.h
 ExtForce.o: ExtForce.h D3vector.h
 ExtForce.o: D3vector.h
@@ -408,9 +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 InteractionPotential.h
-HSEFunctional.o: ExponentialIntegral.h
-HSEFunctional.o: XCFunctional.h InteractionPotential.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
@@ -471,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
@@ -539,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
@@ -710,16 +720,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 InteractionPotential.h HSEFunctional.h
-XCOperator.o: XCFunctional.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: HSEFunctional.h InteractionPotential.h B3LYPFunctional.h
-XCPotential.o: 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
@@ -738,16 +748,16 @@ qb.o: Control.h Timer.h AngleCmd.h AtomCmd.h ComputeMLWFCmd.h MLWFTransform.h
 qb.o: BasisMapping.h ConstraintCmd.h DistanceCmd.h ExtForceCmd.h
 qb.o: FoldInWsCmd.h HelpCmd.h KpointCmd.h ListAtomsCmd.h ListSpeciesCmd.h
 qb.o: LoadCmd.h MoveCmd.h PartialChargeCmd.h PlotCmd.h PrintCmd.h QuitCmd.h
-qb.o: RandomizeRCmd.h RandomizeVCmd.h RandomizeWfCmd.h ResetRotationCmd.h
-qb.o: ResetVcmCmd.h RescaleVCmd.h RseedCmd.h RunCmd.h SaveCmd.h SetCmd.h
-qb.o: SetVelocityCmd.h SpeciesCmd.h StatusCmd.h StrainCmd.h TorsionCmd.h
-qb.o: BisectionCmd.h Bisection.h SlaterDet.h Basis.h Matrix.h AlphaPBE0.h
-qb.o: AtomsDyn.h BlHF.h BtHF.h Cell.h CellDyn.h CellLock.h CellMass.h
-qb.o: ChargeMixCoeff.h ChargeMixNdim.h ChargeMixRcut.h Debug.h Dspin.h Ecut.h
-qb.o: Ecutprec.h Ecuts.h Efield.h Polarization.h Emass.h ExtStress.h
-qb.o: FermiTemp.h IterCmd.h IterCmdPeriod.h Dt.h Nempty.h NetCharge.h
-qb.o: Nrowmax.h Nspin.h RefCell.h ScfTol.h Stress.h Thermostat.h ThTemp.h
-qb.o: ThTime.h ThWidth.h WfDiag.h WfDyn.h Xc.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 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/Xc.h

@@ -51,11 +51,12 @@ class Xc : public Var
             v == "HF" ||