integrator.cc

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

EulerMaclaurinRadialIntegrator::print(ostream &o) const
{
  o << indent
    << scprintf("%s: nr = %d", class_name(), nr()) << endl;
}

//////////////////////////////////////////////////////////////////////////
// LebedevLaikovIntegrator

static ClassDesc LebedevLaikovIntegrator_cd(
  typeid(LebedevLaikovIntegrator),"LebedevLaikovIntegrator",1,"public AngularIntegrator",
  0, create<LebedevLaikovIntegrator>, create<LebedevLaikovIntegrator>);

LebedevLaikovIntegrator::LebedevLaikovIntegrator(StateIn& s):
  SavableState(s),
  AngularIntegrator(s)
{
  s.get(npoint_);
  init(npoint_);
}

LebedevLaikovIntegrator::LebedevLaikovIntegrator()
{
  init(302);
}

LebedevLaikovIntegrator::LebedevLaikovIntegrator(int npoint)
{
  init(npoint);
}

LebedevLaikovIntegrator::LebedevLaikovIntegrator(const Ref<KeyVal>& keyval)
{
  KeyValValueint defnpoint(302);
  
  init(keyval->intvalue("n", defnpoint));
}

LebedevLaikovIntegrator::~LebedevLaikovIntegrator()
{
  delete [] x_;
  delete [] y_;
  delete [] z_;
  delete [] w_;
}

void
LebedevLaikovIntegrator::save_data_state(StateOut& s)
{
  AngularIntegrator::save_data_state(s);
  s.put(npoint_);
}

extern "C" {
    int Lebedev_Laikov_sphere (int N, double *X, double *Y, double *Z,
                               double *W);
    int Lebedev_Laikov_npoint (int lvalue);
}

int
LebedevLaikovIntegrator::nw(void) const
{
  return npoint_;
}

void
LebedevLaikovIntegrator::init(int n)
{
  // ExEnv::outn() << " LebedevLaikovIntegrator::init -> before x_, y_, z_, and w_ malloc's " << endl;
  // ExEnv::outn() << " n = " << n << endl;
  
  x_ = new double[n];
  y_ = new double[n];
  z_ = new double[n];
  w_ = new double[n];

  // ExEnv::outn() << " LebedevLaikovIntegrator::init -> nw_points = " << n << endl;
  
  npoint_ = Lebedev_Laikov_sphere(n, x_, y_, z_, w_);
  if (npoint_ != n) {
      ExEnv::outn() << class_name() << ": bad number of points given: " << n << endl;
      abort();
    }
}

int
LebedevLaikovIntegrator::num_angular_points(double r_value, int ir)
{
  if (ir == 0) return 1;
  return npoint_;
}

double
LebedevLaikovIntegrator
::angular_point_cartesian(int iangular, double r,
                          SCVector3 &integration_point) const
{
  integration_point.x() = r*x_[iangular];
  integration_point.y() = r*y_[iangular];
  integration_point.z() = r*z_[iangular];

  return 4.0*M_PI*w_[iangular];
}

void
LebedevLaikovIntegrator::print(ostream &o) const
{
  o << indent
    << scprintf("%s:  n = %d", class_name(), npoint_) << endl;
}

/////////////////////////////////
//  GaussLegendreAngularIntegrator

static ClassDesc GaussLegendreAngularIntegrator_cd(
  typeid(GaussLegendreAngularIntegrator),"GaussLegendreAngularIntegrator",1,"public AngularIntegrator",
  0, create<GaussLegendreAngularIntegrator>, create<GaussLegendreAngularIntegrator>);

GaussLegendreAngularIntegrator::GaussLegendreAngularIntegrator(StateIn& s):
  SavableState(s),
  AngularIntegrator(s)
{
  s.get(ntheta_);
  s.get(nphi_);
  s.get(Ktheta_);
  theta_quad_weights_ = new double[ntheta_];
  theta_quad_points_ = new double[ntheta_];
}

GaussLegendreAngularIntegrator::GaussLegendreAngularIntegrator()
{
  set_ntheta(16);
  set_nphi(32);
  set_Ktheta(5);
  int ntheta = get_ntheta();
  theta_quad_weights_ = new double [ntheta];
  theta_quad_points_ = new double [ntheta];
}

GaussLegendreAngularIntegrator::GaussLegendreAngularIntegrator(const Ref<KeyVal>& keyval)
{
  set_ntheta( keyval->intvalue("ntheta") );
  if (keyval->error() != KeyVal::OK) set_ntheta(16);
  set_nphi( keyval->intvalue("nphi") );
  if (keyval->error() != KeyVal::OK) set_nphi(2*get_ntheta());
  set_Ktheta( keyval->intvalue("Ktheta") );
  if (keyval->error() != KeyVal::OK) set_Ktheta(5);

  int ntheta = get_ntheta();
  theta_quad_weights_ = new double [ntheta];
  theta_quad_points_ = new double [ntheta];
}

GaussLegendreAngularIntegrator::~GaussLegendreAngularIntegrator()
{
  delete [] theta_quad_points_;
  delete [] theta_quad_weights_;
}

void
GaussLegendreAngularIntegrator::save_data_state(StateOut& s)
{
  AngularIntegrator::save_data_state(s);
  s.put(ntheta_);
  s.put(nphi_);
  s.put(Ktheta_);
}

int
GaussLegendreAngularIntegrator::get_ntheta(void) const
{
  return ntheta_;
}

void
GaussLegendreAngularIntegrator::set_ntheta(int i)
{
  ntheta_ = i;
}

int
GaussLegendreAngularIntegrator::get_nphi(void) const
{
  return nphi_;
}

void
GaussLegendreAngularIntegrator::set_nphi(int i)
{
  nphi_ = i;
}

int
GaussLegendreAngularIntegrator::get_Ktheta(void) const
{
  return Ktheta_;
}

void
GaussLegendreAngularIntegrator::set_Ktheta(int i)
{
  Ktheta_ = i;
}

int
GaussLegendreAngularIntegrator::get_ntheta_r(void) const
{
  return ntheta_r_;
}

void
GaussLegendreAngularIntegrator::set_ntheta_r(int i)
{
  ntheta_r_ = i;
}

int
GaussLegendreAngularIntegrator::get_nphi_r(void) const
{
  return nphi_r_;
}

void
GaussLegendreAngularIntegrator::set_nphi_r(int i)
{
  nphi_r_ = i;
}

int
GaussLegendreAngularIntegrator::get_Ktheta_r(void) const
{
  return Ktheta_r_;
}

void
GaussLegendreAngularIntegrator::set_Ktheta_r(int i)
{
  Ktheta_r_ = i;
}

int
GaussLegendreAngularIntegrator::nw(void) const
{
  return nphi_*ntheta_;
}

double
GaussLegendreAngularIntegrator::sin_theta(SCVector3 &point) const
{
  return sin(point.theta());
}

int
GaussLegendreAngularIntegrator::num_angular_points(double r_value,
                                                   int ir)
{
  int Ktheta, ntheta, ntheta_r;
  
  if (ir == 0) {
      set_ntheta_r(1);
      set_nphi_r(1);
    }
  else {
      Ktheta = get_Ktheta();
      ntheta = get_ntheta();
      ntheta_r= (int) (r_value*Ktheta*ntheta);
      set_ntheta_r(ntheta_r);
      if (ntheta_r > ntheta) set_ntheta_r(ntheta);
      if (ntheta_r < 6) set_ntheta_r(6);
      set_nphi_r(2*get_ntheta_r());
    }

  gauleg(0.0, M_PI, get_ntheta_r());

  return get_ntheta_r()*get_nphi_r();
}

void
GaussLegendreAngularIntegrator::gauleg(double x1, double x2, int n)
{
  int m,j,i;
  double z1,z,xm,xl,pp,p3,p2,p1;
  const double EPS = 10.0 * DBL_EPSILON;

  m=(n+1)/2;
  xm=0.5*(x2+x1);
  xl=0.5*(x2-x1);
  for (i=1;i<=m;i++)  {
      z=cos(M_PI*(i-0.25)/(n+0.5));
      do {
          p1=1.0;
          p2=0.0;
          for (j=1;j<=n;j++) {
              p3=p2;
              p2=p1;
              p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
	    }
          pp=n*(z*p1-p2)/(z*z-1.0);
          z1=z;
          z=z1-p1/pp;
	} while (fabs(z-z1) > EPS);
      theta_quad_points_[i-1]=xm-xl*z;
      theta_quad_points_[n-i]=xm+xl*z;
      theta_quad_weights_[i-1]=2.0*xl/((1.0-z*z)*pp*pp);
      theta_quad_weights_[n-i]=theta_quad_weights_[i-1];
    }
}

double
GaussLegendreAngularIntegrator
::angular_point_cartesian(int iangular, double r,
                          SCVector3 &integration_point) const
{
  int itheta, iphi, nphi_r;

  nphi_r = get_nphi_r();
  itheta = iangular/nphi_r;
  iphi = iangular - itheta*nphi_r;
  SCVector3 point;
  point.theta() = theta_quad_points_[itheta];
  point.phi() = (double) iphi/ (double) nphi_r * 2.0 * M_PI;
  point.r() = r;
  point.spherical_to_cartesian(integration_point);
  return ( sin_theta(point)*theta_quad_weights_[itheta]*2.0*M_PI/(double)nphi_r );
}

void
GaussLegendreAngularIntegrator::print(ostream &o) const
{
  o << indent << class_name() << ":" << endl;
  o << incindent;
  o << indent << scprintf("ntheta   = %5d", get_ntheta()) << endl;
  o << indent << scprintf("nphi     = %5d", get_nphi()) << endl;
  o << indent << scprintf("Ktheta   = %5d", get_Ktheta()) << endl;
  o << decindent;
}

//////////////////////////////////////////////
//  RadialAngularIntegratorThread

class RadialAngularIntegratorThread: public DenIntegratorThread {
  protected:
    SCVector3 *centers_;
    int *nr_;
    double *atomic_radius_;
    Molecule *mol_;
    RadialAngularIntegrator *ra_integrator_;
    IntegrationWeight *weight_;
    int point_count_total_;
    double total_density_;
  public:
    RadialAngularIntegratorThread(int ithread, int nthread,
                                  RadialAngularIntegrator *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);
    ~RadialAngularIntegratorThread();
    void run();
    double total_density() { return total_density_; }
    int point_count() { return point_count_total_; }
};

RadialAngularIntegratorThread
::RadialAngularIntegratorThread(int ithread, int nthread,
                                RadialAngularIntegrator *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):
  DenIntegratorThread(ithread,nthread,
                      integrator, func,
                      alpha_dmat, beta_dmat,
                      dmat_bound,
                      linear_scaling, use_dmat_bound,
                      extent, accuracy,
                      compute_potential_integrals,
                      need_nuclear_gradient)
{
  int icenter;
  ra_integrator_ = integrator;
  int deriv_order = (need_nuclear_gradient==0?0:1);

  mol_ = integrator_->wavefunction()->molecule().pointer();

  weight_ = ra_integrator_->weight().pointer();

  nr_ = new int[natom_];
  
  for (icenter=0; icenter<natom_; icenter++)
      nr_[icenter]
	= ra_integrator_->get_radial_grid(mol_->Z(icenter),deriv_order)->nr();

  centers_ = new SCVector3[natom_];
  for (icenter=0; icenter<natom_; icenter++) {
      centers_[icenter].x() = mol_->r(icenter,0);
      centers_[icenter].y() = mol_->r(icenter,1);
      centers_[icenter].z() = mol_->r(icenter,2);
    }

  atomic_radius_ = new double[natom_];
  for (icenter=0; icenter<natom_; icenter++) {
      atomic_radius_[icenter]
          = mol_->atominfo()->maxprob_radius(mol_->Z(icenter));
    }

  point_count_total_ = 0;
  total_density_ = 0.0;
}

RadialAngularIntegratorThread::~RadialAngularIntegratorThread()
{
  delete[] centers_;
  delete[] atomic_radius_;
  delete[] nr_;
}

void
RadialAngularIntegratorThread::run()
{
  int icenter;
  int nangular;
  int ir, iangular;           // Loop indices for diff. integration dim
  int point_count;            // Counter for # integration points per center
  int nr;

  SCVector3 center;           // Cartesian position of center
  SCVector3 integration_point;

  double w,radial_multiplier,angular_multiplier;
  int deriv_order = (nuclear_gradient_==0?0:1);
        
  int parallel_counter = 0;

  for (icenter=0; icenter < natom_; icenter++) {
      point_count=0;