integrator.cc

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

      center = centers_[icenter];
      // get current radial grid: depends on convergence threshold
      RadialIntegrator *radial
          = ra_integrator_->get_radial_grid(mol_->Z(icenter), deriv_order);
      nr = radial->nr();
      for (ir=0; ir < nr; ir++) {
          if (! (parallel_counter++%nthread_ == ithread_)) continue;
          double r = radial->radial_value(ir, nr, atomic_radius_[icenter],
                                          radial_multiplier);
          // get current angular grid: depends on radial point and threshold
          AngularIntegrator *angular
              = ra_integrator_->get_angular_grid(r, atomic_radius_[icenter],
                                                 mol_->Z(icenter),
						 deriv_order);
          nangular = angular->num_angular_points(r/atomic_radius_[icenter],ir);
          for (iangular=0; iangular<nangular; iangular++) {
              angular_multiplier
                  = angular->angular_point_cartesian(iangular,r,
                                                     integration_point);
              integration_point += center;
              w=weight_->w(icenter, integration_point, w_gradient_);
              point_count++;
              double multiplier = angular_multiplier * radial_multiplier;
              total_density_
                  += w * multiplier
                  * do_point(icenter, integration_point,
                             w, multiplier,
                             nuclear_gradient_, f_gradient_, w_gradient_);
            }
        }
      point_count_total_ += point_count;
    }
}

//////////////////////////////////////////////
//  RadialAngularIntegrator

static ClassDesc RadialAngularIntegrator_cd(
  typeid(RadialAngularIntegrator),"RadialAngularIntegrator",1,"public DenIntegrator",
  0, create<RadialAngularIntegrator>, create<RadialAngularIntegrator>);

RadialAngularIntegrator::RadialAngularIntegrator(StateIn& s):
  SavableState(s),
  DenIntegrator(s)
{
  s.get(natomic_rows_);
  s.get(max_gridtype_);
  s.get(prune_grid_);
  s.get(gridtype_);
  s.get(npruned_partitions_);
  s.get(dynamic_grids_);

//  ExEnv::outn() << "natomic_rows_ = " << natomic_rows_ << endl;
//  ExEnv::outn() << "max_gridtype_ = " << max_gridtype_ << endl;
//  ExEnv::outn() << "prune_grid_ = " << prune_grid_ << endl;
//  ExEnv::outn() << "gridtype_ = " << gridtype_ << endl;
//  ExEnv::outn() << "npruned_partitions_ = " << npruned_partitions_ << endl;
//  ExEnv::outn() << "dynamic_grids_ = " << dynamic_grids_ << endl;
  

//  ExEnv::outn() << "In StateIn Constructor!" << endl;
  weight_ = new BeckeIntegrationWeight;

  int i;
  grid_accuracy_ = new double[max_gridtype_];
  grid_accuracy_[0] = 1e-4;
  for (i=1; i<max_gridtype_; i++) grid_accuracy_[i] = grid_accuracy_[i-1]*1e-1;

  Alpha_coeffs_ = new_c_array2(natomic_rows_,npruned_partitions_-1,
                                       double(0));
  s.get_array_double(Alpha_coeffs_[0], natomic_rows_*(npruned_partitions_-1));

  radial_user_ << SavableState::restore_state(s);
  angular_user_ << SavableState::restore_state(s);

  init_default_grids();    
  set_grids();

}

RadialAngularIntegrator::RadialAngularIntegrator()
{
  weight_  = new BeckeIntegrationWeight;
  
  init_parameters();
  init_default_grids();
  set_grids();

}

RadialAngularIntegrator::RadialAngularIntegrator(const Ref<KeyVal>& keyval):
  DenIntegrator(keyval)
{
  
  radial_user_ << keyval->describedclassvalue("radial");
  angular_user_ << keyval->describedclassvalue("angular");

  weight_ << keyval->describedclassvalue("weight");
  if (weight_.null()) weight_ = new BeckeIntegrationWeight;
//  ExEnv::outn() << "In Ref<KeyVal> Constructor" << endl;
  
  init_parameters(keyval);
  init_default_grids();
  set_grids();

}

RadialAngularIntegrator::~RadialAngularIntegrator()
{
  delete_c_array2(Alpha_coeffs_);
  delete_c_array2(radial_grid_);
  delete_c_array3(angular_grid_);
  delete_c_array2(nr_points_);
  delete[] xcoarse_l_;
  delete[] grid_accuracy_;
}

void
RadialAngularIntegrator::save_data_state(StateOut& s)
{
  DenIntegrator::save_data_state(s);
  s.put(natomic_rows_);
  s.put(max_gridtype_);
  s.put(prune_grid_);
  s.put(gridtype_);
//  s.put(nr_points_[0], natomic_rows_*max_gridtype_);
//  s.put(xcoarse_l_, natomic_rows_);
  s.put(npruned_partitions_);
  s.put(dynamic_grids_);
  s.put_array_double(Alpha_coeffs_[0],natomic_rows_*(npruned_partitions_-1));
  
//  ExEnv::outn() << "natomic_rows_ = " << natomic_rows_ << endl;
//  ExEnv::outn() << "max_gridtype_ = " << max_gridtype_ << endl;
//  ExEnv::outn() << "prune_grid_ = " << prune_grid_ << endl;
//  ExEnv::outn() << "gridtype_ = " << gridtype_ << endl;
//  ExEnv::outn() << "npruned_partitions_ = " << npruned_partitions_ << endl;
//  ExEnv::outn() << "dynamic_grids_ = " << dynamic_grids_ << endl;
  
  SavableState::save_state(radial_user_.pointer(),s);
  SavableState::save_state(angular_user_.pointer(),s);
}

void
RadialAngularIntegrator::init_parameters(void)
{

  prune_grid_ = 1;
  gridtype_ = 3;
  npruned_partitions_ = 5;
  dynamic_grids_ = 1;
  max_gridtype_ = 6;
  natomic_rows_ = 5;
  grid_accuracy_ = new double[max_gridtype_];
  
  int i;
  grid_accuracy_[0] = 1e-4;
  for (i=1; i<max_gridtype_; i++) grid_accuracy_[i] = grid_accuracy_[i-1]*1e-1;
    
  init_pruning_coefficients();

  // ExEnv::outn() << "gridtype_ = " << gridtype_ << endl;
  
}

void
RadialAngularIntegrator::init_parameters(const Ref<KeyVal>& keyval)
{
  char *grid = 0;
  
  max_gridtype_ = 6;
  natomic_rows_ = 5;
  
  grid = keyval->pcharvalue("grid");

  if (grid) {
      if (!strcmp(grid,"xcoarse"))        gridtype_ = 0;
      else if (!strcmp(grid,"coarse"))    gridtype_ = 1;
      else if (!strcmp(grid,"medium"))    gridtype_ = 2;
      else if (!strcmp(grid,"fine"))      gridtype_ = 3;
      else if (!strcmp(grid,"xfine"))     gridtype_ = 4;
      else if (!strcmp(grid,"ultrafine")) gridtype_ = 5;
      else {
          ExEnv::out0()
                       << indent
                       << "ERROR: grid = \"" << grid << "\" not recognized."
                       << endl
                       << indent
                       << "Choices are: xcoarse, coarse, medium, fine, xfine."
                       << endl
                       << indent
                       << "The default is grid = fine."
                       << endl;
          abort();
        }

    }
  else {
      gridtype_ = 3;
    }


  
  //ExEnv::outn() << " gridtype = " << gridtype_ << endl;
  //ExEnv::outn() << " max_gridtype = " << max_gridtype_ << endl;
  dynamic_grids_ = keyval->intvalue("dynamic");
  if (keyval->error() != KeyVal::OK) dynamic_grids_ = 1;
  grid_accuracy_ = new double[max_gridtype_];
  //ExEnv::outn() << "init_parameters:: max_gridtype_ = " << max_gridtype_;
  
  int i;
  grid_accuracy_[0] = 1e-4;
  for (i=1; i<max_gridtype_; i++) grid_accuracy_[i] = grid_accuracy_[i-1]*1e-1;
  
  init_pruning_coefficients(keyval);

  delete[] grid;
}


void
RadialAngularIntegrator::set_grids(void)
{
  int i, j, k;

  radial_grid_ = new_c_array2(natomic_rows_,gridtype_+1,
                              Ref<RadialIntegrator>());
  angular_grid_ = new_c_array3(natomic_rows_, npruned_partitions_,
                               gridtype_+1, Ref<AngularIntegrator>());
  
  int prune_formula_1[5] = {26, 18, 12, 0, 12};  // for H to Ne
  int prune_formula_2[5] = {36, 24, 12, 0, 12};  // for Na and up

  double prune_factor[5] = {0.1, 0.4, 0.6, 1.0, 0.6};

  if (npruned_partitions_ == 1) {
      prune_factor[0] = 1.0;
      prune_formula_1[0] = 0;
      prune_formula_2[0] = 0;
    }
  else if (npruned_partitions_ != 5) {
      ExEnv::errn() << "RadialAngularIntegrator::set_grids: "
                   << "npruned_partitations must be 1 or 5" << endl;
      abort();
    }

  for (i=0; i<natomic_rows_; i++) {
      for (j=0; j<=gridtype_; j++)
          radial_grid_[i][j]
              = new EulerMaclaurinRadialIntegrator(nr_points_[i][j]);

      for (j=0; j<npruned_partitions_; j++) {
          for (k=0; k<=gridtype_; k++) {
              int grid_l = xcoarse_l_[i] + angular_grid_offset(k);
              int use_l;

              // the constant shift depends on the row and the partition
              int prune_formula;
              if (i<=1) prune_formula = prune_formula_1[j]; // H to Ne
              else prune_formula = prune_formula_2[j];      // Na and up

              // Compute the l to be used from both the constant shift and
              // the multiplicative factor method.  Use the method giving
              // the highest l.
              int use_l_formula = grid_l - prune_formula;
              int use_l_factor = int(grid_l*prune_factor[j] + 0.5);
              if (use_l_formula > use_l_factor) use_l = use_l_formula;
              else use_l = use_l_factor;

              angular_grid_[i][j][k]
                  = new LebedevLaikovIntegrator(Lebedev_Laikov_npoint(use_l));
//                ExEnv::outn() << " angular_grid_["
//                             << i << "]["
//                             << j << "]["
//                             << k
//                             << "]->nw = " << angular_grid_[i][j][k]->nw()
//                             << " xc_l = " << xcoarse_l_[i]
//                             << " off = " << angular_grid_offset(k)
//                             << " grid_l = " << grid_l
//                             << " use_l = " << use_l
//                             << endl;
            }
        }
    }
}

void
RadialAngularIntegrator::init_pruning_coefficients(void)
{
  // Set up Alpha arrays for pruning
  //ExEnv::outn() << "npruned_partitions = " << npruned_partitions_ << endl;
  //ExEnv::outn() << "natomic_rows = " << natomic_rows_ << endl;
  int num_boundaries = npruned_partitions_-1;
  Alpha_coeffs_ = new_zero_c_array2(natomic_rows_, num_boundaries,
                                            double(0));

  // set pruning cutoff variables - Alpha -> radial shell cutoffs
  init_alpha_coefficients();
   
}

void
RadialAngularIntegrator::init_pruning_coefficients(const Ref<KeyVal>& keyval)
{
  int i, j;

  prune_grid_ = keyval->booleanvalue("prune_grid");
  if (keyval->error() != KeyVal::OK) prune_grid_ = 1;

  // ExEnv::outn() << "prune_grid = " << prune_grid_ << endl;
  
  // Need to generalize this to parse input for any number of grids
  if (prune_grid_) {
      npruned_partitions_ = keyval->count("angular_points");
      if (keyval->error() != KeyVal::OK) npruned_partitions_ = 5;

      // set pruning cutoff variables - Alpha -> radial shell cutoffs
      int num_boundaries = npruned_partitions_-1;
      Alpha_coeffs_ = new_zero_c_array2(natomic_rows_, num_boundaries,
                                                double(0));
      int alpha_rows = keyval->count("alpha_coeffs");
      if (keyval->error() != KeyVal::OK) {
          if (npruned_partitions_ != 5) {
              ExEnv::outn() << " RadialAngularIntegrator:: Need to supply alpha coefficients "
                   << "for the " << num_boundaries << " partition boundaries " << endl;
              abort();
            }
          init_alpha_coefficients();
        }
      else { // alpha coefficients defined in input
          int check;
          for (i=0; i<alpha_rows; i++) {
              check = keyval->count("alpha_coeffs", i);
              if (check != num_boundaries) {
                  ExEnv::outn() << "RadialAngularIntegrator:: Number of alpha coefficients does "
                       << "not match the number of boundaries (" << check << " != "
                       << num_boundaries << ")" << endl;
                  abort();
                }
              for (j=0; j<num_boundaries; j++) 
                  Alpha_coeffs_[i][j] = keyval->doublevalue("alpha_coeffs", i, j);
            }
        }
    }
  else {
      npruned_partitions_ = 1;
      Alpha_coeffs_ = new_zero_c_array2(natomic_rows_,0,
                                                double(0));
    }
}

void
RadialAngularIntegrator::init_alpha_coefficients(void)
{
  // assumes Alpha_coeffs_ is allocated and zeroed.

  Alpha_coeffs_[0][0] = 0.25;   Alpha_coeffs_[0][1] = 0.5;
  Alpha_coeffs_[0][2] = 0.9;    Alpha_coeffs_[0][3] = 4.5;
  Alpha_coeffs_[1][0] = 0.1667; Alpha_coeffs_[1][1] = 0.5;
  Alpha_coeffs_[1][2] = 0.8;    Alpha_coeffs_[1][3] = 4.5;
  Alpha_coeffs_[2][0] = 0.1;    Alpha_coeffs_[2][1] = 0.4;
  Alpha_coeffs_[2][2] = 0.7;    Alpha_coeffs_[2][3] = 2.5;

  //  No pruning for atoms past second row
  int i;
  for (i=3; i<natomic_rows_; i++) Alpha_coeffs_[i][2] = DBL_MAX;

}

void
RadialAngularIntegrator::init_default_grids(void)
{
  xcoarse_l_ = new int[natomic_rows_];

  nr_points_ = new_c_array2(natomic_rows_,max_gridtype_,int(0));

  // Set angular momentum level of reference xcoarse grids for each atomic row
  xcoarse_l_[0] = 17; xcoarse_l_[1] = 17; xcoarse_l_[2] = 21;
  xcoarse_l_[3] = 25; xcoarse_l_[4] = 31;

  // Set number of radial points for each atomic row and grid type
  nr_points_[0][0] = 30; nr_points_[0][1] = 50; nr_points_[0][2] = 75;
  nr_points_[0][3] = 85; nr_points_[0][4] = 100; nr_points_[0][5] = 140;

  nr_points_[1][0] = 30; nr_points_[1][1] = 50; nr_points_[1][2] = 75;
  nr_points_[1][3] = 85; nr_points_[1][4] = 100; nr_points_[1][5] = 140;

  nr_points_[2][0] = 45; nr_points_[2][1] = 75; nr_points_[2][2] = 95;
  nr_points_[2][3] = 110; nr_points_[2][4] = 125; nr_points_[2][5] = 175;

  nr_points_[3][0] = 75; nr_points_[3][1] = 95; nr_points_[3][2] = 110;
  nr_points_[3][3] = 130; nr_points_[3][4] = 160; nr_points_[3][5] = 210;

  nr_points_[4][0] = 105; nr_points_[4][1] = 130; nr_points_[4][2] = 155;
  nr_points_[4][3] = 180; nr_points_[4][4] = 205; nr_points_[4][5] = 235;

  // prune_grid_ = 1; npruned_partitions_ = 5; gridtype_ = 2;
}

int
RadialAngularIntegrator::angular_grid_offset(int gridtype)
{
  switch (gridtype) {
  case 0: return 0;
  case 1: return 6;
  case 2: return 12;
  case 3: return 18;
  case 4: return 30;
  case 5: return 42;
  default: abort();
    }
  return 0;
}

RadialIntegrator *
RadialAngularIntegrator::get_radial_grid(int charge, int deriv_order)
{
  if (radial_user_.null()) {

      int select_grid;
      
      if (dynamic_grids_ && deriv_order == 0)
	  select_grid = select_dynamic_grid();
      else select_grid = gridtype_;
      //ExEnv::out << "RAI::get_radial_grid -> select_grid = " << select_grid;
      
      if (charge<3) return radial_grid_[0][select_grid].pointer();
      else if (charge<11) return radial_grid_[1][select_grid].pointer();
      else if (charge<19) return radial_grid_[2][select_grid].pointer();
      else if (charge<37) return radial_grid_[3][select_grid].pointer();
      else if (charge<55) return radial_grid_[4][select_grid].pointer();
      else {
          ExEnv::outn() << " No default radial grids for atomic charge " << charge << endl;
          abort();
        }
    }

  return radial_user_.pointer();

}

int
RadialAngularIntegrato