ossscf.cc

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

        reset_density();
      }
    }

    memcpy(ndocc_,newdocc,sizeof(int)*nirrep_);
    
    delete[] newdocc;
  }
}

void
OSSSCF::symmetry_changed()
{
  SCF::symmetry_changed();
  cl_fock_.result_noupdate()=0;
  op_focka_.result_noupdate()=0;
  op_fockb_.result_noupdate()=0;
  nirrep_ = molecule()->point_group()->char_table().ncomp();
  set_occupations(0);
}

//////////////////////////////////////////////////////////////////////////////
//
// scf things
//

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

  // 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);

  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
  cl_gmat_ = basis()->matrixkit()->symmmatrix(basis()->basisdim());
  cl_gmat_.assign(0.0);

  op_gmata_ = cl_gmat_.clone();
  op_gmata_.assign(0.0);

  op_gmatb_ = cl_gmat_.clone();
  op_gmatb_.assign(0.0);

  // test to see if we need a guess vector.
  if (cl_fock_.result_noupdate().null()) {
    cl_fock_ = hcore_.clone();
    cl_fock_.result_noupdate().assign(0.0);
    op_focka_ = hcore_.clone();
    op_focka_.result_noupdate().assign(0.0);
    op_fockb_ = hcore_.clone();
    op_fockb_.result_noupdate().assign(0.0);
  }

  // set up trial vector
  initial_vector(1);

  oso_scf_vector_ = oso_eigenvectors_.result_noupdate();
}

void
OSSSCF::done_vector()
{
  done_threads();
  
  cl_gmat_ = 0;
  cl_dens_ = 0;
  cl_dens_diff_ = 0;
  op_gmata_ = 0;
  op_densa_ = 0;
  op_densa_diff_ = 0;
  op_gmatb_ = 0;
  op_densb_ = 0;
  op_densb_diff_ = 0;

  oso_scf_vector_ = 0;
}

RefSymmSCMatrix
OSSSCF::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, 1.0);
    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
OSSSCF::alpha_density()
{
  RefSymmSCMatrix dens1(so_dimension(), basis_matrixkit());
  RefSymmSCMatrix dens2(so_dimension(), basis_matrixkit());

  so_density(dens1, 2.0);
  so_density(dens2, 1.0);
  dynamic_cast<BlockedSymmSCMatrix*>(dens2.pointer())->block(osb_)->assign(0.0);

  dens1.accumulate(dens2);
  dens2=0;

  return dens1;
}

RefSymmSCMatrix
OSSSCF::beta_density()
{
  RefSymmSCMatrix dens1(so_dimension(), basis_matrixkit());
  RefSymmSCMatrix dens2(so_dimension(), basis_matrixkit());

  so_density(dens1, 2.0);
  so_density(dens2, 1.0);
  dynamic_cast<BlockedSymmSCMatrix*>(dens2.pointer())->block(osa_)->assign(0.0);

  dens1.accumulate(dens2);
  dens2=0;

  return dens1;
}

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

  op_gmata_.assign(0.0);
  op_densa_diff_.assign(op_densa_);

  op_gmatb_.assign(0.0);
  op_densb_diff_.assign(op_densb_);
}

double
OSSSCF::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_, 1.0);

  cl_dens_.accumulate(op_densa_);

  op_densb_.assign(op_densa_);
  dynamic_cast<BlockedSymmSCMatrix*>(op_densa_.pointer())->block(osb_)->assign(0.0);
  dynamic_cast<BlockedSymmSCMatrix*>(op_densb_.pointer())->block(osa_)->assign(0.0);
  
  cl_dens_diff_.accumulate(cl_dens_);
  op_densa_diff_.accumulate(op_densa_);
  op_densb_diff_.accumulate(op_densb_);

  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
OSSSCF::scf_energy()
{
  RefSymmSCMatrix t = cl_fock_.result_noupdate().copy();
  t.accumulate(hcore_);

  RefSymmSCMatrix ga = op_focka_.result_noupdate().copy();
  ga.scale(-1.0);
  ga.accumulate(cl_fock_.result_noupdate());
  
  RefSymmSCMatrix gb = op_fockb_.result_noupdate().copy();
  gb.scale(-1.0);
  gb.accumulate(cl_fock_.result_noupdate());
  
  SCFEnergy *eop = new SCFEnergy;
  eop->reference();
  Ref<SCElementOp2> op = eop;
  t.element_op(op, cl_dens_);

  double cl_e = eop->result();
  
  eop->reset();
  ga.element_op(op, op_densa_);
  double opa_e = eop->result();

  eop->reset();
  gb.element_op(op, op_densb_);
  double opb_e = eop->result();

  op=0;
  eop->dereference();
  delete eop;

  return cl_e-opa_e-opb_e;
}

////////////////////////////////////////////////////////////////////////////
    
Ref<SCExtrapData>
OSSSCF::extrap_data()
{
  RefSymmSCMatrix *m = new RefSymmSCMatrix[3];
  m[0] = cl_fock_.result_noupdate();
  m[1] = op_focka_.result_noupdate();
  m[2] = op_fockb_.result_noupdate();
  
  Ref<SCExtrapData> data = new SymmSCMatrixNSCExtrapData(3, m);
  delete[] m;
  
  return data;
}

RefSymmSCMatrix
OSSSCF::effective_fock()
{
  // 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);

  RefSymmSCMatrix mofocka(oso_dimension(), basis_matrixkit());
  mofocka.assign(0.0);
  
  RefSymmSCMatrix mofockb(oso_dimension(), basis_matrixkit());
  mofockb.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_;
  }

  mofock.accumulate_transform(vec, fock(0),
                              SCMatrix::TransposeTransform);
  mofocka.accumulate_transform(vec, fock(1),
                               SCMatrix::TransposeTransform);
  mofockb.accumulate_transform(vec, fock(2),
                               SCMatrix::TransposeTransform);
  
  dynamic_cast<BlockedSymmSCMatrix*>(mofocka.pointer())->block(osb_)->assign(0.0);
  dynamic_cast<BlockedSymmSCMatrix*>(mofockb.pointer())->block(osa_)->assign(0.0);
  
  mofocka.accumulate(mofockb);
  mofockb=0;
  
  Ref<SCElementOp2> op = new GSGeneralEffH(this);
  mofock.element_op(op, mofocka);

  return mofock;
}

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

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

void
OSSSCF::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
OSSSCF::lagrangian()
{
  RefSCMatrix vec = so_to_orthog_so().t() * oso_scf_vector_;

  RefSymmSCMatrix mofock(oso_dimension(), basis_matrixkit());
  mofock.assign(0.0);
  mofock.accumulate_transform(vec, cl_fock_.result_noupdate(),
                              SCMatrix::TransposeTransform);

  RefSymmSCMatrix mofocka(oso_dimension(), basis_matrixkit());
  mofocka.assign(0.0);
  mofocka.accumulate_transform(vec, op_focka_.result_noupdate(),
                               SCMatrix::TransposeTransform);
  
  RefSymmSCMatrix mofockb(oso_dimension(), basis_matrixkit());
  mofockb.assign(0.0);
  mofockb.accumulate_transform(vec, op_fockb_.result_noupdate(),
                               SCMatrix::TransposeTransform);
  
  dynamic_cast<BlockedSymmSCMatrix*>(mofocka.pointer())->block(osb_)->assign(0.0);
  dynamic_cast<BlockedSymmSCMatrix*>(mofockb.pointer())->block(osa_)->assign(0.0);
  
  mofocka.accumulate(mofockb);
  mofockb=0;
  
  mofock.scale(2.0);
  
  Ref<SCElementOp2> op = new MOLagrangian(this);
  mofock.element_op(op, mofocka);
  mofocka=0;

  // transform MO lagrangian to SO basis
  RefSymmSCMatrix so_lag(so_dimension(), basis_matrixkit());
  so_lag.assign(0.0);
  so_lag.accumulate_transform(vec, 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(-1.0);

  return ao_lag;
}

RefSymmSCMatrix
OSSSCF::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_, 1.0);
  op_densb_.assign(op_densa_);

  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);
  op_densb_.scale(2.0);
  
  return tdens;
}

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

void
OSSSCF::init_hessian()
{
}

void
OSSSCF::done_hessian()
{
}

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

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