uscf.cc

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

            oso_eigenvectors_beta_ = guess_wfn_->oso_beta_eigenvectors().copy();
            eigenvalues_beta_ = guess_wfn_->beta_eigenvalues().copy();
          } else if (ug) {
            oso_eigenvectors_ = ug->oso_eigenvectors_.result_noupdate().copy();
            eigenvalues_ = ug->eigenvalues_.result_noupdate().copy();
            oso_eigenvectors_beta_
              = ug->oso_eigenvectors_beta_.result_noupdate().copy();
            eigenvalues_beta_ = ug->eigenvalues_beta_.result_noupdate().copy();
          }
          ExEnv::out0() << decindent << decindent;
        } else {
          ExEnv::out0() << indent
               << "Projecting guess wavefunction into the present basis set"
               << endl;

          // indent output of projected_eigenvectors() call if there is any
          ExEnv::out0() << incindent << incindent;
          oso_eigenvectors_ = projected_eigenvectors(guess_wfn_, 1);
          eigenvalues_ = projected_eigenvalues(guess_wfn_, 1);
          oso_eigenvectors_beta_ = projected_eigenvectors(guess_wfn_, 0);
          eigenvalues_beta_ = projected_eigenvalues(guess_wfn_, 0);
          ExEnv::out0() << decindent << decindent;
        }

        // we should only have to do this once, so free up memory used
        // for the old wavefunction, unless told otherwise
        if (!keep_guess_wfn_) guess_wfn_=0;

        ExEnv::out0() << endl;
      
      } else {
        ExEnv::out0() << indent << "Starting from core Hamiltonian guess\n"
             << endl;
        oso_eigenvectors_ = hcore_guess(eigenvalues_.result_noupdate());
        oso_eigenvectors_beta_ = oso_eigenvectors_.result_noupdate().copy();
        eigenvalues_beta_ = eigenvalues_.result_noupdate().copy();
      }
    } else {
      // this is just an old vector
    }
  }

  need_vec_=needv;
}

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

void
UnrestrictedSCF::set_occupations(const RefDiagSCMatrix& ev)
{
  abort();
}

void
UnrestrictedSCF::set_occupations(const RefDiagSCMatrix& eva,
                                 const RefDiagSCMatrix& evb)
{
  if (user_occupations_ || (initial_nalpha_ && eva.null())) {
    if (form_occupations(nalpha_, initial_nalpha_)) {
      form_occupations(nbeta_, initial_nbeta_);
      most_recent_pg_ = new PointGroup(molecule()->point_group());
      return;
    }
    ExEnv::out0() << indent
         << "UnrestrictedSCF: WARNING: reforming occupation vector from scratch" << endl;
  }
  
  if (nirrep_==1) {
    delete[] nalpha_;
    nalpha_=new int[1];
    nalpha_[0] = tnalpha_;
    delete[] nbeta_;
    nbeta_=new int[1];
    nbeta_[0] = tnbeta_;
    if (!initial_nalpha_ && initial_pg_->equiv(molecule()->point_group())) {
      initial_nalpha_=new int[1];
      initial_nalpha_[0] = tnalpha_;
    }
    if (!initial_nbeta_ && initial_pg_->equiv(molecule()->point_group())) {
      initial_nbeta_=new int[1];
      initial_nbeta_[0] = tnbeta_;
    }
    return;
  }
  
  int i,j;
  
  RefDiagSCMatrix evalsa, evalsb;
  
  if (eva.null()) {
    initial_vector(0);
    evalsa = eigenvalues_.result_noupdate();
    evalsb = eigenvalues_beta_.result_noupdate();
  }
  else {
    evalsa = eva;
    evalsb = evb;
  }

  // first convert evals to something we can deal with easily
  BlockedDiagSCMatrix *bevalsa = require_dynamic_cast<BlockedDiagSCMatrix*>(evalsa,
                                "UnrestrictedSCF::set_occupations");
  BlockedDiagSCMatrix *bevalsb = require_dynamic_cast<BlockedDiagSCMatrix*>(evalsb,
                                "UnrestrictedSCF::set_occupations");
  
  double **valsa = new double*[nirrep_];
  double **valsb = new double*[nirrep_];
  for (i=0; i < nirrep_; i++) {
    int nf=oso_dimension()->blocks()->size(i);
    if (nf) {
      valsa[i] = new double[nf];
      valsb[i] = new double[nf];
      bevalsa->block(i)->convert(valsa[i]);
      bevalsb->block(i)->convert(valsb[i]);
    } else {
      valsa[i] = 0;
      valsb[i] = 0;
    }
  }

  // now loop to find the tnalpha_ lowest eigenvalues and populate those
  // MO's
  int *newalpha = new int[nirrep_];
  memset(newalpha,0,sizeof(int)*nirrep_);

  for (i=0; i < tnalpha_; i++) {
    // find lowest eigenvalue
    int lir=0,ln=0;
    double lowest=999999999;

    for (int ir=0; ir < nirrep_; ir++) {
      int nf=oso_dimension()->blocks()->size(ir);
      if (!nf)
        continue;
      for (j=0; j < nf; j++) {
        if (valsa[ir][j] < lowest) {
          lowest=valsa[ir][j];
          lir=ir;
          ln=j;
        }
      }
    }
    valsa[lir][ln]=999999999;
    newalpha[lir]++;
  }

  int *newbeta = new int[nirrep_];
  memset(newbeta,0,sizeof(int)*nirrep_);

  for (i=0; i < tnbeta_; i++) {
    // find lowest eigenvalue
    int lir=0,ln=0;
    double lowest=999999999;

    for (int ir=0; ir < nirrep_; ir++) {
      int nf=oso_dimension()->blocks()->size(ir);
      if (!nf)
        continue;
      for (j=0; j < nf; j++) {
        if (valsb[ir][j] < lowest) {
          lowest=valsb[ir][j];
          lir=ir;
          ln=j;
        }
      }
    }
    valsb[lir][ln]=999999999;
    newbeta[lir]++;
  }

  // get rid of vals
  for (i=0; i < nirrep_; i++) {
    if (valsa[i])
      delete[] valsa[i];
    if (valsb[i])
      delete[] valsb[i];
  }
  delete[] valsa;
  delete[] valsb;

  if (!nalpha_) {
    nalpha_=newalpha;
    nbeta_=newbeta;
  } else if (most_recent_pg_.nonnull()
             && most_recent_pg_->equiv(molecule()->point_group())) {
    // test to see if newocc is different from nalpha_
    for (i=0; i < nirrep_; i++) {
      if (nalpha_[i] != newalpha[i]) {
        ExEnv::err0() << indent << "UnrestrictedSCF::set_occupations:  WARNING!!!!\n"
             << incindent << indent
             << scprintf("occupations for irrep %d have changed\n",i+1)
             << indent
             << scprintf("nalpha was %d, changed to %d", nalpha_[i], newalpha[i])
             << endl << decindent;
      }
      if (nbeta_[i] != newbeta[i]) {
        ExEnv::err0() << indent << "UnrestrictedSCF::set_occupations:  WARNING!!!!\n"
             << incindent << indent
             << scprintf("occupations for irrep %d have changed\n",i+1)
             << indent
             << scprintf("nbeta was %d, changed to %d", nbeta_[i], newbeta[i])
             << endl << decindent;
      }
    }

    memcpy(nalpha_,newalpha,sizeof(int)*nirrep_);
    memcpy(nbeta_,newbeta,sizeof(int)*nirrep_);
    delete[] newalpha;
    delete[] newbeta;
  }

  if (initial_pg_->equiv(molecule()->point_group())) {
    delete[] initial_nalpha_;
    initial_nalpha_ = new int[nirrep_];
    memcpy(initial_nalpha_,nalpha_,sizeof(int)*nirrep_);
  }

  if (initial_pg_->equiv(molecule()->point_group())) {
    delete[] initial_nbeta_;
    initial_nbeta_ = new int[nirrep_];
    memcpy(initial_nbeta_,nbeta_,sizeof(int)*nirrep_);
  }

  most_recent_pg_ = new PointGroup(molecule()->point_group());
}

void
UnrestrictedSCF::symmetry_changed()
{
  SCF::symmetry_changed();
  nirrep_ = molecule()->point_group()->char_table().ncomp();
  oso_eigenvectors_beta_.result_noupdate() = 0;
  eigenvalues_beta_.result_noupdate() = 0;
  focka_.result_noupdate() = 0;
  fockb_.result_noupdate() = 0;
  set_occupations(0,0);
}

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

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

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

  densb_ = hcore_.clone();
  densb_.assign(0.0);
  
  diff_densb_ = hcore_.clone();
  diff_densb_.assign(0.0);

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

  gmatb_ = gmata_.clone();
  gmatb_.assign(0.0);

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

  // set up trial vector
  initial_vector(1);
    
  oso_scf_vector_ = oso_eigenvectors_.result_noupdate();
  oso_scf_vector_beta_ = oso_eigenvectors_beta_.result_noupdate();
}

void
UnrestrictedSCF::done_vector()
{
  done_threads();
  
  hcore_ = 0;
  gmata_ = 0;
  densa_ = 0;
  diff_densa_ = 0;
  gmatb_ = 0;
  densb_ = 0;
  diff_densb_ = 0;

  oso_scf_vector_ = 0;
  oso_scf_vector_beta_ = 0;
}

RefSymmSCMatrix
UnrestrictedSCF::alpha_density()
{
  RefSymmSCMatrix dens(so_dimension(), basis_matrixkit());
  so_density(dens, 1.0, 1);
  return dens;
}

RefSymmSCMatrix
UnrestrictedSCF::beta_density()
{
  RefSymmSCMatrix dens(so_dimension(), basis_matrixkit());
  so_density(dens, 1.0, 0);
  return dens;
}

void
UnrestrictedSCF::reset_density()
{
  gmata_.assign(0.0);
  diff_densa_.assign(densa_);

  gmatb_.assign(0.0);
  diff_densb_.assign(densb_);
}

double
UnrestrictedSCF::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.
  diff_densa_.assign(densa_);
  diff_densa_.scale(-1.0);

  diff_densb_.assign(densb_);
  diff_densb_.scale(-1.0);

  so_density(densa_, 1.0, 1);
  so_density(densb_, 1.0, 0);

  diff_densa_.accumulate(densa_);
  diff_densb_.accumulate(densb_);

  RefSymmSCMatrix d = diff_densa_ + diff_densb_;

  Ref<SCElementScalarProduct> sp(new SCElementScalarProduct);
  d.element_op(sp.pointer(), d);
  d=0;
  
  double delta = sp->result();
  delta = sqrt(delta/i_offset(diff_densa_.n()));

  return delta;
}

RefSymmSCMatrix
UnrestrictedSCF::density()
{
  if (!density_.computed()) {
    RefSymmSCMatrix densa(so_dimension(), basis_matrixkit());
    RefSymmSCMatrix densb(so_dimension(), basis_matrixkit());
    so_density(densa, 1.0, 1);
    so_density(densb, 1.0, 0);
    densa.accumulate(densb);
    densb=0;
    
    density_ = densa;
    // only flag the density as computed if the calc is converged
    if (!value_needed()) density_.computed() = 1;
  }

  return density_.result_noupdate();
}

double
UnrestrictedSCF::scf_energy()
{
  SCFEnergy *eop = new SCFEnergy;
  eop->reference();
  Ref<SCElementOp2> op = eop;
  focka_.result_noupdate().element_op(op, densa_);
  double ea = eop->result();
  
  eop->reset();
  fockb_.result_noupdate().element_op(op, densb_);
  double eb = eop->result();

  RefSymmSCMatrix denst = densa_+densb_;
  eop->reset();
  hcore_.element_op(op, denst);
  double ec = eop->result();
  denst=0;
  
  op=0;
  eop->dereference();
  delete eop;

  return ec+ea+eb;
}

RefSymmSCMatrix
UnrestrictedSCF::effective_fock()
{
  abort();
  return 0;
}

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

class UAExtrapErrorOp : public BlockedSCElementOp {
  private:
    UnrestrictedSCF *scf_;

  public:
    UAExtrapErrorOp(UnrestrictedSCF *s) : scf_(s) {}
    ~UAExtrapErrorOp() {}

    int has_side_effects() { return 1; }

    void process(SCMatrixBlockIter& bi) {
      int ir=current_block();
      
      for (bi.reset(); bi; bi++) {
        int i=bi.i();
        int j=bi.j();
        if (scf_->alpha_occupation(ir,i) == scf_->alpha_occupation(ir,j))
          bi.set(0.0);
      }
    }
};

class UBExtrapErrorOp : public BlockedSCElementOp {
  private:
    UnrestrictedSCF *scf_;

  public: