clscf.cc

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

void
CLSCF::init_vector()
{
  init_threads();

  // initialize the two electron integral classes
  ExEnv::out0() << indent
       << "integral intermediate storage = " << integral()->storage_used()
       << " bytes" << endl;
  ExEnv::out0() << indent
       << "integral cache = " << integral()->storage_unused()
       << " bytes" << endl;

  // allocate storage for other temp matrices
  cl_dens_ = hcore_.clone();
  cl_dens_.assign(0.0);
  
  cl_dens_diff_ = hcore_.clone();
  cl_dens_diff_.assign(0.0);

  // gmat is in AO basis
  cl_gmat_ = basis()->matrixkit()->symmmatrix(basis()->basisdim());
  cl_gmat_.assign(0.0);

  if (cl_fock_.result_noupdate().null()) {
    cl_fock_ = hcore_.clone();
    cl_fock_.result_noupdate().assign(0.0);
  }

  // set up trial vector
  initial_vector(1);

  oso_scf_vector_ = oso_eigenvectors_.result_noupdate();

  if (accumddh_.nonnull()) accumddh_->init(this);
}

void
CLSCF::done_vector()
{
  done_threads();

  if (accumddh_.nonnull()) {
      accumddh_->print_summary();
      accumddh_->done();
  }

  cl_gmat_ = 0;
  cl_dens_ = 0;
  cl_dens_diff_ = 0;

  oso_scf_vector_ = 0;
}

void
CLSCF::reset_density()
{
  cl_gmat_.assign(0.0);
  cl_dens_diff_.assign(cl_dens_);
}

RefSymmSCMatrix
CLSCF::density()
{
  if (!density_.computed()) {
    RefSymmSCMatrix dens(so_dimension(), basis_matrixkit());
    so_density(dens, 2.0);
    dens.scale(2.0);

    density_ = dens;
    // only flag the density as computed if the calc is converged
    if (!value_needed()) density_.computed() = 1;
  }

  return density_.result_noupdate();
}

double
CLSCF::new_density()
{
  // copy current density into density diff and scale by -1.  later we'll
  // add the new density to this to get the density difference.
  cl_dens_diff_.assign(cl_dens_);
  cl_dens_diff_.scale(-1.0);
  
  so_density(cl_dens_, 2.0);
  cl_dens_.scale(2.0);

  cl_dens_diff_.accumulate(cl_dens_);
  
  Ref<SCElementScalarProduct> sp(new SCElementScalarProduct);
  cl_dens_diff_.element_op(sp.pointer(), cl_dens_diff_);
  
  double delta = sp->result();
  delta = sqrt(delta/i_offset(cl_dens_diff_.n()));

  return delta;
}

double
CLSCF::scf_energy()
{
  RefSymmSCMatrix t = cl_fock_.result_noupdate().copy();
  t.accumulate(hcore_);

  SCFEnergy *eop = new SCFEnergy;
  eop->reference();
  Ref<SCElementOp2> op = eop;
  t.element_op(op,cl_dens_);
  op=0;
  eop->dereference();

  double eelec = eop->result();

  delete eop;

  return eelec;
}

Ref<SCExtrapData>
CLSCF::extrap_data()
{
  Ref<SCExtrapData> data =
    new SymmSCMatrixSCExtrapData(cl_fock_.result_noupdate());
  return data;
}

RefSymmSCMatrix
CLSCF::effective_fock()
{
  // just return MO fock matrix.  use fock() instead of cl_fock_ just in
  // case this is called from someplace outside SCF::compute_vector()
  RefSymmSCMatrix mofock(oso_dimension(), basis_matrixkit());
  mofock.assign(0.0);

  if (debug_ > 1) {
    fock(0).print("CL Fock matrix in SO basis");
  }

  // use eigenvectors if scf_vector_ is null
  if (oso_scf_vector_.null())
    mofock.accumulate_transform(eigenvectors(), fock(0),
                                SCMatrix::TransposeTransform);
  else
    mofock.accumulate_transform(so_to_orthog_so().t() * oso_scf_vector_,
                                fock(0),
                                SCMatrix::TransposeTransform);

  return mofock;
}

//////////////////////////////////////////////////////////////////////////////

void
CLSCF::init_gradient()
{
  // presumably the eigenvectors have already been computed by the time
  // we get here
  oso_scf_vector_ = oso_eigenvectors_.result_noupdate();
}

void
CLSCF::done_gradient()
{
  cl_dens_=0;
  oso_scf_vector_ = 0;
}

/////////////////////////////////////////////////////////////////////////////

class CLLag : public BlockedSCElementOp {
  private:
    CLSCF *scf_;

  public:
    CLLag(CLSCF* s) : scf_(s) {}
    ~CLLag() {}

    int has_side_effects() { return 1; }

    void process(SCMatrixBlockIter& bi) {
      int ir=current_block();

      for (bi.reset(); bi; bi++) {
        double occi = scf_->occupation(ir,bi.i());

        if (occi==0.0)
          bi.set(0.0);
      }
    }
};

RefSymmSCMatrix
CLSCF::lagrangian()
{
  // the MO lagrangian is just the eigenvalues of the occupied MO's
  RefSymmSCMatrix mofock = effective_fock();

  Ref<SCElementOp> op = new CLLag(this);
  mofock.element_op(op);
  
  // transform MO lagrangian to SO basis
  RefSymmSCMatrix so_lag(so_dimension(), basis_matrixkit());
  so_lag.assign(0.0);
  so_lag.accumulate_transform(so_to_orthog_so().t() * oso_scf_vector_, mofock);
  
  // and then from SO to AO
  Ref<PetiteList> pl = integral()->petite_list();
  RefSymmSCMatrix ao_lag = pl->to_AO_basis(so_lag);
  ao_lag->scale(-2.0);

  return ao_lag;
}

RefSymmSCMatrix
CLSCF::gradient_density()
{
  cl_dens_ = basis_matrixkit()->symmmatrix(so_dimension());
  cl_dens_.assign(0.0);
  
  so_density(cl_dens_, 2.0);
  cl_dens_.scale(2.0);
  
  Ref<PetiteList> pl = integral()->petite_list(basis());
  
  cl_dens_ = pl->to_AO_basis(cl_dens_);

  return cl_dens_;
}

/////////////////////////////////////////////////////////////////////////////

void
CLSCF::init_hessian()
{
}

void
CLSCF::done_hessian()
{
}

/////////////////////////////////////////////////////////////////////////////

void
CLSCF::two_body_deriv_hf(double * tbgrad, double exchange_fraction)
{
  int i;
  int na3 = molecule()->natom()*3;
  int nthread = threadgrp_->nthread();

  tim_enter("setup");
  Ref<SCElementMaxAbs> m = new SCElementMaxAbs;
  cl_dens_.element_op(m.pointer());
  double pmax = m->result();
  m=0;

  double **grads = new double*[nthread];
  Ref<TwoBodyDerivInt> *tbis = new Ref<TwoBodyDerivInt>[nthread];
  for (i=0; i < nthread; i++) { 
    tbis[i] = integral()->electron_repulsion_deriv();
    grads[i] = new double[na3];
    memset(grads[i], 0, sizeof(double)*na3);
  }
  
  Ref<PetiteList> pl = integral()->petite_list();
  
  tim_change("contribution");
  
  // now try to figure out the matrix specialization we're dealing with.
  // if we're using Local matrices, then there's just one subblock, or
  // see if we can convert P to a local matrix
  if (local_ || local_dens_) {
    double *pmat;
    RefSymmSCMatrix ptmp = get_local_data(cl_dens_, pmat, SCF::Read);

    LocalCLHFGradContribution l(pmat);
    LocalTBGrad<LocalCLHFGradContribution> **tblds =
      new LocalTBGrad<LocalCLHFGradContribution>*[nthread];

    for (i=0; i < nthread; i++) {
      tblds[i] = new LocalTBGrad<LocalCLHFGradContribution>(
        l, tbis[i], pl, basis(), scf_grp_, grads[i], pmax,
        desired_gradient_accuracy(), nthread, i, exchange_fraction);
      threadgrp_->add_thread(i, tblds[i]);
    }

    tim_enter("start thread");
    if (threadgrp_->start_threads() < 0) {
      ExEnv::err0() << indent
           << "CLSCF: error starting threads" << endl;
      abort();
    }
    tim_exit("start thread");

    tim_enter("stop thread");
    if (threadgrp_->wait_threads() < 0) {
      ExEnv::err0() << indent
           << "CLSCF: error waiting for threads" << endl;
      abort();
    }
    tim_exit("stop thread");

    for (i=0; i < nthread; i++) {
      if (i) {
        for (int j=0; j < na3; j++)
          grads[0][j] += grads[i][j];

        delete[] grads[i];
      }

      delete tblds[i];
    }

    if (scf_grp_->n() > 1)
      scf_grp_->sum(grads[0], na3);

    for (i=0; i < na3; i++)
      tbgrad[i] += grads[0][i];
    
    delete[] grads[0];
    delete[] tblds;
    delete[] grads;
  }

  // for now quit
  else {
    ExEnv::err0() << indent
         << "CLHF::two_body_deriv: can't do gradient yet\n";
    abort();
  }

  for (i=0; i < nthread; i++)
    tbis[i] = 0;
  delete[] tbis;
  
  tim_exit("contribution");
}

/////////////////////////////////////////////////////////////////////////////

}

// Local Variables:
// mode: c++
// c-file-style: "ETS"
// End: