scf.cc

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

//
// scf.cc --- implementation of the SCF abstract base class
//
// Copyright (C) 1996 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/misc/formio.h>
#include <util/state/stateio.h>
#include <util/group/mstate.h>

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

#include <math/optimize/diis.h>

#include <chemistry/qc/basis/petite.h>
#include <chemistry/qc/scf/scf.h>

using namespace std;
using namespace sc;

///////////////////////////////////////////////////////////////////////////
// SCF

static ClassDesc SCF_cd(
  typeid(SCF),"SCF",6,"public OneBodyWavefunction",
  0, 0, 0);

SCF::SCF(StateIn& s) :
  SavableState(s),
  OneBodyWavefunction(s)
{
  need_vec_ = 1;
  compute_guess_ = 0;

  s.get(maxiter_,"maxiter");
  s.get(dens_reset_freq_);
  s.get(reset_occ_);
  s.get(local_dens_);
  if (s.version(::class_desc<SCF>()) >= 3) {
    double dstorage;
    s.get(dstorage);
    storage_ = size_t(dstorage);
  }
  else {
    unsigned int istorage;
    s.get(istorage);
    storage_ = istorage;
  }
  if (s.version(::class_desc<SCF>()) >= 2) {
    s.get(print_all_evals_);
    s.get(print_occ_evals_);
  }
  else {
    print_all_evals_ = 0;
    print_occ_evals_ = 0;
  }
  s.get(level_shift_);
  if (s.version(::class_desc<SCF>()) >= 5) {
    s.get(keep_guess_wfn_);
    guess_wfn_ << SavableState::restore_state(s);
  }
  else keep_guess_wfn_ = 0;
  if (s.version(::class_desc<SCF>()) >= 6) {
    s.get(always_use_guess_wfn_);
  }
  else always_use_guess_wfn_ = 0;

  extrap_ << SavableState::restore_state(s);
  accumdih_ << SavableState::restore_state(s);
  accumddh_ << SavableState::restore_state(s);

  scf_grp_ = basis()->matrixkit()->messagegrp();
  threadgrp_ = ThreadGrp::get_default_threadgrp();
}

SCF::SCF(const Ref<KeyVal>& keyval) :
  OneBodyWavefunction(keyval),
  need_vec_(1),
  compute_guess_(0),
  maxiter_(40),
  dens_reset_freq_(10),
  reset_occ_(0),
  local_dens_(1),
  storage_(0),
  level_shift_(0)
{
  if (keyval->exists("maxiter"))
    maxiter_ = keyval->intvalue("maxiter");

  if (keyval->exists("density_reset_frequency"))
    dens_reset_freq_ = keyval->intvalue("density_reset_frequency");

  if (keyval->exists("reset_occupations"))
    reset_occ_ = keyval->booleanvalue("reset_occupations");

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

  extrap_ << keyval->describedclassvalue("extrap");
  if (extrap_.null())
    extrap_ = new DIIS;

  accumdih_ << keyval->describedclassvalue("accumdih");
  if (accumdih_.null())
    accumdih_ = new AccumHNull;
  
  accumddh_ << keyval->describedclassvalue("accumddh");
  if (accumddh_.null())
    accumddh_ = new AccumHNull;
  
  KeyValValuesize defaultmem(DEFAULT_SC_MEMORY);
  storage_ = keyval->sizevalue("memory",defaultmem);
  
  if (keyval->exists("local_density"))
    local_dens_ = keyval->booleanvalue("local_density");
    
  print_all_evals_ = keyval->booleanvalue("print_evals");
  print_occ_evals_ = keyval->booleanvalue("print_occupied_evals");
  
  scf_grp_ = basis()->matrixkit()->messagegrp();
  threadgrp_ = ThreadGrp::get_default_threadgrp();
  
  keep_guess_wfn_ = keyval->booleanvalue("keep_guess_wavefunction");

  always_use_guess_wfn_
    = keyval->booleanvalue("always_use_guess_wavefunction");

  // first see if guess_wavefunction is a wavefunction, then check to
  // see if it's a string.
  if (keyval->exists("guess_wavefunction")) {
    ExEnv::out0() << incindent << incindent;
    guess_wfn_ << keyval->describedclassvalue("guess_wavefunction");
    compute_guess_=1;
    if (guess_wfn_.null()) {
      compute_guess_=0;
      char *path = keyval->pcharvalue("guess_wavefunction");
      struct stat sb;
      if (path && stat(path, &sb)==0 && sb.st_size) {
        BcastStateInBin s(scf_grp_, path);

        // reset the default matrixkit so that the matrices in the guess
        // wavefunction will match those in this wavefunction
        Ref<SCMatrixKit> oldkit = SCMatrixKit::default_matrixkit();
        SCMatrixKit::set_default_matrixkit(basis()->matrixkit());

        guess_wfn_ << SavableState::restore_state(s);

        // go back to the original default matrixkit
        SCMatrixKit::set_default_matrixkit(oldkit);
        delete[] path;
      }
    }
    ExEnv::out0() << decindent << decindent;
  }
}

SCF::~SCF()
{
}

void
SCF::save_data_state(StateOut& s)
{
  OneBodyWavefunction::save_data_state(s);
  s.put(maxiter_);
  s.put(dens_reset_freq_);
  s.put(reset_occ_);
  s.put(local_dens_);
  double dstorage = storage_;
  s.put(dstorage);
  s.put(print_all_evals_);
  s.put(print_occ_evals_);
  s.put(level_shift_);
  s.put(keep_guess_wfn_);
  SavableState::save_state(guess_wfn_.pointer(),s);
  s.put(always_use_guess_wfn_);
  SavableState::save_state(extrap_.pointer(),s);
  SavableState::save_state(accumdih_.pointer(),s);
  SavableState::save_state(accumddh_.pointer(),s);
}

RefSCMatrix
SCF::oso_eigenvectors()
{
  return oso_eigenvectors_.result();
}

RefDiagSCMatrix
SCF::eigenvalues()
{
  return eigenvalues_.result();
}

int
SCF::spin_unrestricted()
{
  return 0;
}

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

void
SCF::print(ostream&o) const
{
  OneBodyWavefunction::print(o);
  o << indent << "SCF Parameters:\n" << incindent
    << indent << "maxiter = " << maxiter_ << endl
    << indent << "density_reset_frequency = " << dens_reset_freq_ << endl
    << indent << scprintf("level_shift = %f\n",level_shift_) 
    << decindent << endl;
}

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

void
SCF::compute()
{
  local_ = (dynamic_cast<LocalSCMatrixKit*>(basis()->matrixkit().pointer()) ||
            dynamic_cast<ReplSCMatrixKit*>(basis()->matrixkit().pointer())) ? 1:0;
  
  const double hess_to_grad_acc = 1.0/100.0;
  if (hessian_needed())
    set_desired_gradient_accuracy(desired_hessian_accuracy()*hess_to_grad_acc);

  const double grad_to_val_acc = 1.0/100.0;
  if (gradient_needed())
    set_desired_value_accuracy(desired_gradient_accuracy()*grad_to_val_acc);

  double delta;
  if (value_needed()) {
    ExEnv::out0() << endl << indent
         << scprintf("SCF::compute: energy accuracy = %10.7e\n",
                     desired_value_accuracy())
         << endl;

    double eelec;
    delta = compute_vector(eelec);
      
    // this will be done elsewhere eventually
    double nucrep = molecule()->nuclear_repulsion_energy();
    double eother = 0.0;
    if (accumddh_.nonnull()) eother = accumddh_->e();
    ExEnv::out0() << endl << indent
         << scprintf("total scf energy = %15.10f", eelec+eother+nucrep)
         << endl;

    set_energy(eelec+eother+nucrep);
    set_actual_value_accuracy(delta);
  }
  else {
    delta = actual_value_accuracy();
  }

  if (gradient_needed()) {
    RefSCVector gradient = matrixkit()->vector(moldim());

    ExEnv::out0() << endl << indent
         << scprintf("SCF::compute: gradient accuracy = %10.7e\n",
                     desired_gradient_accuracy())
         << endl;

    compute_gradient(gradient);
    print_natom_3(gradient,"Total Gradient:");
    set_gradient(gradient);

    set_actual_gradient_accuracy(delta/grad_to_val_acc);
  }
  
  if (hessian_needed()) {
    RefSymmSCMatrix hessian = matrixkit()->symmmatrix(moldim());
    
    ExEnv::out0() << endl << indent
         << scprintf("SCF::compute: hessian accuracy = %10.7e\n",
                     desired_hessian_accuracy())
         << endl;

    compute_hessian(hessian);
    set_hessian(hessian);

    set_actual_hessian_accuracy(delta/grad_to_val_acc/hess_to_grad_acc);
  }
}

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

signed char *
SCF::init_pmax(double *pmat_data)
{
  double l2inv = 1.0/log(2.0);
  double tol = pow(2.0,-126.0);
  
  GaussianBasisSet& gbs = *basis().pointer();
  
  signed char * pmax = new signed char[i_offset(gbs.nshell())];

  int ish, jsh, ij;
  for (ish=ij=0; ish < gbs.nshell(); ish++) {
    int istart = gbs.shell_to_function(ish);
    int iend = istart + gbs(ish).nfunction();
    
    for (jsh=0; jsh <= ish; jsh++,ij++) {
      int jstart = gbs.shell_to_function(jsh);
      int jend = jstart + gbs(jsh).nfunction();
      
      double maxp=0, tmp;

      for (int i=istart; i < iend; i++) {
        int ijoff = i_offset(i) + jstart;
        for (int j=jstart; j < ((ish==jsh) ? i+1 : jend); j++,ijoff++)
          if ((tmp=fabs(pmat_data[ijoff])) > maxp)
            maxp=tmp;
      }

      if (maxp <= tol)
        maxp=tol;

      long power = long(ceil(log(maxp)*l2inv));
      if (power < SCHAR_MIN) pmax[ij] = SCHAR_MIN;
      else if (power > SCHAR_MAX) pmax[ij] = SCHAR_MAX;
      else pmax[ij] = (signed char) power;
    }
  }

  return pmax;
}

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

RefSymmSCMatrix
SCF::get_local_data(const RefSymmSCMatrix& m, double*& p, Access access)
{
  RefSymmSCMatrix l = m;
  
  if (!dynamic_cast<LocalSymmSCMatrix*>(l.pointer())
      && !dynamic_cast<ReplSymmSCMatrix*>(l.pointer())) {
    Ref<SCMatrixKit> k = new ReplSCMatrixKit;
    l = k->symmmatrix(m.dim());
    l->convert(m);

    if (access == Accum)