functional.cc

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

          if (testgamma[j] > testrho[i]) continue;
          id.a.gamma = testgamma[j];
          id.compute_derived(0, need_density_gradient());
          ret += test(id);
        }
    }

  set_spin_polarized(1);
  ExEnv::out0() << "Testing with rho_a != rho_b" << endl;
  for (i=0; testrho[i] != -1.0; i++) {
      id.a.rho=testrho[i];
      for (j=0; testrho[j] != -1.0; j++) {
          id.b.rho=testrho[j];
          if (testrho[i]+testrho[j] == 0.0) continue;
          for (k=0; testgamma[k] != -1.0; k++) {
              if (testgamma[k] > testrho[i]) continue;
              id.a.gamma = testgamma[k];
              double sqrt_gamma_a = sqrt(id.a.gamma);
              for (l=0; testgamma[l] != -1.0; l++) {
                  if (testgamma[l] > testrho[j]) continue;
                  id.b.gamma = testgamma[l];
                  double sqrt_gamma_b = sqrt(id.b.gamma);
                  for (m=0; testgammaab[m] != -1.0; m++) {
                      // constrain gamma_ab to values allowed by the
                      // current gamma_a and gamma_b
                      id.gamma_ab = testgammaab[m];
                      if (id.gamma_ab > sqrt_gamma_a*sqrt_gamma_b) {
                          id.gamma_ab = sqrt_gamma_a*sqrt_gamma_b;
                        }
                      if (id.gamma_ab < -0.5*(id.a.gamma+id.b.gamma)) {
                          id.gamma_ab = -0.5*(id.a.gamma+id.b.gamma);
                        }
                      if (id.gamma_ab < -sqrt_gamma_a*sqrt_gamma_b) {
                          id.gamma_ab = -sqrt_gamma_a*sqrt_gamma_b;
                        }
                      id.compute_derived(1, need_density_gradient());
                      ret += test(id);
                    }
                }
            }
        }
    }
  return ret;
}

/////////////////////////////////////////////////////////////////////////////
// NElFunctional

static ClassDesc NElFunctional_cd(
  typeid(NElFunctional),"NElFunctional",1,"public DenFunctional",
  0, create<NElFunctional>, create<NElFunctional>);

NElFunctional::NElFunctional(StateIn& s):
  SavableState(s),
  DenFunctional(s)
{
}

NElFunctional::NElFunctional(const Ref<KeyVal>& keyval):
  DenFunctional(keyval)
{
}

NElFunctional::~NElFunctional()
{
}

void
NElFunctional::save_data_state(StateOut& s)
{
  DenFunctional::save_data_state(s);
}

void
NElFunctional::point(const PointInputData &id,
                     PointOutputData &od)
{
  od.zero();
  od.energy = id.a.rho + id.b.rho;
}

/////////////////////////////////////////////////////////////////////////////
// SumDenFunctional

static ClassDesc SumDenFunctional_cd(
  typeid(SumDenFunctional),"SumDenFunctional",1,"public DenFunctional",
  0, create<SumDenFunctional>, create<SumDenFunctional>);

SumDenFunctional::SumDenFunctional(StateIn& s):
  SavableState(s),
  DenFunctional(s),
  n_(0),
  funcs_(0),
  coefs_(0)
{
  s.get(n_);
  if (n_) {
      s.get(coefs_);
      funcs_ = new Ref<DenFunctional>[n_];
      for (int i=0; i < n_; i++)
          funcs_[i] << SavableState::restore_state(s);
    }
}

SumDenFunctional::SumDenFunctional() :
  n_(0),
  funcs_(0),
  coefs_(0)
{
}

SumDenFunctional::SumDenFunctional(const Ref<KeyVal>& keyval):
  DenFunctional(keyval),
  n_(0),
  funcs_(0),
  coefs_(0)
{
  int ncoef = keyval->count("coefs");
  int nfunc = keyval->count("funcs");
  if (ncoef != nfunc && ncoef != 0) {
      ExEnv::out0() << "SumDenFunctional: number of coefs and funcs differ" << endl;
      abort();
    }
  
  n_ = nfunc;
  coefs_ = new double[n_];
  funcs_ = new Ref<DenFunctional>[n_];
  for (int i=0; i < n_; i++) {
      if (ncoef)
          coefs_[i] = keyval->doublevalue("coefs", i);
      else
          coefs_[i] = 1.0;
      funcs_[i] << keyval->describedclassvalue("funcs", i);
    }
}

SumDenFunctional::~SumDenFunctional()
{
  if (n_) {
      for (int i=0; i < n_; i++) funcs_[i] = 0; // just in case
      delete[] funcs_;
      delete[] coefs_;
    }
  n_=0;
  funcs_=0;
  coefs_=0;
}

void
SumDenFunctional::save_data_state(StateOut& s)
{
  DenFunctional::save_data_state(s);
  s.put(n_);
  if (n_) {
      s.put(coefs_, n_);
      for (int i=0; i < n_; i++) 
          SavableState::save_state(funcs_[i].pointer(),s);
    }
}

int
SumDenFunctional::need_density_gradient()
{
  for (int i=0; i < n_; i++)
      if (funcs_[i]->need_density_gradient())
          return 1;

  return 0;
}

void
SumDenFunctional::set_spin_polarized(int p)
{
  spin_polarized_ = p;
  for (int i=0; i < n_; i++)
      funcs_[i]->set_spin_polarized(p);
}

void
SumDenFunctional::set_compute_potential(int val)
{
  compute_potential_ = val;
  for (int i=0; i < n_; i++)
      funcs_[i]->set_compute_potential(val);
}

void
SumDenFunctional::point(const PointInputData &id,
                        PointOutputData &od)
{
  od.zero();
  PointOutputData tmpod;
  for (int i=0; i < n_; i++) {
      funcs_[i]->point(id, tmpod);
      
      od.energy += coefs_[i] * tmpod.energy;
      if (compute_potential_) {
          od.df_drho_a += coefs_[i] * tmpod.df_drho_a;
          od.df_drho_b += coefs_[i] * tmpod.df_drho_b;
          od.df_dgamma_aa += coefs_[i] * tmpod.df_dgamma_aa;
          od.df_dgamma_ab += coefs_[i] * tmpod.df_dgamma_ab;
          od.df_dgamma_bb += coefs_[i] * tmpod.df_dgamma_bb;
        }
    }
}

void
SumDenFunctional::print(ostream& o) const
{
  o
    << indent << "Sum of Functionals:" << endl;
  o << incindent;
  for (int i=0; i<n_; i++) {
      o << indent << scprintf("%+18.16f",coefs_[i]) << endl;
      o << incindent;
      funcs_[i]->print(o);
      o << decindent;
    }
  o << decindent;
}

/////////////////////////////////////////////////////////////////////////////
// StdDenFunctional

static ClassDesc StdDenFunctional_cd(
  typeid(StdDenFunctional),"StdDenFunctional",1,"public SumDenFunctional",
  0, create<StdDenFunctional>, create<StdDenFunctional>);

StdDenFunctional::StdDenFunctional(StateIn& s):
  SavableState(s),
  SumDenFunctional(s),
  name_(0)
{
  s.getstring(name_);
}

StdDenFunctional::StdDenFunctional():
  name_(0)
{
}

void
StdDenFunctional::init_arrays(int n)
{
  n_ = n;
  funcs_ = new Ref<DenFunctional>[n_];
  coefs_ = new double[n_];
  for (int i=0; i<n_; i++) coefs_[i] = 1.0;
}

StdDenFunctional::StdDenFunctional(const Ref<KeyVal>& keyval)
{
  name_ = keyval->pcharvalue("name");
  if (name_) {
      if (!strcmp(name_,"HFK")) {
          n_ = 0;
          a0_ = 1.0;
        }
      else if (!strcmp(name_,"XALPHA")) {
          init_arrays(1);
          funcs_[0] = new XalphaFunctional;
        }
      else if (!strcmp(name_,"HFS")) {
          init_arrays(1);
          funcs_[0] = new SlaterXFunctional;
        }
      else if (!strcmp(name_,"HFB")) {
          init_arrays(2);
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new Becke88XFunctional;
        }
      else if (!strcmp(name_,"HFG96")) {
          init_arrays(1);
          funcs_[0] = new G96XFunctional;
        }
      else if (!strcmp(name_,"G96LYP")) {
          init_arrays(2);
          funcs_[0] = new G96XFunctional;
          funcs_[1] = new LYPCFunctional;
        }
      else if (!strcmp(name_,"BLYP")) {
          init_arrays(3);
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new Becke88XFunctional;
          funcs_[2] = new LYPCFunctional;
        }
      else if (!strcmp(name_,"SVWN1")) {
          init_arrays(2);
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new VWN1LCFunctional;
        }
      else if (!strcmp(name_,"SVWN1RPA")) {
          init_arrays(2);
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new VWN1LCFunctional(1);
        }
      else if (!strcmp(name_,"SVWN2")) {
          init_arrays(2);
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new VWN2LCFunctional;
        }
      else if (!strcmp(name_,"SVWN3")) {
          init_arrays(2);
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new VWN3LCFunctional;
        }
      else if (!strcmp(name_,"SVWN4")) {
          init_arrays(2);
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new VWN4LCFunctional;
        }
      else if (!strcmp(name_,"SVWN5")) {
          init_arrays(2);
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new VWN5LCFunctional;
        }
      else if (!strcmp(name_,"SPZ81")) {
          init_arrays(2);
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new PZ81LCFunctional;
        }
      else if (!strcmp(name_,"SPW92")) {
          init_arrays(2);
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new PW92LCFunctional;
        }
      else if (!strcmp(name_,"BPW91")) {
          init_arrays(3);
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new Becke88XFunctional;
          funcs_[2] = new PW91CFunctional;
        }
      else if (!strcmp(name_,"BP86")) {
          init_arrays(4);
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new Becke88XFunctional;
          funcs_[2] = new P86CFunctional;
          funcs_[3] = new PZ81LCFunctional;
        }
      else if (!strcmp(name_,"B3LYP")) {
          init_arrays(4);
          a0_ = 0.2;
          coefs_[0] = 0.8;
          coefs_[1] = 0.72;
          coefs_[2] = 0.19;
          coefs_[3] = 0.81;
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new Becke88XFunctional;
          funcs_[2] = new VWN1LCFunctional(1);
          funcs_[3] = new LYPCFunctional;
        }
      else if (!strcmp(name_,"B3PW91")) {
          init_arrays(4);
          a0_ = 0.2;
          coefs_[0] = 0.8;
          coefs_[1] = 0.72;
          coefs_[2] = 0.81;
          coefs_[3] = 0.19;
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new Becke88XFunctional;
          funcs_[2] = new PW91CFunctional;
          funcs_[3] = new PW92LCFunctional;
        }
      else if (!strcmp(name_,"B3P86")) {
          init_arrays(4);
          a0_ = 0.2;
          coefs_[0] = 0.8;
          coefs_[1] = 0.72;
          coefs_[2] = 0.81;
          coefs_[3] = 1.0;
          funcs_[0] = new SlaterXFunctional;
          funcs_[1] = new Becke88XFunctional;
          funcs_[2] = new P86CFunctional;
          funcs_[3] = new VWN1LCFunctional(1);
        }
      else if (!strcmp(name_,"PBE")) {
          init_arrays(2);
          funcs_[0] = new PBEXFunctional;
          funcs_[1] = new PBECFunctional;
        }
      else if (!strcmp(name_,"PW91")) {
          init_arrays(2);
          funcs_[0] = new PW91XFunctional;
          funcs_[1] = new PW91CFunctional;
        }
      else if (!strcmp(name_,"mPW(PW91)PW91")) {
          init_arrays(2);
          funcs_[0] = new mPW91XFunctional(mPW91XFunctional::PW91);
          funcs_[1] = new PW91CFunctional;
        }
      else if (!strcmp(name_,"mPWPW91")) {
          init_arrays(2);
          funcs_[0] = new mPW91XFunctional(mPW91XFunctional::mPW91);
          funcs_[1] = new PW91CFunctional;
        }
      else if (!strcmp(name_,"mPW1PW91")) {
          init_arrays(2);
          a0_ = 0.16;
          coefs_[0] = 0.84;
          coefs_[1] = 1.0;
          funcs_[0] = new mPW91XFunctional(mPW91XFunctional::mPW91);
          funcs_[1] = new PW91CFunctional;
        }
      else {
          ExEnv::out0() << "StdDenFunctional: bad name: " << name_ << endl;
          abort();
        }
    }
}

StdDenFunctional::~StdDenFunctional()
{
  delete[] name_;
}

void
StdDenFunctional::save_data_state(StateOut& s)
{
  SumDenFunctional::save_data_state(s);
  s.putstring(name_);
}

void
StdDenFunctional::print(ostream& o) const
{
  const char *n = name_;
  if (!n) n = "Null";

  o