tcscf.cc

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

  cl_dens_ = hcore_.clone();
  cl_dens_.assign(0.0);
  
  cl_dens_diff_ = hcore_.clone();
  cl_dens_diff_.assign(0.0);

  op_densa_ = hcore_.clone();
  op_densa_.assign(0.0);
  
  op_densa_diff_ = hcore_.clone();
  op_densa_diff_.assign(0.0);

  op_densb_ = hcore_.clone();
  op_densb_.assign(0.0);
  
  op_densb_diff_ = hcore_.clone();
  op_densb_diff_.assign(0.0);

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

  ao_gmatb_ = ao_gmata_.clone();
  ao_gmatb_.assign(0.0);

  ao_ka_ = ao_gmata_.clone();
  ao_ka_.assign(0.0);

  ao_kb_ = ao_gmata_.clone();
  ao_kb_.assign(0.0);

  // test to see if we need a guess vector
  if (focka_.result_noupdate().null()) {
    focka_ = hcore_.clone();
    focka_.result_noupdate().assign(0.0);
    fockb_ = hcore_.clone();
    fockb_.result_noupdate().assign(0.0);
    ka_ = hcore_.clone();
    ka_.result_noupdate().assign(0.0);
    kb_ = hcore_.clone();
    kb_.result_noupdate().assign(0.0);
  }

  // set up trial vector
  initial_vector(1);

  oso_scf_vector_ = oso_eigenvectors_.result_noupdate();
}

void
TCSCF::done_vector()
{
  done_threads();
  
  cl_dens_ = 0;
  cl_dens_diff_ = 0;
  op_densa_ = 0;
  op_densa_diff_ = 0;
  op_densb_ = 0;
  op_densb_diff_ = 0;

  ao_gmata_ = 0;
  ao_gmatb_ = 0;
  ao_ka_ = 0;
  ao_kb_ = 0;

  oso_scf_vector_ = 0;
}

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

RefSymmSCMatrix
TCSCF::density()
{
  if (!density_.computed()) {
    RefSymmSCMatrix dens(so_dimension(), basis_matrixkit());
    RefSymmSCMatrix dens1(so_dimension(), basis_matrixkit());

    so_density(dens, 2.0);
    dens.scale(2.0);
  
    so_density(dens1, occa_);
    dens1.scale(occa_);
    dens.accumulate(dens1);
  
    so_density(dens1, occb_);
    dens1.scale(occb_);
    dens.accumulate(dens1);
    
    dens1=0;
    
    density_ = dens;
    // only flag the density as computed if the calc is converged
    if (!value_needed()) density_.computed() = 1;
  }

  return density_.result_noupdate();
}

RefSymmSCMatrix
TCSCF::alpha_density()
{
  RefSymmSCMatrix dens(so_dimension(), basis_matrixkit());
  RefSymmSCMatrix dens1(so_dimension(), basis_matrixkit());

  so_density(dens, 2.0);
  so_density(dens1, occa_);
  dens.accumulate(dens1);

  dens.scale(2.0);
  return dens;
}

RefSymmSCMatrix
TCSCF::beta_density()
{
  RefSymmSCMatrix dens(so_dimension(), basis_matrixkit());
  RefSymmSCMatrix dens1(so_dimension(), basis_matrixkit());

  so_density(dens, 2.0);
  so_density(dens1, occb_);
  dens.accumulate(dens1);

  dens.scale(2.0);
  return dens;
}

void
TCSCF::reset_density()
{
  cl_dens_diff_.assign(cl_dens_);
  
  ao_gmata_.assign(0.0);
  op_densa_diff_.assign(op_densa_);

  ao_gmatb_.assign(0.0);
  op_densb_diff_.assign(op_densb_);

  ao_ka_.assign(0.0);
  ao_kb_.assign(0.0);
}

double
TCSCF::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);

  op_densa_diff_.assign(op_densa_);
  op_densa_diff_.scale(-1.0);

  op_densb_diff_.assign(op_densb_);
  op_densb_diff_.scale(-1.0);

  so_density(cl_dens_, 2.0);
  cl_dens_.scale(2.0);

  so_density(op_densa_, occa_);
  dynamic_cast<BlockedSymmSCMatrix*>(op_densa_.pointer())->block(osb_)->assign(0.0);
  op_densa_.scale(2.0);

  so_density(op_densb_, occb_);
  dynamic_cast<BlockedSymmSCMatrix*>(op_densb_.pointer())->block(osa_)->assign(0.0);
  op_densb_.scale(2.0);

  cl_dens_diff_.accumulate(cl_dens_);
  op_densa_diff_.accumulate(op_densa_);
  op_densb_diff_.accumulate(op_densb_);

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

  return delta;
}

double
TCSCF::scf_energy()
{
  // first calculate the elements of the CI matrix
  SCFEnergy *eop = new SCFEnergy;
  eop->reference();
  Ref<SCElementOp2> op = eop;

  RefSymmSCMatrix t = focka_.result_noupdate().copy();
  t.accumulate(hcore_);

  RefSymmSCMatrix d = cl_dens_.copy();
  d.accumulate(op_densa_);

  t.element_op(op, d);
  double h11 = eop->result();

  t.assign(fockb_.result_noupdate().copy());
  t.accumulate(hcore_);

  d.assign(cl_dens_);
  d.accumulate(op_densb_);

  eop->reset();
  t.element_op(op, d);
  double h22 = eop->result();

  //t = ka_.result_noupdate();
  //eop->reset();
  //t.element_op(op, op_densb_);
  //double h21 = eop->result();

  t = kb_.result_noupdate();
  eop->reset();
  t.element_op(op, op_densa_);
  double h12 = eop->result();
  
  op=0;
  eop->dereference();
  delete eop;

  // now diagonalize the CI matrix to get the coefficients
  RefSCDimension l2 = new SCDimension(2);
  Ref<SCMatrixKit> lkit = new LocalSCMatrixKit;
  RefSymmSCMatrix h = lkit->symmmatrix(l2);
  RefSCMatrix hv = lkit->matrix(l2,l2);
  RefDiagSCMatrix hl = lkit->diagmatrix(l2);
  
  h.set_element(0,0,h11);
  h.set_element(1,1,h22);
  h.set_element(1,0,h12);
  h.diagonalize(hl,hv);

  ci1_ = hv.get_element(0,0);
  ci2_ = hv.get_element(1,0);
  double c1c2 = ci1_*ci2_;

  ExEnv::out0() << indent
               << scprintf("c1 = %10.7f c2 = %10.7f", ci1_, ci2_)
               << endl;
  
  occa_ = 2*ci1_*ci1_;
  occb_ = 2*ci2_*ci2_;
  
  double eelec = 0.5*occa_*h11 + 0.5*occb_*h22 + 2.0*c1c2*h12;
  
  return eelec;
}

Ref<SCExtrapData>
TCSCF::extrap_data()
{
  RefSymmSCMatrix *m = new RefSymmSCMatrix[4];
  m[0] = focka_.result_noupdate();
  m[1] = fockb_.result_noupdate();
  m[2] = ka_.result_noupdate();
  m[3] = kb_.result_noupdate();
  
  Ref<SCExtrapData> data = new SymmSCMatrixNSCExtrapData(4, m);
  delete[] m;
  
  return data;
}

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

  RefSymmSCMatrix mofockb(oso_dimension(), basis_matrixkit());
  mofockb.assign(0.0);

  RefSymmSCMatrix moka = mofocka.clone();
  moka.assign(0.0);

  RefSymmSCMatrix mokb = mofocka.clone();
  mokb.assign(0.0);

  // use eigenvectors if oso_scf_vector_ is null
  RefSCMatrix vec;
  if (oso_scf_vector_.null()) {
    vec = eigenvectors();
  } else {
    vec = so_to_orthog_so().t() * oso_scf_vector_;
  }
  mofocka.accumulate_transform(vec, fock(0),
                               SCMatrix::TransposeTransform);
  mofockb.accumulate_transform(vec, fock(1),
                               SCMatrix::TransposeTransform);
  moka.accumulate_transform(vec, fock(2),
                            SCMatrix::TransposeTransform);
  mokb.accumulate_transform(vec, fock(3),
                            SCMatrix::TransposeTransform);
  
  mofocka.scale(ci1_*ci1_);
  mofockb.scale(ci2_*ci2_);
  moka.scale(ci1_*ci2_);
  mokb.scale(ci1_*ci2_);

  RefSymmSCMatrix mofock = mofocka.copy();
  mofock.accumulate(mofockb);

  BlockedSymmSCMatrix *F = dynamic_cast<BlockedSymmSCMatrix*>(mofock.pointer());
  BlockedSymmSCMatrix *Fa = dynamic_cast<BlockedSymmSCMatrix*>(mofocka.pointer());
  BlockedSymmSCMatrix *Fb = dynamic_cast<BlockedSymmSCMatrix*>(mofockb.pointer());
  BlockedSymmSCMatrix *Ka = dynamic_cast<BlockedSymmSCMatrix*>(moka.pointer());
  BlockedSymmSCMatrix *Kb = dynamic_cast<BlockedSymmSCMatrix*>(mokb.pointer());
  
  double scalea = (fabs(ci1_) < fabs(ci2_)) ? 1.0/(ci1_*ci1_ + 0.05) : 1.0;
  double scaleb = (fabs(ci2_) < fabs(ci1_)) ? 1.0/(ci2_*ci2_ + 0.05) : 1.0;

  for (int b=0; b < Fa->nblocks(); b++) {
    if (b==osa_) {
      RefSymmSCMatrix f = F->block(b);
      RefSymmSCMatrix fa = Fa->block(b);
      RefSymmSCMatrix fb = Fb->block(b);
      RefSymmSCMatrix kb = Kb->block(b);

      int i,j;

      i=ndocc_[b];
      for (j=0; j < ndocc_[b]; j++) 
        f->set_element(i,j,
                       scaleb*(fb->get_element(i,j)-kb->get_element(i,j)));

      j=ndocc_[b];
      for (i=ndocc_[b]+1; i < f->n(); i++)
        f->set_element(i,j,
                       scalea*(fa->get_element(i,j)+kb->get_element(i,j)));
      
    } else if (b==osb_) {
      RefSymmSCMatrix f = F->block(b);
      RefSymmSCMatrix fa = Fa->block(b);
      RefSymmSCMatrix fb = Fb->block(b);
      RefSymmSCMatrix ka = Ka->block(b);

      int i,j;

      i=ndocc_[b];
      for (j=0; j < ndocc_[b]; j++) 
        f->set_element(i,j,
                       scalea*(fa->get_element(i,j)-ka->get_element(i,j)));

      j=ndocc_[b];
      for (i=ndocc_[b]+1; i < f->n(); i++)
        f->set_element(i,j,
                       scaleb*(fb->get_element(i,j)+ka->get_element(i,j)));
    }
  }

  return mofock;
}

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

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

void
TCSCF::done_gradient()
{
  cl_dens_=0;
  op_densa_=0;
  op_densb_=0;
  oso_scf_vector_ = 0;
}

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

// MO lagrangian
//       c    o   v
//  c  |2*FC|2*FC|0|
//     -------------
//  o  |2*FC| FO |0|
//     -------------
//  v  | 0  |  0 |0|
//
RefSymmSCMatrix
TCSCF::lagrangian()
{
  RefSCMatrix vec = so_to_orthog_so().t() * oso_scf_vector_;

  RefSymmSCMatrix mofocka = focka_.result_noupdate().clone();
  mofocka.assign(0.0);
  mofocka.accumulate_transform(vec, focka_.result_noupdate(),
                               SCMatrix::TransposeTransform);
  mofocka.scale(ci1_*ci1_);

  RefSymmSCMatrix mofockb = mofocka.clone();
  mofockb.assign(0.0);
  mofockb.accumulate_transform(vec, fockb_.result_noupdate(),
                               SCMatrix::TransposeTransform);
  mofockb.scale(ci2_*ci2_);

  // FOa = c1^2*Fa + c1c2*Kb
  RefSymmSCMatrix moka = mofocka.clone();
  moka.assign(0.0);
  moka.accumulate_transform(vec, kb_.result_noupdate(),
                            SCMatrix::TransposeTransform);
  moka.scale(ci1_*ci2_);
  moka.accumulate(mofocka);

  // FOb = c1^2*Fb + c1c2*Ka
  RefSymmSCMatrix mokb = mofocka.clone();
  mokb.assign(0.0);
  mokb.accumulate_transform(vec, ka_.result_noupdate(),
                            SCMatrix::TransposeTransform);
  mokb.scale(ci1_*ci2_);
  mokb.accumulate(mofockb);

  dynamic_cast<BlockedSymmSCMatrix*>(moka.pointer())->block(osb_)->assign(0.0);
  dynamic_cast<BlockedSymmSCMatrix*>(mokb.pointer())->block(osa_)->assign(0.0);
  
  moka.accumulate(mokb);
  mokb=0;

  // FC = c1^2*Fa + c2^2*Fb
  mofocka.accumulate(mofockb);
  mofockb=0;
  
  Ref<SCElementOp2> op = new MOLagrangian(this);
  mofocka.element_op(op, moka);
  moka=0;
  mofocka.scale(2.0);

  // transform MO lagrangian to SO basis
  RefSymmSCMatrix so_lag(so_dimension(), basis_matrixkit());
  so_lag.assign(0.0);
  so_lag.accumulate_transform(vec, mofocka);
  
  // and then from SO to AO
  Ref<PetiteList> pl = integral()->petite_list();
  RefSymmSCMatrix ao_lag = pl->to_AO_basis(so_lag);

  ao_lag.scale(-1.0);

  return ao_lag;
}

RefSymmSCMatrix
TCSCF::gradient_density()
{
  cl_dens_ = basis_matrixkit()->symmmatrix(so_dimension());
  op_densa_ = cl_dens_.clone();
  op_densb_ = cl_dens_.clone();
  
  so_density(cl_dens_, 2.0);
  cl_dens_.scale(2.0);
  
  so_density(op_densa_, occa_);
  op_densa_.scale(occa_);
  
  so_density(op_densb_, occb_);
  op_densb_.scale(occb_);
  
  dynamic_cast<BlockedSymmSCMatrix*>(op_densa_.pointer())->block(osb_)->assign(0.0);
  dynamic_cast<BlockedSymmSCMatrix*>(op_densb_.pointer())->block(osa_)->assign(0.0);
  
  Ref<PetiteList> pl = integral()->petite_list(basis());
  
  cl_dens_ = pl->to_AO_basis(cl_dens_);
  op_densa_ = pl->to_AO_basis(op_densa_);
  op_densb_ = pl->to_AO_basis(op_densb_);

  RefSymmSCMatrix tdens = cl_dens_.copy();
  tdens.accumulate(op_densa_);
  tdens.accumulate(op_densb_);

  op_densa_.scale(2.0/occa_);
  op_densb_.scale(2.0/occb_);
  
  return tdens;
}

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

void
TCSCF::init_hessian()
{
}

void
TCSCF::done_hessian()
{
}

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

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