integrator.cc

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

//
// integrator.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/timer.h>
#include <util/misc/formio.h>
#include <util/state/stateio.h>
#include <util/container/carray.h>
#include <chemistry/qc/dft/integrator.h>

using namespace std;
using namespace sc;

namespace sc {

#ifdef EXPLICIT_TEMPLATE_INSTANTIATION
template void delete_c_array2<Ref<RadialIntegrator> >(Ref<RadialIntegrator>**);
template Ref<RadialIntegrator>**
  new_c_array2<Ref<RadialIntegrator> >(int,int,Ref<RadialIntegrator>);
template void delete_c_array3<Ref<RadialIntegrator> >(Ref<RadialIntegrator>***);
template Ref<RadialIntegrator>***
  new_c_array3<Ref<RadialIntegrator> >(int,int,int,Ref<RadialIntegrator>);

template void delete_c_array2<Ref<AngularIntegrator> >(Ref<AngularIntegrator>**);
template Ref<AngularIntegrator>**
  new_c_array2<Ref<AngularIntegrator> >(int,int,Ref<AngularIntegrator>);
template void delete_c_array3<Ref<AngularIntegrator> >(Ref<AngularIntegrator>***);
template Ref<AngularIntegrator>***
  new_c_array3<Ref<AngularIntegrator> >(int,int,int,Ref<AngularIntegrator>);
#endif

//#define CHECK_ALIGN(v) if(int(&v)&7)ExEnv::outn()<<"Bad Alignment: "<< ## v <<endl;
#define CHECK_ALIGN(v)

///////////////////////////////////////////////////////////////////////////
// DenIntegratorThread

//ThreadLock *tlock;

class DenIntegratorThread: public Thread {
  protected:
    // data common to all threads
    int nthread_;
    double *alpha_dmat_;
    double *beta_dmat_;
    double *dmat_bound_;
    int nshell_;
    int nbasis_;
    int natom_;
    DenIntegrator *integrator_;
    int spin_polarized_;
    int need_hessian_;
    int need_gradient_;
    GaussianBasisSet *basis_;
    int linear_scaling_;
    int use_dmat_bound_;
    double value_;
    DenFunctional *func_;

    // not yet initialized
    ShellExtent *extent_;
    double accuracy_;
    int compute_potential_integrals_;

    // data local to thread
    int ithread_;
    int ncontrib_;
    int *contrib_;
    int ncontrib_bf_;
    int *contrib_bf_;
    double *bs_values_;
    double *bsg_values_;
    double *bsh_values_;
    double *alpha_vmat_; // lower triangle of xi_i(r) v(r) xi_j(r) integrals
    double *beta_vmat_; // lower triangle of xi_i(r) v(r) xi_j(r) integrals
    double *nuclear_gradient_;
    double *w_gradient_;
    double *f_gradient_;
    GaussianBasisSet::ValueData *valdat_;

  public:
    DenIntegratorThread(int ithread, int nthread,
                        DenIntegrator *integrator,
                        DenFunctional *func,
                        double *alpha_dmat,
                        double *beta_dmat,
                        double *dmat_bound,
                        int linear_scaling, int use_dmat_bound,
                        ShellExtent *extent,
                        double accuracy,
                        int compute_potential_integrals,
                        int need_nuclear_gradient);
    virtual ~DenIntegratorThread();
    void get_density(double *dmat, PointInputData::SpinData &d);
    double do_point(int acenter, const SCVector3 &r,
                    double weight, double multiplier,
                    double *nuclear_gradient,
                    double *f_gradient, double *w_gradient);
    double *nuclear_gradient() { return nuclear_gradient_; }
    double *alpha_vmat() { return alpha_vmat_; }
    double *beta_vmat() { return beta_vmat_; }
    double value() { return value_; }
};

DenIntegratorThread::DenIntegratorThread(int ithread, int nthread,
                                         DenIntegrator *integrator,
                                         DenFunctional *func,
                                         double *alpha_dmat,
                                         double *beta_dmat,
                                         double *dmat_bound,
                                         int linear_scaling,
                                         int use_dmat_bound,
                                         ShellExtent *extent,
                                         double accuracy,
                                         int compute_potential_integrals,
                                         int need_nuclear_gradient)
{
  value_ = 0.0;

  ithread_ = ithread;
  nthread_ = nthread;
  integrator_ = integrator;
  spin_polarized_ = integrator->wavefunction()->spin_polarized();
  nshell_ = integrator->wavefunction()->basis()->nshell();
  nbasis_ = integrator->wavefunction()->basis()->nbasis();
  natom_ = integrator->wavefunction()->molecule()->natom();
  need_gradient_ = func->need_density_gradient();
  need_hessian_ = func->need_density_hessian();
  basis_ = integrator->wavefunction()->basis().pointer();
  linear_scaling_ = linear_scaling;
  use_dmat_bound_ = use_dmat_bound;
  func_ = func;
  extent_ = extent;
  accuracy_ = accuracy;
  compute_potential_integrals_ = compute_potential_integrals;

  alpha_dmat_ = alpha_dmat;
  beta_dmat_ = beta_dmat;
  dmat_bound_ = dmat_bound;

  valdat_ = new GaussianBasisSet::ValueData(
      integrator->wavefunction()->basis(),
      integrator->wavefunction()->integral());

  if (need_nuclear_gradient) {
      nuclear_gradient_ = new double[3*natom_];
      memset(nuclear_gradient_, 0, 3*natom_*sizeof(double));
    }
  else nuclear_gradient_ = 0;

  contrib_ = new int[nshell_];

  contrib_bf_ = new int[nbasis_];

  bs_values_ = new double[nbasis_];

  alpha_vmat_ = 0;
  beta_vmat_ = 0;
  if (compute_potential_integrals_) {
      int ntri = (nbasis_*(nbasis_+1))/2;
      alpha_vmat_ = new double[ntri];
      memset(alpha_vmat_, 0, sizeof(double)*ntri);
      if (spin_polarized_) {
          beta_vmat_ = new double[ntri];
          memset(beta_vmat_, 0, sizeof(double)*ntri);
        }
    }

  bsg_values_=0;
  bsh_values_=0;
  if (need_gradient_ || nuclear_gradient_) bsg_values_ = new double[3*nbasis_];
  if (need_hessian_ || (need_gradient_ && nuclear_gradient_))
      bsh_values_ = new double[6*nbasis_];

  w_gradient_ = 0;
  f_gradient_ = 0;
  if (nuclear_gradient_) {
      w_gradient_ = new double[natom_*3];
      f_gradient_ = new double[natom_*3];
    }
}

DenIntegratorThread::~DenIntegratorThread()
{
  delete[] contrib_;
  delete[] contrib_bf_;
  delete[] bs_values_;
  delete[] bsg_values_;
  delete[] bsh_values_;
  delete[] alpha_vmat_;
  delete[] beta_vmat_;
  delete[] nuclear_gradient_;
  delete[] f_gradient_;
  delete[] w_gradient_;
  delete valdat_;
}

///////////////////////////////////////////////////////////////////////////
// DenIntegrator

static ClassDesc DenIntegrator_cd(
  typeid(DenIntegrator),"DenIntegrator",1,"public SavableState",
  0, 0, 0);

DenIntegrator::DenIntegrator(StateIn& s):
  SavableState(s)
{
  init_object();
  s.get(linear_scaling_);
  s.get(use_dmat_bound_);
}

DenIntegrator::DenIntegrator()
{
  init_object();
}

DenIntegrator::DenIntegrator(const Ref<KeyVal>& keyval)
{
  init_object();

  linear_scaling_ = keyval->booleanvalue("linear_scaling",
                                         KeyValValueboolean(linear_scaling_));
  use_dmat_bound_ = keyval->booleanvalue("use_dmat_bound",
                                         KeyValValueboolean(use_dmat_bound_));
}

DenIntegrator::~DenIntegrator()
{
  delete[] dmat_bound_;
  delete[] alpha_dmat_;
  delete[] beta_dmat_;
}

void
DenIntegrator::save_data_state(StateOut& s)
{
  s.put(linear_scaling_);
  s.put(use_dmat_bound_);
}

void
DenIntegrator::init_object()
{
  threadgrp_ = ThreadGrp::get_default_threadgrp();
  messagegrp_ = MessageGrp::get_default_messagegrp();
  compute_potential_integrals_ = 0;
  accuracy_ = DBL_EPSILON;
  linear_scaling_ = 1;
  use_dmat_bound_ = 1;
  dmat_bound_ = 0;
  alpha_dmat_ = 0;
  beta_dmat_ = 0;
  alpha_vmat_ = 0;
  beta_vmat_ = 0;
}

void
DenIntegrator::set_compute_potential_integrals(int i)
{
  compute_potential_integrals_=i;
}

void
DenIntegrator::init(const Ref<Wavefunction> &wfn)
{
  wfn_ = wfn;
  if (linear_scaling_) {
      ExEnv::out0() << indent << "Initializing ShellExtent" << endl;
      extent_ = new ShellExtent;
      extent_->init(wfn_->basis());
      ExEnv::out0() << indent
           << "  nshell = " << wfn_->basis()->nshell() << endl;
      int ncell = extent_->n(0)*extent_->n(1)*extent_->n(2);
      ExEnv::out0() << indent << "  ncell = " << ncell << endl;
      int maxval = 0;
      double ave = 0;
      for (int i=0; i<extent_->n(0); i++) {
          for (int j=0; j<extent_->n(1); j++) {
              for (int k=0; k<extent_->n(2); k++) {
                  const std::vector<ExtentData> &e
                      = extent_->contributing_shells(i,j,k);
                  int val = e.size();
                  if (val > maxval) maxval = val;
                  ave += val;
                }
            }
        }
      ave /= ncell;
      ExEnv::out0() << indent << "  ave nsh/cell = " << ave << endl;
      ExEnv::out0() << indent << "  max nsh/cell = " << maxval << endl;
    }
}

void
DenIntegrator::done()
{
  wfn_ = 0;
  extent_ = 0;
}

void
DenIntegrator::init_integration(const Ref<DenFunctional> &func,
                                const RefSymmSCMatrix& densa,
                                const RefSymmSCMatrix& densb,
                                double *nuclear_gradient)
{
  int i;
  value_ = 0.0;

  func->set_compute_potential(
      compute_potential_integrals_ || nuclear_gradient != 0);

  spin_polarized_ = wfn_->spin_polarized();
  func->set_spin_polarized(spin_polarized_);

  natom_ = wfn_->molecule()->natom();

  nshell_ = wfn_->basis()->nshell();
  nbasis_ = wfn_->basis()->nbasis();
   
  delete[] alpha_dmat_;
  RefSymmSCMatrix adens = densa;
  if (adens.null())
      adens = wfn_->alpha_ao_density();
  alpha_dmat_ = new double[(nbasis_*(nbasis_+1))/2];
  adens->convert(alpha_dmat_);

  delete[] beta_dmat_;
  beta_dmat_ = 0;
  if (spin_polarized_) {
      RefSymmSCMatrix bdens = densb;
      if (bdens.null())
          bdens = wfn_->beta_ao_density();
      beta_dmat_ = new double[(nbasis_*(nbasis_+1))/2];
      bdens->convert(beta_dmat_);
    }

  delete[] alpha_vmat_;
  delete[] beta_vmat_;
  alpha_vmat_ = 0;
  beta_vmat_ = 0;
  if (compute_potential_integrals_) {
      int ntri = (nbasis_*(nbasis_+1))/2;
      alpha_vmat_ = new double[ntri];
      memset(alpha_vmat_, 0, sizeof(double)*ntri);
      if (spin_polarized_) {
          beta_vmat_ = new double[ntri];
          memset(beta_vmat_, 0, sizeof(double)*ntri);
        }
    }

  delete[] dmat_bound_;
  dmat_bound_ = new double[(nshell_*(nshell_+1))/2];
  int ij = 0;
  for (i=0; i<nshell_; i++) {
      int ni = wfn_->basis()->shell(i).nfunction();
      for (int j=0; j<=i; j++,ij++) {
          int nj = wfn_->basis()->shell(j).nfunction();
          double bound = 0.0;
          int ibf = wfn_->basis()->shell_to_function(i);
          for (int k=0; k<ni; k++,ibf++) {
              int lmax = nj-1;
              if (i==j) lmax = k;
              int jbf = wfn_->basis()->shell_to_function(j);
              int ijbf = (ibf*(ibf+1))/2 + jbf;
              for (int l=0; l<=lmax; l++,ijbf++) {
                  double a = fabs(alpha_dmat_[ijbf]);
                  if (a > bound) bound = a;
                  if (beta_dmat_) {
                      double b = fabs(beta_dmat_[ijbf]);
                      if (b > bound) bound = b;
                    }
                }
            }
          dmat_bound_[ij] = bound;
        }
    }
}

void
DenIntegrator::done_integration()
{
  messagegrp_->sum(value_);
  if (compute_potential_integrals_) {
      int ntri = (nbasis_*(nbasis_+1))/2;
      messagegrp_->sum(alpha_vmat_,ntri);
      if (spin_polarized_) {
          messagegrp_->sum(beta_vmat_,ntri);
        }
    }
  delete[] alpha_dmat_; alpha_dmat_ = 0;
  delete[] beta_dmat_;  beta_dmat_ = 0;
  delete[] dmat_bound_; dmat_bound_ = 0;
}

inline static double
norm(double v[3])
{
  double x,y,z;
  return sqrt((x=v[0])*x + (y=v[1])*y + (z=v[2])*z);
}

inline static double
dot(double v[3], double w[3])
{
  return v[0]*w[0] + v[1]*w[1] + v[2]*w[2];
}

void
DenIntegratorThread::get_density(double *dmat, PointInputData::SpinData &d)
{
  int i, j;

  double tmp = 0.0;