obwfn.cc

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

  }

#if DEBUG
  val.print("projected values");
#endif

  return val;
}

RefSCMatrix
OneBodyWavefunction::so_to_mo()
{
  return orthog_so_to_mo() * so_to_orthog_so();
}

RefSCMatrix
OneBodyWavefunction::orthog_so_to_mo()
{
  return oso_eigenvectors().t();
}

RefSCMatrix
OneBodyWavefunction::mo_to_so()
{
  return so_to_orthog_so_inverse() * mo_to_orthog_so();
}

RefSCMatrix
OneBodyWavefunction::mo_to_orthog_so()
{
  return oso_eigenvectors();
}

RefSCMatrix
OneBodyWavefunction::eigenvectors()
{
  return so_to_orthog_so().t() * oso_eigenvectors();
}

RefSCMatrix
OneBodyWavefunction::hcore_guess()
{
  RefDiagSCMatrix val;
  return hcore_guess(val);
}

RefSCMatrix
OneBodyWavefunction::hcore_guess(RefDiagSCMatrix &val)
{
  RefSCMatrix vec(oso_dimension(), oso_dimension(), basis_matrixkit());
  val = basis_matrixkit()->diagmatrix(oso_dimension());

  // I'm about to do something strange, but it will only work
  // if the SO and orthogonal SO dimensions are equivalent.  This
  // is not the case for canonical orthogonalization when there
  // are linear dependencies.
  if (so_dimension()->equiv(oso_dimension())) {
    // Yes, this is diagonalizing Hcore in a nonorthogonal basis
    // and does not really make any sense--except it seems to
    // always give a better initial guess.  I don't understand
    // why it works better.
    core_hamiltonian().diagonalize(val,vec);
  }
  else {
    RefSymmSCMatrix hcore_oso(oso_dimension(), basis_matrixkit());
    hcore_oso->assign(0.0);
    hcore_oso->accumulate_transform(so_to_orthog_so(), core_hamiltonian());

    if (debug_ > 1) {
      hcore_oso.print("hcore in ortho SO basis");
    }

    hcore_oso.diagonalize(val,vec);

    if (debug_ > 1) {
      val.print("hcore eigenvalues in ortho SO basis");
      vec.print("hcore eigenvectors in ortho SO basis");
    }
  }

  return vec;
}

// Function for returning an orbital value at a point
double
OneBodyWavefunction::orbital(const SCVector3& r, int iorb)
{
  return Wavefunction::orbital(r,iorb,eigenvectors());
}

// Function for returning an orbital value at a point
double
OneBodyWavefunction::orbital_density(const SCVector3& r,
                                            int iorb,
                                            double* orbval)
{
  return Wavefunction::orbital_density(r,iorb,eigenvectors(),orbval);
}

void
OneBodyWavefunction::print(ostream&o) const
{
  Wavefunction::print(o);
}

void
OneBodyWavefunction::init_sym_info()
{
  RefSCDimension d = oso_dimension();
  nirrep_ = d->blocks()->nblock();
  nvecperirrep_ = new int[nirrep_];
  occupations_ = new double[d->n()];
  alpha_occupations_ = new double[d->n()];
  beta_occupations_ = new double[d->n()];

  int ij=0;
  for (int i=0; i < nirrep_; i++) {
    nvecperirrep_[i] = d->blocks()->size(i);

    for (int j=0; j < nvecperirrep_[i]; j++, ij++) {
      if (!spin_unrestricted()) occupations_[ij] = occupation(i,j);
      else occupations_[ij] = 0.0;
      alpha_occupations_[ij] = alpha_occupation(i,j);
      beta_occupations_[ij] = beta_occupation(i,j);
    }
  }
}

double
OneBodyWavefunction::occupation(int vectornum)
{
  if (spin_unrestricted()) {
    ExEnv::errn() << "OneBodyWavefunction::occupation: called for USCF case"
                 << endl;
    abort();
  }
  if (!nirrep_) init_sym_info();
  return occupations_[vectornum];
}

double
OneBodyWavefunction::alpha_occupation(int vectornum)
{
  if (!nirrep_) init_sym_info();
  return alpha_occupations_[vectornum];
}

double
OneBodyWavefunction::beta_occupation(int vectornum)
{
  if (!nirrep_) init_sym_info();
  return beta_occupations_[vectornum];
}

double
OneBodyWavefunction::alpha_occupation(int irrep, int vectornum)
{
  if (!spin_polarized())
    return 0.5*occupation(irrep, vectornum);
  
  ExEnv::errn() << class_name() << "::alpha_occupation not implemented" << endl;
  abort();
  return 0;
}

double
OneBodyWavefunction::beta_occupation(int irrep, int vectornum)
{
  if (!spin_polarized())
    return 0.5*occupation(irrep, vectornum);
  
  ExEnv::errn() << class_name() << "::beta_occupation not implemented" << endl;
  abort();
  return 0;
}

RefSCMatrix
OneBodyWavefunction::oso_alpha_eigenvectors()
{
  if (!spin_unrestricted())
    return oso_eigenvectors().copy();

  ExEnv::errn() << class_name() << "::oso_alpha_eigenvectors not implemented" << endl;
  abort();
  return 0;
}

RefSCMatrix
OneBodyWavefunction::oso_beta_eigenvectors()
{
  if (!spin_unrestricted())
    return oso_eigenvectors().copy();

  ExEnv::errn() << class_name() << "::oso_beta_eigenvectors not implemented" << endl;
  abort();
  return 0;
}

RefSCMatrix
OneBodyWavefunction::alpha_eigenvectors()
{
  if (!spin_unrestricted())
    return eigenvectors().copy();

  ExEnv::errn() << class_name() << "::alpha_eigenvectors not implemented" << endl;
  abort();
  return 0;
}

RefSCMatrix
OneBodyWavefunction::beta_eigenvectors()
{
  if (!spin_unrestricted())
    return eigenvectors().copy();

  ExEnv::errn() << class_name() << "::beta_eigenvectors not implemented" << endl;
  abort();
  return 0;
}

RefDiagSCMatrix
OneBodyWavefunction::alpha_eigenvalues()
{
  if (!spin_unrestricted())
    return eigenvalues().copy();

  ExEnv::errn() << class_name() << "::alpha_eigenvalues not implemented" << endl;
  abort();
  return 0;
}

RefDiagSCMatrix
OneBodyWavefunction::beta_eigenvalues()
{
  if (!spin_unrestricted())
    return eigenvalues().copy();

  ExEnv::errn() << class_name() << "::beta_eigenvalues not implemented" << endl;
  abort();
  return 0;
}

int
OneBodyWavefunction::nelectron()
{
  int noso = oso_dimension()->n();
  double tocc = 0.0;
  if (!spin_polarized()) {
    for (int i=0; i<noso; i++) {
      tocc += occupation(i);
    }
  }
  else {
    for (int i=0; i<noso; i++) {
      tocc += alpha_occupation(i) + beta_occupation(i);
    }
  }
  return int(tocc+0.5);
}

void
OneBodyWavefunction::symmetry_changed()
{
  Wavefunction::symmetry_changed();

  // junk the old occupation information
  delete[] nvecperirrep_;
  delete[] occupations_;
  delete[] alpha_occupations_;
  delete[] beta_occupations_;
  nirrep_ = 0;
  nvecperirrep_=0;
  occupations_=0;
  alpha_occupations_=0;
  beta_occupations_=0;
  
  // for now, delete old eigenvectors...later we'll transform to new
  // pointgroup
  oso_eigenvectors_.result_noupdate() = 0;
}

int
OneBodyWavefunction::form_occupations(int *&newocc, const int *oldocc)
{
  delete[] newocc;
  newocc = 0;

  CorrelationTable corrtab;
  if (corrtab.initialize_table(initial_pg_, molecule()->point_group()))
    return 0;

  newocc = new int[corrtab.subn()];
  memset(newocc,0,sizeof(int)*corrtab.subn());

  for (int i=0; i<corrtab.n(); i++) {
      for (int j=0; j<corrtab.ngamma(i); j++) {
          int gam = corrtab.gamma(i,j);
          newocc[gam] += (corrtab.subdegen(gam)*oldocc[i])/corrtab.degen(i);
        }
    }

  return 1;
}
  
/////////////////////////////////////////////////////////////////////////////

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