functional.cc

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

//
// functional.cc
//
// Copyright (C) 1997 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <[e-mail]>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//

#ifdef __GNUC__
#pragma implementation
#endif

#include <math.h>

#include <util/misc/formio.h>
#include <util/state/stateio.h>
#include <chemistry/qc/dft/functional.h>

using namespace std;
using namespace sc;

#define MIN_DENSITY 1.e-14
#define MIN_GAMMA 1.e-24
#define MIN_SQRTGAMMA 1.e-12
#define MAX_ZETA 1.-1.e-12
#define MIN_ZETA -(1.-1.e-12)

///////////////////////////////////////////////////////////////////////////
// PointInputData

void
PointInputData::compute_derived(int spin_polarized, int need_gradient)
{
  a.rho_13 = pow(a.rho, 1.0/3.0);
  if (spin_polarized) {
      b.rho_13 = pow(b.rho, 1.0/3.0);
    }
  else {
      b = a;
      if (need_gradient) gamma_ab = a.gamma;
    }
}

///////////////////////////////////////////////////////////////////////////
// DenFunctional

static ClassDesc DenFunctional_cd(
  typeid(DenFunctional),"DenFunctional",1,"public SavableState",
  0, 0, 0);

DenFunctional::DenFunctional(StateIn& s):
  SavableState(s)
{
  s.get(a0_);
  s.get(spin_polarized_);
  s.get(compute_potential_);
}

DenFunctional::DenFunctional()
{
  a0_ = 0;
  spin_polarized_ = 0;
  compute_potential_ = 0;
}

DenFunctional::DenFunctional(const Ref<KeyVal>& keyval)
{
  // a0 is usually zero, except for ACM functionals. 
  a0_ = keyval->doublevalue("a0");
  spin_polarized_ = 0;
  compute_potential_ = 0;
}

DenFunctional::~DenFunctional()
{
}

void
DenFunctional::save_data_state(StateOut& s)
{
  s.put(a0_);
  s.put(spin_polarized_);
  s.put(compute_potential_);
}

int
DenFunctional::need_density_gradient()
{
  return 0;
}

int
DenFunctional::need_density_hessian()
{
  return 0;
}

void
DenFunctional::set_spin_polarized(int i)
{
  spin_polarized_ = i;
}

void
DenFunctional::set_compute_potential(int i)
{
  compute_potential_ = i;
}

void
DenFunctional::gradient(const PointInputData& id, PointOutputData& od,
                        double *grad_f, int acenter,
                        GaussianBasisSet *basis,
                        const double *dmat_a, const double *dmat_b,
                        int ncontrib, const int *contrib,
                        int ncontrib_bf_, const int *contrib_bf_,
                        const double *bs, const double *bsg,
                        const double *bsh)
{
  int need_gamma_terms = need_density_gradient();
  point(id, od);
  memset(grad_f, 0, sizeof(double)*basis->ncenter()*3);
#if 0
  ExEnv::outn() << scprintf("gradient: rho_a= %12.8f rho_b= %12.8f need_gamma = %d",
                   id.a.rho, id.b.rho, need_gamma_terms) << endl;
  ExEnv::outn() << scprintf("  gamma_aa= %12.8f gamma_bb= %12.8f gamma_ab= % 12.8f",
                   id.a.gamma, id.b.gamma, id.gamma_ab) << endl;
  ExEnv::outn() << scprintf("  df_drho_a= % 12.8f df_drho_b= % 12.8f",
                   od.df_drho_a, od.df_drho_b) << endl;
  ExEnv::outn() << scprintf("  df_dg_aa= % 12.8f df_dg_bb= % 12.8f df_dg_ab= % 12.8f",
                   od.df_dgamma_aa,od.df_dgamma_bb,od.df_dgamma_ab) << endl;
#endif

  if (need_gamma_terms) {
      double drhoa = od.df_drho_a;
      double drhob = od.df_drho_b;
      for (int nu=0; nu<ncontrib_bf_; nu++) {
          int nut = contrib_bf_[nu];
          int nuatom = basis->shell_to_center(basis->function_to_shell(nut));
          double dfa_phi_nu = drhoa * bs[nu];
          double dfb_phi_nu = drhob * bs[nu];
          for (int mu=0; mu<ncontrib_bf_; mu++) {
              int mut = contrib_bf_[mu];
              int muatom
                  = basis->shell_to_center(basis->function_to_shell(mut));
              if (muatom!=acenter) {
                  int nutmut
                      = (nut>mut?((nut*(nut+1))/2+mut):((mut*(mut+1))/2+nut));
                  double rho_a = dmat_a[nutmut];
                  double rho_b = dmat_b[nutmut];
                  int ixyz;
                  for (ixyz=0; ixyz<3; ixyz++) {
                      double contrib = -2.0*bsg[mu*3+ixyz]
                                     * (rho_a*dfa_phi_nu + rho_b*dfb_phi_nu);
#define hoff(i,j) ((j)<(i)?((i)*((i)+1))/2+(j):((j)*((j)+1))/2+(i))
                      // gamma_aa contrib
                      if (need_gamma_terms) {
                          contrib
                              += 4.0 * od.df_dgamma_aa * rho_a
                              * ( - bsg[mu*3+ixyz]
                                  * ( id.a.del_rho[0]*bsg[nu*3+0]
                                      +id.a.del_rho[1]*bsg[nu*3+1]
                                      +id.a.del_rho[2]*bsg[nu*3+2])
                                  - bs[nu]
                                  * ( id.a.del_rho[0]*bsh[mu*6+hoff(0,ixyz)]
                                      +id.a.del_rho[1]*bsh[mu*6+hoff(1,ixyz)]
                                      +id.a.del_rho[2]*bsh[mu*6+hoff(2,ixyz)]
                                      )
                                  );
                        }
                      // gamma_ab contrib
                      if (need_gamma_terms) {
                          contrib
                              += 2.0 * od.df_dgamma_ab * rho_a
                              * ( - bsg[mu*3+ixyz]
                                  * ( id.b.del_rho[0]*bsg[nu*3+0]
                                      +id.b.del_rho[1]*bsg[nu*3+1]
                                      +id.b.del_rho[2]*bsg[nu*3+2])
                                  - bs[nu]
                                  * ( id.b.del_rho[0]*bsh[mu*6+hoff(0,ixyz)]
                                      +id.b.del_rho[1]*bsh[mu*6+hoff(1,ixyz)]
                                      +id.b.del_rho[2]*bsh[mu*6+hoff(2,ixyz)]
                                      )
                                  );
                          contrib
                              += 2.0 * od.df_dgamma_ab * rho_b
                              * ( - bsg[mu*3+ixyz]
                                  * ( id.a.del_rho[0]*bsg[nu*3+0]
                                      +id.a.del_rho[1]*bsg[nu*3+1]
                                      +id.a.del_rho[2]*bsg[nu*3+2])
                                  - bs[nu]
                                  * ( id.a.del_rho[0]*bsh[mu*6+hoff(0,ixyz)]
                                      +id.a.del_rho[1]*bsh[mu*6+hoff(1,ixyz)]
                                      +id.a.del_rho[2]*bsh[mu*6+hoff(2,ixyz)]
                                      )
                                  );
                        }
                      // gamma_bb contrib
                      if (need_gamma_terms) {
                          contrib
                              += 4.0 * od.df_dgamma_bb * rho_b
                              * ( - bsg[mu*3+ixyz]
                                  * ( id.b.del_rho[0]*bsg[nu*3+0]
                                      +id.b.del_rho[1]*bsg[nu*3+1]
                                      +id.b.del_rho[2]*bsg[nu*3+2])
                                  - bs[nu]
                                  * ( id.b.del_rho[0]*bsh[mu*6+hoff(0,ixyz)]
                                      +id.b.del_rho[1]*bsh[mu*6+hoff(1,ixyz)]
                                      +id.b.del_rho[2]*bsh[mu*6+hoff(2,ixyz)]
                                      )
                                  );
                        }
                      grad_f[3*muatom+ixyz] += contrib;
                      grad_f[3*acenter+ixyz] -= contrib;
                    }
                }
            }
        }
    }
  else {
      double drhoa = od.df_drho_a;
      double drhob = od.df_drho_b;
      for (int nu=0; nu<ncontrib_bf_; nu++) {
          int nut = contrib_bf_[nu];
          int nuatom = basis->shell_to_center(basis->function_to_shell(nut));
          double dfa_phi_nu = drhoa * bs[nu];
          double dfb_phi_nu = drhob * bs[nu];
          for (int mu=0; mu<ncontrib_bf_; mu++) {
              int mut = contrib_bf_[mu];
              int muatom
                  = basis->shell_to_center(basis->function_to_shell(mut));
              if (muatom!=acenter) {
                  int nutmut
                      = (nut>mut?((nut*(nut+1))/2+mut):((mut*(mut+1))/2+nut));
                  double rho_a = dmat_a[nutmut];
                  double rho_b = dmat_b[nutmut];
                  int ixyz;
                  for (ixyz=0; ixyz<3; ixyz++) {
                      double contrib = -2.0*bsg[mu*3+ixyz]
                                     * (rho_a*dfa_phi_nu + rho_b*dfb_phi_nu);
                      grad_f[3*muatom+ixyz] += contrib;
                      grad_f[3*acenter+ixyz] -= contrib;
                    }
                }
            }
        }
    }
}

void
DenFunctional::do_fd_point(PointInputData&id,
                           double&in,double&out,
                           double lower_bound, double upper_bound)
{
  double delta = 0.0000000001;
  PointOutputData tod;
  double insave = in;

  point(id,tod);
  double outsave = tod.energy;

  int spin_polarized_save = spin_polarized_;
  set_spin_polarized(1);

  if (insave-delta>=lower_bound && insave+delta<=upper_bound) {
      in = insave+delta;
      id.compute_derived(1, need_density_gradient());
      point(id,tod);
      double plus = tod.energy;

      in = insave-delta;
      id.compute_derived(1, need_density_gradient());
      point(id,tod);
      double minu = tod.energy;
      out = 0.5*(plus-minu)/delta;
    }
  else if (insave+2*delta<=upper_bound) {
      in = insave+delta;
      id.compute_derived(1, need_density_gradient());
      point(id,tod);
      double plus = tod.energy;

      in = insave+2*delta;
      id.compute_derived(1, need_density_gradient());
      point(id,tod);
      double plus2 = tod.energy;
      out = 0.5*(4.0*plus-plus2-3.0*outsave)/delta;
    }
  else if (insave-2*delta>=lower_bound) {
      in = insave-delta;
      id.compute_derived(1, need_density_gradient());
      point(id,tod);
      double minu = tod.energy;

      in = insave-2*delta;
      id.compute_derived(1, need_density_gradient());
      point(id,tod);
      double minu2 = tod.energy;
      out = -0.5*(4.0*minu-minu2-3.0*outsave)/delta;
    }
  else {
      // the derivative is not well defined for this case
      out = -135711.;
    }
  in = insave;
  id.compute_derived(1, need_density_gradient());

  set_spin_polarized(spin_polarized_save);
}

void
DenFunctional::fd_point(const PointInputData&id, PointOutputData&od)
{
  PointInputData tid(id);

  // fill in the energy at the initial density values
  point(id,od);

  ExEnv::out0() << scprintf("ra=%7.5f rb=%7.5f gaa=%7.5f gbb=%7.5f gab= % 9.7f",
                   id.a.rho, id.b.rho, id.a.gamma, id.b.gamma, id.gamma_ab)
       << endl;

  double ga = tid.a.gamma;
  double gb = tid.b.gamma;
  double gab = tid.gamma_ab;
  double sga = sqrt(ga);
  double sgb = sqrt(gb);

  double g_a_lbound = -2*gab - gb;
  if (gb > 0 && gab*gab/gb > g_a_lbound) g_a_lbound = gab*gab/gb;
  if (g_a_lbound < 0) g_a_lbound = 0.0;

  double g_b_lbound = -2*gab - ga;
  if (ga > 0 && gab*gab/ga > g_b_lbound) g_b_lbound = gab*gab/ga;
  if (g_b_lbound < 0) g_b_lbound = 0.0;

  double g_ab_lbound = -0.5*(ga+gb);
  if (-sga*sgb > g_ab_lbound) g_ab_lbound = -sga*sgb;
  // if (-sga*sgb < g_ab_lbound) g_ab_lbound = -sga*sgb;
  double g_ab_ubound = sga*sgb;

  do_fd_point(tid, tid.a.rho, od.df_drho_a, 0.0, 10.0);
  do_fd_point(tid, tid.b.rho, od.df_drho_b, 0.0, 10.0);
  do_fd_point(tid, tid.a.gamma, od.df_dgamma_aa, g_a_lbound, 10.0);
  do_fd_point(tid, tid.b.gamma, od.df_dgamma_bb, g_b_lbound, 10.0);
  do_fd_point(tid, tid.gamma_ab, od.df_dgamma_ab, g_ab_lbound, g_ab_ubound);
}

static int
check(const char *name, double fd, double an, const char *class_name)
{
  // -135711. flags an undefined FD
  if (fd == -135711.) return 0;

  double err = fabs(fd - an);
  ExEnv::out0() << scprintf("%20s: fd = % 12.8f an = % 12.8f", name, fd, an)
       << endl;
  if ((fabs(an) > 0.03 && err/fabs(an) > 0.03)
      || ((fabs(an) <= 0.03) && err > 0.03)
      || isnan(an)) {
      ExEnv::out0() << scprintf("Error: %12s: fd = % 12.8f an = % 12.8f (%s)",
                       name, fd, an, class_name)
           << endl;
      return 1;
    }
  return 0;
}

int
DenFunctional::test(const PointInputData &id)
{
  PointOutputData fd_od;
  fd_point(id,fd_od);
  PointOutputData an_od;
  point(id,an_od);
  int r = 0;
  r+=check("df_drho_a", fd_od.df_drho_a, an_od.df_drho_a, class_name());
  r+=check("df_drho_b", fd_od.df_drho_b, an_od.df_drho_b, class_name());
  r+=check("df_dgamma_aa",fd_od.df_dgamma_aa,an_od.df_dgamma_aa,class_name());
  r+=check("df_dgamma_ab",fd_od.df_dgamma_ab,an_od.df_dgamma_ab,class_name());
  r+=check("df_dgamma_bb",fd_od.df_dgamma_bb,an_od.df_dgamma_bb,class_name());
  return r;
}

int
DenFunctional::test()
{
  int i, j, k, l, m;
  set_compute_potential(1);
  set_spin_polarized(0);
  SCVector3 r = 0.0;
  PointInputData id(r);

  for (i=0; i<6; i++) id.a.hes_rho[i] = 0.0;
  id.a.lap_rho = 0.0;

  // del rho should not be used by any of the functionals
  for (i=0; i<3; i++) id.a.del_rho[i] = id.b.del_rho[i] = 0.0;

  double testrho[] = { 0.0, 0.001, 0.5, -1 };
  double testgamma[] = { 0.0, 0.001, 0.5, -1 };
  double testgammaab[] = { -0.5, 0.0, 0.5, -1 };

  int ret = 0;

  ExEnv::out0() << "Testing with rho_a == rho_b" << endl;
  for (i=0; testrho[i] != -1.0; i++) {
      if (testrho[i] == 0.0) continue;
      id.a.rho=testrho[i];
      for (j=0; testgamma[j] != -1.0; j++) {