uscf.cc

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

//
// uscf.cc --- implementation of the UnrestrictedSCF abstract base class
//
// Copyright (C) 1997 Limit Point Systems, Inc.
//
// Author: Edward Seidl <[e-mail]>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//

#ifdef __GNUC__
#pragma implementation
#endif

#include <math.h>
#include <limits.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

#include <util/state/stateio.h>

#include <util/misc/timer.h>
#include <util/misc/formio.h>

#include <math/scmat/local.h>
#include <math/scmat/repl.h>
#include <math/scmat/offset.h>
#include <math/scmat/block.h>
#include <math/scmat/blocked.h>
#include <math/scmat/blkiter.h>
#include <math/scmat/local.h>

#include <math/optimize/scextrapmat.h>
#include <math/optimize/diis.h>

#include <chemistry/qc/basis/petite.h>
#include <chemistry/qc/scf/scfops.h>
#include <chemistry/qc/scf/scflocal.h>
#include <chemistry/qc/scf/uscf.h>
#include <chemistry/qc/scf/ltbgrad.h>
#include <chemistry/qc/scf/uhftmpl.h>

using namespace std;
using namespace sc;

namespace sc {

///////////////////////////////////////////////////////////////////////////
// UnrestrictedSCF

static ClassDesc UnrestrictedSCF_cd(
  typeid(UnrestrictedSCF),"UnrestrictedSCF",2,"public SCF",
  0, 0, 0);

UnrestrictedSCF::UnrestrictedSCF(StateIn& s) :
  SavableState(s),
  SCF(s),
  oso_eigenvectors_beta_(this),
  eigenvalues_beta_(this),
  focka_(this),
  fockb_(this)
{
  need_vec_ = 1;
  compute_guess_ = 0;

  oso_eigenvectors_beta_.result_noupdate() =
    basis_matrixkit()->matrix(so_dimension(), oso_dimension());
  oso_eigenvectors_beta_.restore_state(s);
  oso_eigenvectors_beta_.result_noupdate().restore(s);

  eigenvalues_beta_.result_noupdate() =
    basis_matrixkit()->diagmatrix(oso_dimension());
  eigenvalues_beta_.restore_state(s);
  eigenvalues_beta_.result_noupdate().restore(s);

  focka_.result_noupdate() =
    basis_matrixkit()->symmmatrix(so_dimension());
  focka_.restore_state(s);
  focka_.result_noupdate().restore(s);

  fockb_.result_noupdate() =
    basis_matrixkit()->symmmatrix(so_dimension());
  fockb_.restore_state(s);
  fockb_.result_noupdate().restore(s);

  s.get(user_occupations_);
  s.get(tnalpha_);
  s.get(tnbeta_);
  s.get(nirrep_);
  s.get(nalpha_);
  s.get(nbeta_);

  if (s.version(::class_desc<UnrestrictedSCF>()) >= 2) {
    s.get(initial_nalpha_);
    s.get(initial_nbeta_);
    most_recent_pg_ << SavableState::restore_state(s);
  } else {
    initial_nalpha_ = new int[nirrep_];
    memcpy(initial_nalpha_, nalpha_, sizeof(int)*nirrep_);
    initial_nbeta_ = new int[nirrep_];
    memcpy(initial_nbeta_, nbeta_, sizeof(int)*nirrep_);
  }

  init_mem(4);
}

UnrestrictedSCF::UnrestrictedSCF(const Ref<KeyVal>& keyval) :
  SCF(keyval),
  oso_eigenvectors_beta_(this),
  eigenvalues_beta_(this),
  focka_(this),
  fockb_(this)
{
  int i;
  
  double acc = oso_eigenvectors_.desired_accuracy();
  oso_eigenvectors_beta_.set_desired_accuracy(acc);
  eigenvalues_beta_.set_desired_accuracy(acc);

  if (oso_eigenvectors_beta_.desired_accuracy() < DBL_EPSILON) {
    oso_eigenvectors_beta_.set_desired_accuracy(DBL_EPSILON);
    eigenvalues_beta_.set_desired_accuracy(DBL_EPSILON);
  }

  focka_.compute()=0;
  focka_.computed()=0;
  fockb_.compute()=0;
  fockb_.computed()=0;

  // calculate the total nuclear charge
  double Znuc=molecule()->nuclear_charge();

  // check to see if this is to be a charged molecule
  double charge = keyval->doublevalue("total_charge");
  int nelectrons = (int)(Znuc-charge+1.0e-4);

  // first let's try to figure out how many open shells there are
  if (keyval->exists("multiplicity")) {
    int mult = keyval->intvalue("multiplicity");
    if (mult < 1) {
      ExEnv::err0() << endl << indent
           << "USCF::init: bad value for multiplicity: " << mult << endl
           << indent << "assuming singlet" << endl;
      mult=1;
    }
    
    // for singlet, triplet, etc. we need an even number of electrons
    // for doublet, quartet, etc. we need an odd number of electrons
    if ((mult%2 && nelectrons%2) || (!(mult%2) && !(nelectrons%2))) {
      ExEnv::err0() << endl << indent
           << "USCF::init: Warning, there's a leftover electron..."
           << " I'm going to get rid of it" << endl
           << incindent << indent << "total_charge = " << charge << endl
           << indent << "total nuclear charge = " << Znuc << endl
           << indent << "multiplicity = " << mult << endl << decindent;
      nelectrons--;
    }
    if (mult%2)
      tnalpha_ = nelectrons/2 + (mult-1)/2;
    else
      tnalpha_ = nelectrons/2 + mult/2;

  } else {
    // if there's an odd number of electrons, then do a doublet, otherwise
    // do a triplet
    tnalpha_=nelectrons/2+1;
  }

  tnbeta_ = nelectrons-tnalpha_;
  
  ExEnv::out0() << endl << indent
       << "USCF::init: total charge = " << Znuc-tnalpha_-tnbeta_
       << endl << endl;

  nirrep_ = molecule()->point_group()->char_table().ncomp();

  nalpha_ = read_occ(keyval, "alpha", nirrep_);
  nbeta_ = read_occ(keyval, "beta", nirrep_);
  if (nalpha_ && nbeta_) {
    tnalpha_ = 0;
    tnbeta_ = 0;
    user_occupations_=1;
    for (i=0; i < nirrep_; i++) {
      tnalpha_ += nalpha_[i];
      tnbeta_ += nbeta_[i];
    }
    initial_nalpha_ = new int[nirrep_];
    memcpy(initial_nalpha_, nalpha_, sizeof(int)*nirrep_);
    initial_nbeta_ = new int[nirrep_];
    memcpy(initial_nbeta_, nbeta_, sizeof(int)*nirrep_);
  }
  else if (nalpha_ && !nbeta_ || !nalpha_ && nbeta_) {
    ExEnv::out0() << "ERROR: USCF: only one of alpha and beta specified: "
                 << "give both or none" << endl;
    abort();
  }
  else {
    initial_nalpha_=0;
    initial_nbeta_=0;
    nalpha_=0;
    nbeta_=0;
    user_occupations_=0;
    set_occupations(0,0);
  }

  ExEnv::out0() << indent << "alpha = [";
  for (i=0; i < nirrep_; i++)
    ExEnv::out0() << " " << nalpha_[i];
  ExEnv::out0() << " ]\n";

  ExEnv::out0() << indent << "beta  = [";
  for (i=0; i < nirrep_; i++)
    ExEnv::out0() << " " << nbeta_[i];
  ExEnv::out0() << " ]\n";

  // check to see if this was done in SCF(keyval)
  if (!keyval->exists("maxiter"))
    maxiter_ = 100;

  if (!keyval->exists("level_shift"))
    level_shift_ = 0.25;

  // now take care of memory stuff
  init_mem(4);
}

UnrestrictedSCF::~UnrestrictedSCF()
{
  if (nalpha_) {
    delete[] nalpha_;
    nalpha_=0;
  }
  if (nbeta_) {
    delete[] nbeta_;
    nbeta_=0;
  }
  delete[] initial_nalpha_;
  delete[] initial_nbeta_;
}

void
UnrestrictedSCF::save_data_state(StateOut& s)
{
  SCF::save_data_state(s);
  oso_eigenvectors_beta_.save_data_state(s);
  oso_eigenvectors_beta_.result_noupdate().save(s);
  eigenvalues_beta_.save_data_state(s);
  eigenvalues_beta_.result_noupdate().save(s);
  focka_.save_data_state(s);
  focka_.result_noupdate().save(s);
  fockb_.save_data_state(s);
  fockb_.result_noupdate().save(s);

  s.put(user_occupations_);
  s.put(tnalpha_);
  s.put(tnbeta_);
  s.put(nirrep_);
  s.put(nalpha_, nirrep_);
  s.put(nbeta_, nirrep_);

  s.put(initial_nalpha_,initial_pg_->char_table().ncomp());
  s.put(initial_nbeta_,initial_pg_->char_table().ncomp());
  SavableState::save_state(most_recent_pg_.pointer(),s);
}

double
UnrestrictedSCF::occupation(int ir, int i)
{
  abort();
  return 0;
}

double
UnrestrictedSCF::alpha_occupation(int ir, int i)
{
  if (i < nalpha_[ir]) return 1.0;
  return 0.0;
}

double
UnrestrictedSCF::beta_occupation(int ir, int i)
{
  if (i < nbeta_[ir]) return 1.0;
  return 0.0;
}

RefSCMatrix
UnrestrictedSCF::eigenvectors()
{
  abort();
  return 0;
}

RefDiagSCMatrix
UnrestrictedSCF::eigenvalues()
{
  abort();
  return 0;
}

RefSCMatrix
UnrestrictedSCF::oso_alpha_eigenvectors()
{
  return oso_eigenvectors_.result();
}

RefSCMatrix
UnrestrictedSCF::alpha_eigenvectors()
{
  return so_to_orthog_so().t() * oso_eigenvectors_.result();
}

RefDiagSCMatrix
UnrestrictedSCF::alpha_eigenvalues()
{
  return eigenvalues_.result();
}

RefSCMatrix
UnrestrictedSCF::oso_beta_eigenvectors()
{
  return oso_eigenvectors_beta_.result();
}

RefSCMatrix
UnrestrictedSCF::beta_eigenvectors()
{
  return so_to_orthog_so().t() * oso_eigenvectors_beta_.result();
}

RefDiagSCMatrix
UnrestrictedSCF::beta_eigenvalues()
{
  return eigenvalues_beta_.result();
}

int
UnrestrictedSCF::spin_polarized()
{
  return 1;
}

int
UnrestrictedSCF::spin_unrestricted()
{
  return 1;
}

int
UnrestrictedSCF::n_fock_matrices() const
{
  return 2;
}

RefSymmSCMatrix
UnrestrictedSCF::fock(int n)
{
  if (n > 1) {
    ExEnv::err0() << indent
         << "USCF::fock: there are only two fock matrices, "
         << scprintf("but fock(%d) was requested\n",n);
    abort();
  }

  if (n==0)
    return focka_.result();
  else
    return fockb_.result();
}

void
UnrestrictedSCF::print(ostream&o) const
{
  int i;

  SCF::print(o);
  o << indent << "UnrestrictedSCF Parameters:\n" << incindent
    << indent << "charge = " << molecule()->nuclear_charge()
                                - tnalpha_ - tnbeta_ << endl
    << indent << "nalpha = " << tnalpha_ << endl
    << indent << "nbeta = " << tnbeta_ << endl
    << indent << "alpha = [";

  for (i=0; i < nirrep_; i++)
    o << " " << nalpha_[i];
  o << " ]" << endl;

  o << indent << "beta  = [";
  for (i=0; i < nirrep_; i++)
    o << " " << nbeta_[i];
  o << " ]" << endl << decindent << endl;
}

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

void
UnrestrictedSCF::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.
      // right now the check is crude...there should be an equiv member in
      // GaussianBasisSet
      if (guess_wfn_.nonnull()) {
        if (guess_wfn_->basis()->nbasis() == basis()->nbasis()) {
          ExEnv::out0() << indent
               << "Using guess wavefunction as starting vector" << endl;

          // indent output of eigenvectors() call if there is any
          ExEnv::out0() << incindent << incindent;
          UnrestrictedSCF *ug =
            dynamic_cast<UnrestrictedSCF*>(guess_wfn_.pointer());
          if (!ug || compute_guess_) {
            oso_eigenvectors_ = guess_wfn_->oso_alpha_eigenvectors().copy();
            eigenvalues_ = guess_wfn_->alpha_eigenvalues().copy();