mbpt.cc

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

{
  // Solaris 2.7 workshop 5.0 is causing this routine to
  // be incorrectly called in a base class CTOR.  Thus
  // reference_ might be null and it must be tested.
  if (reference_.nonnull()) reference_->obsolete();
  Wavefunction::obsolete();
}

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

int
MBPT2::gradient_implemented() const
{
  int nb = reference_->oso_dimension()->n();
  int n = 0;

  for (int i=0; i<nb; i++) {
    if (reference_->occupation(i) == 1.0) n++;
    }

  if (n) return 0;
  return 1;
}

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

int
MBPT2::value_implemented() const
{
  return 1;
}

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

void
MBPT2::eigen(RefDiagSCMatrix &vals, RefSCMatrix &vecs, RefDiagSCMatrix &occs)
{
  int i, j;
  if (nsocc) {
      if (reference_->n_fock_matrices() != 2) {
          ExEnv::errn() << "MBPT2: given open reference with"
               << " wrong number of Fock matrices" << endl;
          abort();
        }

      // Notation: oo = orthonormal symmetry orbital basis
      //           ao = atomic orbital basis
      //           so = symmetrized atomic orbital basis
      //           mo1 = SCF molecular orbital basis
      //           mo2 = MBPT molecular orbital basis

      // get the closed shell and open shell AO fock matrices
      RefSymmSCMatrix fock_c_so = reference_->fock(0);
      RefSymmSCMatrix fock_o_so = reference_->fock(1);

      // transform the AO fock matrices to the MO basis
      RefSymmSCMatrix fock_c_mo1
          = basis_matrixkit()->symmmatrix(oso_dimension());
      RefSymmSCMatrix fock_o_mo1
          = basis_matrixkit()->symmmatrix(oso_dimension());
      RefSCMatrix vecs_so_mo1 = reference_->eigenvectors();

      fock_c_mo1.assign(0.0);
      fock_o_mo1.assign(0.0);
      fock_c_mo1.accumulate_transform(vecs_so_mo1.t(), fock_c_so);
      fock_o_mo1.accumulate_transform(vecs_so_mo1.t(), fock_o_so);
      fock_c_so = 0;
      fock_o_so = 0;

     /* Convert to the Guest & Saunders general form.
        This is the form used for an OPT2 calculation.

            C        O         V
        ----------
        |        |
     C  |   fc   |
        |        |
        -------------------
        |        |        |
     O  | 2fc-fo |   fc   |
        |        |        |
        ----------------------------
        |        |        |        |
     V  |   fc   |   fo   |   fc   |
        |        |        |        |
        ----------------------------
      */
      RefSymmSCMatrix fock_eff_mo1 = fock_c_mo1.clone();
      fock_eff_mo1.assign(fock_c_mo1);
      for (i=0; i<oso_dimension()->n(); i++) {
          double occi = reference_->occupation(i);
          for (j=0; j<=i; j++) {
              double occj = reference_->occupation(j);
              if (occi == 2.0 && occj == 1.0
                  || occi == 1.0 && occj == 2.0) {
                  fock_eff_mo1.accumulate_element(i,j,
                                                  fock_c_mo1(i,j)-fock_o_mo1(i,j));
                }
              else if (occi == 0.0 && occj == 1.0
                  || occi == 1.0 && occj == 0.0) {
                  fock_eff_mo1.accumulate_element(i,j,
                                                  fock_o_mo1(i,j)-fock_c_mo1(i,j));
                }
            }
        }

      // diagonalize the new fock matrix
      RefDiagSCMatrix vals_mo2(fock_eff_mo1.dim(), fock_eff_mo1.kit());
      RefSCMatrix vecs_mo1_mo2(fock_eff_mo1.dim(), fock_eff_mo1.dim(),
                               fock_eff_mo1.kit());
      fock_eff_mo1.diagonalize(vals_mo2, vecs_mo1_mo2);
      vals = vals_mo2;

      // compute the AO to new MO scf vector
      RefSCMatrix so_ao = reference_->integral()->petite_list()->sotoao();

      if (debug_ > 1) {
          vecs_mo1_mo2.t().print("vecs_mo1_mo2.t()");
          vecs_so_mo1.t().print("vecs_so_mo1.t()");
          so_ao.print("so_ao");
        }
      
      vecs = vecs_mo1_mo2.t() * vecs_so_mo1.t() * so_ao;
    }
  else {
      if (debug_) ExEnv::out0() << indent << "getting fock matrix" << endl;
      // get the closed shell AO fock matrices
      RefSymmSCMatrix fock_c_so = reference_->fock(0);

      // transform the AO fock matrices to the MO basis
      RefSymmSCMatrix fock_c_mo1
          = basis_matrixkit()->symmmatrix(oso_dimension());
      RefSCMatrix vecs_so_mo1 = reference_->eigenvectors();

      fock_c_mo1.assign(0.0);
      fock_c_mo1.accumulate_transform(vecs_so_mo1.t(), fock_c_so);
      fock_c_so = 0;

      if (debug_) ExEnv::out0() << indent << "diagonalizing" << endl;
      // diagonalize the fock matrix
      vals = fock_c_mo1.eigvals();

      // compute the AO to new MO scf vector
      if (debug_) ExEnv::out0() << indent << "AO to MO" << endl;
      RefSCMatrix so_ao = reference_->integral()->petite_list()->sotoao();
      vecs = vecs_so_mo1.t() * so_ao;
    }
  // fill in the occupations
  occs = matrixkit()->diagmatrix(vals.dim());
  for (i=0; i<oso_dimension()->n(); i++) {
      occs(i) = reference_->occupation(i);
    }
  // allocate storage for symmetry information
  if (!symorb_irrep_) symorb_irrep_ = new int[nbasis];
  if (!symorb_num_) symorb_num_ = new int[nbasis];
  // Check for degenerate eigenvalues.  Use unsorted eigenvalues since it
  // only matters if the degeneracies occur within a given irrep.
  BlockedDiagSCMatrix *bvals = dynamic_cast<BlockedDiagSCMatrix*>(vals.pointer());
  for (i=0; i<bvals->nblocks(); i++) {
      int done = 0;
      RefDiagSCMatrix valsi = bvals->block(i);
      for (j=1; j<valsi.n(); j++) {
          if (fabs(valsi(j)-valsi(j-1)) < 1.0e-7) {
              ExEnv::out0() << indent
                   << "NOTE: There are degenerate orbitals within an irrep."
                   << "  This will make"
                   << endl
                   << indent
                   << "      some diagnostics, such as the largest amplitude,"
                   << " nonunique."
                   << endl;
              done = 1;
              break;
            }
          if (done) break;
        }
    }
  // sort the eigenvectors and values if symmetry is not c1
  if (molecule()->point_group()->char_table().order() != 1) {
      if (debug_) ExEnv::out0() << indent << "sorting eigenvectors" << endl;
      double *evals = new double[noso];
      vals->convert(evals);
      int *indices = new int[noso];
      dquicksort(evals,indices,noso);
      delete[] evals;
      // make sure all nodes see the same indices and evals
      msg_->bcast(indices,noso);
      RefSCMatrix newvecs(vecs.rowdim(), vecs.coldim(), matrixkit());
      RefDiagSCMatrix newvals(vals.dim(), matrixkit());
      RefDiagSCMatrix newoccs(vals.dim(), matrixkit());
      for (i=0; i<noso; i++) {
          newvals(i) = vals(indices[i]);
          newoccs(i) = occs(indices[i]);
          for (j=0; j<nbasis; j++) {
              newvecs(i,j) = vecs(indices[i],j);
            }
        }
      occs = newoccs;
      vecs = newvecs;
      vals = newvals;

      // compute orbital symmetry information
      CharacterTable ct = molecule()->point_group()->char_table();
      int orbnum = 0;
      int *tmp_irrep = new int[noso];
      int *tmp_num = new int[noso];
      for (i=0; i<oso_dimension()->blocks()->nblock(); i++) {
          for (j=0; j<oso_dimension()->blocks()->size(i); j++, orbnum++) {
              tmp_irrep[orbnum] = i;
              tmp_num[orbnum] = j;
            }
        }
      for (i=0; i<noso; i++) {
          symorb_irrep_[i] = tmp_irrep[indices[i]];
          symorb_num_[i] = tmp_num[indices[i]];
        }
      delete[] tmp_irrep;
      delete[] tmp_num;

      delete[] indices;
    }
  else {
      // compute orbital symmetry information for c1
      for (i=0; i<noso; i++) {
          symorb_irrep_[i] = 0;
          symorb_num_[i] = i;
        }
    }
  // check the splitting between frozen and nonfrozen orbitals
  if (nfzc && nfzc < noso) {
      double split = vals(nfzc) - vals(nfzc-1);
      if (split < 0.2) {
          ExEnv::out0() << endl
               << indent << "WARNING: "
               << "MBPT2: gap between frozen and active occupied orbitals is "
               << split << " au" << endl << endl;
        }
    }
  if (nfzv && noso-nfzv-1 >= 0) {
      double split = vals(nbasis-nfzv) - vals(nbasis-nfzv-1);
      if (split < 0.2) {
          ExEnv::out0() << endl
               << indent << "WARNING: "
               << "MBPT2: gap between frozen and active virtual orbitals is "
               << split << " au" << endl << endl;
        }
    }
  if (debug_) ExEnv::out0() << indent << "eigen done" << endl;
}

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

void
MBPT2::init_variables()
{
  nbasis = so_dimension()->n();
  noso = oso_dimension()->n();
//    if (nbasis != noso) {
//        ExEnv::outn() << "MBPT2: Noso != Nbasis: MBPT2 not checked for this case" << endl;
//        abort();
//      }
  nocc = nvir = nsocc = 0;
  for (int i=0; i<noso; i++) {
    if (reference_->occupation(i) == 2.0) nocc++;
    else if (reference_->occupation(i) == 1.0) nsocc++;
    else nvir++;
    }
}

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

void
MBPT2::symmetry_changed()
{
  Wavefunction::symmetry_changed();
  reference_->symmetry_changed();
}

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

int
MBPT2::nelectron()
{
  return reference_->nelectron();
}

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

double
MBPT2::ref_energy()
{
  return reference_->energy();
}

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

double
MBPT2::corr_energy()
{
  return energy() - reference_->energy();
}

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

RefSCVector
MBPT2::ref_energy_gradient()
{
  gradient();
  return hf_gradient_;
}

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

RefSCVector
MBPT2::corr_energy_gradient()
{
  gradient();
  return get_cartesian_gradient() - hf_gradient_;
}

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

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