scf.cc

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

      l->assign(0.0);
  } else if (scf_grp_->n() > 1 && access==Accum) {
    l = m.clone();
    l.assign(0.0);
  }

  if (dynamic_cast<ReplSymmSCMatrix*>(l.pointer()))
    p = dynamic_cast<ReplSymmSCMatrix*>(l.pointer())->get_data();
  else
    p = dynamic_cast<LocalSymmSCMatrix*>(l.pointer())->get_data();

  return l;
}

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

void
SCF::initial_vector(int needv)
{
  if (need_vec_) {
    if (always_use_guess_wfn_ || oso_eigenvectors_.result_noupdate().null()) {
      // if guess_wfn_ is non-null then try to get a guess vector from it.
      // First check that the same basis is used...if not, then project the
      // guess vector into the present basis.
      if (guess_wfn_.nonnull()) {
        if (basis()->equiv(guess_wfn_->basis())
            &&orthog_method() == guess_wfn_->orthog_method()
            &&oso_dimension()->equiv(guess_wfn_->oso_dimension().pointer())) {
          ExEnv::out0() << indent
               << "Using guess wavefunction as starting vector" << endl;

          // indent output of eigenvectors() call if there is any
          ExEnv::out0() << incindent << incindent;
          SCF *g = dynamic_cast<SCF*>(guess_wfn_.pointer());
          if (!g || compute_guess_) {
            oso_eigenvectors_ = guess_wfn_->oso_eigenvectors().copy();
            eigenvalues_ = guess_wfn_->eigenvalues().copy();
          } else {
            oso_eigenvectors_ = g->oso_eigenvectors_.result_noupdate().copy();
            eigenvalues_ = g->eigenvalues_.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_);
          eigenvalues_ = projected_eigenvalues(guess_wfn_);
          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());
      }
    } else {
      // this is just an old vector
    }
  }

  need_vec_=needv;
}

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

void
SCF::init_mem(int nm)
{
  // if local_den_ is already 0, then that means it was set to zero by
  // the user.
  if (!local_dens_) {
    integral()->set_storage(storage_);
    return;
  }
  
  size_t nmem = i_offset(basis()->nbasis())*nm*sizeof(double);

  // if we're actually using local matrices, then there's no choice
  if (dynamic_cast<LocalSCMatrixKit*>(basis()->matrixkit().pointer())
      ||dynamic_cast<ReplSCMatrixKit*>(basis()->matrixkit().pointer())) {
    if (nmem > storage_)
      return;
  } else {
    if (nmem > storage_) {
      local_dens_=0;
      integral()->set_storage(storage_);
      return;
    }
  }

  integral()->set_storage(storage_-nmem);
}

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

void
SCF::so_density(const RefSymmSCMatrix& d, double occ, int alp)
{
  int i,j,k;
  int me=scf_grp_->me();
  int nproc=scf_grp_->n();
  int uhf = spin_unrestricted();
  
  d->assign(0.0);
  
  RefSCMatrix oso_vector;
  if (alp || !uhf) {
    if (oso_scf_vector_.nonnull())
      oso_vector = oso_scf_vector_;
  }
  else {
    if (oso_scf_vector_beta_.nonnull())
      oso_vector = oso_scf_vector_beta_;
  }
      
  if (oso_vector.null()) {
    if (uhf) {
      if (alp)
        oso_vector = oso_alpha_eigenvectors();
      else
        oso_vector = oso_beta_eigenvectors();
    } else
      oso_vector = oso_eigenvectors();
  }

  if (debug_ > 1) oso_vector.print("ortho SO vector");

  RefSCMatrix vector = so_to_orthog_so().t() * oso_vector;
  oso_vector = 0;

  if (debug_ > 1) vector.print("SO vector");
  
  BlockedSCMatrix *bvec = require_dynamic_cast<BlockedSCMatrix*>(
    vector, "SCF::so_density: blocked vector");

  BlockedSymmSCMatrix *bd = require_dynamic_cast<BlockedSymmSCMatrix*>(
    d, "SCF::so_density: blocked density");
  
  for (int ir=0; ir < oso_dimension()->blocks()->nblock(); ir++) {
    RefSCMatrix vir = bvec->block(ir);
    RefSymmSCMatrix dir = bd->block(ir);
    
    if (vir.null() || vir.ncol()==0)
      continue;
    
    int n_orthoSO = oso_dimension()->blocks()->size(ir);
    int n_SO = so_dimension()->blocks()->size(ir);
    
    // figure out which columns of the scf vector we'll need
    int col0 = -1, coln = -1;
    for (i=0; i < n_orthoSO; i++) {
      double occi;
      if (!uhf)
        occi = occupation(ir, i);
      else if (alp)
        occi = alpha_occupation(ir, i);
      else
        occi = beta_occupation(ir, i);

      if (fabs(occi-occ) < 1.0e-8) {
        if (col0 == -1)
          col0 = i;
        continue;
      } else if (col0 != -1) {
        coln = i-1;
        break;
      }
    }

    if (col0 == -1)
      continue;

    if (coln == -1)
      coln = n_orthoSO-1;
    
    if (local_ || local_dens_) {
      RefSymmSCMatrix ldir = dir;

      RefSCMatrix occbits; // holds the occupied bits of the scf vector

      // get local copies of vector and density matrix
      if (!local_) {
        Ref<SCMatrixKit> rk = new ReplSCMatrixKit;
        RefSCMatrix lvir = rk->matrix(vir.rowdim(), vir.coldim());
        lvir->convert(vir);
        occbits = lvir->get_subblock(0, n_SO-1, col0, coln);
        lvir = 0;

        ldir = rk->symmmatrix(dir.dim());
        ldir->convert(dir);

      } else {
        occbits = vir->get_subblock(0, n_SO-1, col0, coln);
      }
    
      double **c;
      double *dens;
      
      if (dynamic_cast<LocalSCMatrix*>(occbits.pointer()))
        c = dynamic_cast<LocalSCMatrix*>(occbits.pointer())->get_rows();
      else if (dynamic_cast<ReplSCMatrix*>(occbits.pointer()))
        c = dynamic_cast<ReplSCMatrix*>(occbits.pointer())->get_rows();
      else
        abort();

      if (dynamic_cast<LocalSymmSCMatrix*>(ldir.pointer()))
        dens = dynamic_cast<LocalSymmSCMatrix*>(ldir.pointer())->get_data();
      else if (dynamic_cast<ReplSymmSCMatrix*>(ldir.pointer()))
        dens = dynamic_cast<ReplSymmSCMatrix*>(ldir.pointer())->get_data();
      else
        abort();

      int ij=0;
      for (i=0; i < n_SO; i++) {
        for (j=0; j <= i; j++, ij++) {
          if (ij%nproc != me)
            continue;
          
          double dv = 0;

          int kk=0;
          for (k=col0; k <= coln; k++, kk++)
            dv += c[i][kk]*c[j][kk];

          dens[ij] = dv;
        }
      }

      if (nproc > 1)
        scf_grp_->sum(dens, i_offset(n_SO));

      if (!local_) {
        dir->convert(ldir);
      }
    }

    // for now quit
    else {
      ExEnv::err0() << indent
           << "Cannot yet use anything but Local matrices"
           << endl;
      abort();
    }
  }

  if (debug_ > 0) {
    ExEnv::out0() << indent
         << "Nelectron = " << 2.0 * (d * overlap()).trace() << endl;
  }

  int use_alternate_density = 0;
  if (use_alternate_density || debug_ > 2) {
    // double check the density with this simpler, slower way to compute
    // the density matrix
    RefSymmSCMatrix occ(oso_dimension(), basis_matrixkit());
    occ.assign(0.0);
    for (i=0; i<oso_dimension()->n(); i++) {
      occ(i,i) = occupation(i);
    }
    occ.scale(0.5);
    RefSymmSCMatrix d2(so_dimension(), basis_matrixkit());
    d2.assign(0.0);
    d2.accumulate_transform(vector, occ);
    if (debug_ > 2) {
      d2.print("d2 density");
      ExEnv::out0() << indent << "d2 Nelectron = "
                   << 2.0 * (d2 * overlap()).trace() << endl;
    }
    if (use_alternate_density) {
      d.assign(d2);
      ExEnv::out0() << indent
                   << "using alternate density algorithm" << endl;
    }
  }

  if (debug_ > 1) {
    d.print("SO Density");
    RefSCMatrix rd(d.dim(), d.dim(), basis_matrixkit());
    rd.assign(0.0);
    rd.accumulate(d);
    (d*overlap()*d-rd).print("SO Density idempotency error");
  }
}

double
SCF::one_body_energy()
{
  RefSymmSCMatrix dens = ao_density().copy();
  RefSymmSCMatrix hcore = dens->clone();
  hcore.assign(0.0);
  Ref<SCElementOp> hcore_op = new OneBodyIntOp(integral()->hcore());
  hcore.element_op(hcore_op);

  dens->scale_diagonal(0.5);
  SCElementScalarProduct *prod = new SCElementScalarProduct;
  prod->reference();
  Ref<SCElementOp2> op = prod;
  hcore->element_op(prod, dens);
  double e = prod->result();
  op = 0;
  prod->dereference();
  delete prod;
  return 2.0 * e;
}

void
SCF::two_body_energy(double &ec, double &ex)
{
  ExEnv::errn() << class_name() << ": two_body_energy not implemented" << endl;
}

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

void
SCF::init_threads()
{
  int nthread = threadgrp_->nthread();
  size_t int_store = integral()->storage_unused()/nthread;
  
  // initialize the two electron integral classes
  tbis_ = new Ref<TwoBodyInt>[nthread];
  for (int i=0; i < nthread; i++) {
    tbis_[i] = integral()->electron_repulsion();
    tbis_[i]->set_integral_storage(int_store);
  }

}

void
SCF::done_threads()
{
  for (int i=0; i < threadgrp_->nthread(); i++) tbis_[i] = 0;
  delete[] tbis_;
  tbis_ = 0;
}

int *
SCF::read_occ(const Ref<KeyVal> &keyval, const char *name, int nirrep)
{
  int *occ = 0;
  if (keyval->exists(name)) {
    if (keyval->count(name) != nirrep) {
      ExEnv::err0() << indent
                   << "ERROR: SCF: have " << nirrep << " irreps but "
                   << name << " vector is length " << keyval->count(name)
                   << endl;
      abort();
    }
    occ = new int[nirrep];
    for (int i=0; i<nirrep; i++) {
      occ[i] = keyval->intvalue(name,i);
    }
  }
  return occ;
}

void
SCF::obsolete()
{
  OneBodyWavefunction::obsolete();
  if (guess_wfn_.nonnull()) guess_wfn_->obsolete();
}

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

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