ossscf.cc

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

//
// ossscf.cc --- implementation of the open shell singlet SCF 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 <util/misc/timer.h>
#include <util/misc/formio.h>
#include <util/state/stateio.h>

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

#include <math/optimize/scextrapmat.h>

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

#include <chemistry/qc/scf/scflocal.h>
#include <chemistry/qc/scf/scfops.h>
#include <chemistry/qc/scf/effh.h>
#include <chemistry/qc/scf/ossscf.h>

using namespace std;
using namespace sc;

///////////////////////////////////////////////////////////////////////////
// OSSSCF

static ClassDesc OSSSCF_cd(
  typeid(OSSSCF),"OSSSCF",1,"public SCF",
  0, 0, 0);

OSSSCF::OSSSCF(StateIn& s) :
  SavableState(s),
  SCF(s),
  cl_fock_(this),
  op_focka_(this),
  op_fockb_(this)
{
  cl_fock_.result_noupdate() =
    basis_matrixkit()->symmmatrix(so_dimension());
  cl_fock_.restore_state(s);
  cl_fock_.result_noupdate().restore(s);
  
  op_focka_.result_noupdate() =
    basis_matrixkit()->symmmatrix(so_dimension());
  op_focka_.restore_state(s);
  op_focka_.result_noupdate().restore(s);
  
  op_fockb_.result_noupdate() =
    basis_matrixkit()->symmmatrix(so_dimension());
  op_fockb_.restore_state(s);
  op_fockb_.result_noupdate().restore(s);
  
  s.get(user_occupations_);
  s.get(tndocc_);
  s.get(nirrep_);
  s.get(ndocc_);
  s.get(osa_);
  s.get(osb_);

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

OSSSCF::OSSSCF(const Ref<KeyVal>& keyval) :
  SCF(keyval),
  cl_fock_(this),
  op_focka_(this),
  op_fockb_(this)
{
  cl_fock_.compute()=0;
  cl_fock_.computed()=0;
  
  op_focka_.compute()=0;
  op_focka_.computed()=0;
  
  op_fockb_.compute()=0;
  op_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);

  // figure out how many doubly occupied shells there are
  if (keyval->exists("ndocc")) {
    tndocc_ = keyval->intvalue("ndocc");
  } else {
    tndocc_ = (nelectrons-2)/2;
    if ((nelectrons-2)%2) {
      ExEnv::err0() << endl << indent
           << "OSSSCF::init: Warning, there's a leftover electron.\n"
           << incindent
           << indent << "total_charge = " << charge << endl
           << indent << "total nuclear charge = " << Znuc << endl
           << indent << "ndocc_ = " << tndocc_ << endl << decindent;
    }
  }

  ExEnv::out0() << endl << indent << "OSSSCF::init: total charge = "
       << Znuc-2*tndocc_-2 << endl << endl;

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

  if (nirrep_==1) {
    ExEnv::err0() << indent << "OSSSCF::init: cannot do C1 symmetry\n";
    abort();
  }

  osa_=-1;
  osb_=-1;
  
  ndocc_ = read_occ(keyval, "docc", nirrep_);
  int *nsocc = read_occ(keyval, "socc", nirrep_);
  if (ndocc_ && nsocc) {
    user_occupations_=1;
    for (int i=0; i < nirrep_; i++) {
      int nsi = nsocc[i];
      if (nsi && osa_<0)
        osa_=i;
      else if (nsi && osb_<0)
        osb_=i;
      else if (nsi) {
        ExEnv::err0() << indent << "OSSSCF::init: too many open shells\n";
        abort();
      }
    }
    delete[] nsocc;
  }
  else if (ndocc_ && !nsocc || !ndocc_ && nsocc) {
    ExEnv::outn() << "ERROR: OSSSCF: only one of docc and socc specified: "
                 << "give both or none" << endl;
    abort();
  }
  else {
    ndocc_=0;
    user_occupations_=0;
    set_occupations(0);
  }

  int i;
  ExEnv::out0() << indent << "docc = [";
  for (i=0; i < nirrep_; i++)
    ExEnv::out0() << " " << ndocc_[i];
  ExEnv::out0() << " ]\n";

  ExEnv::out0() << indent << "socc = [";
  for (i=0; i < nirrep_; i++)
    ExEnv::out0() << " " << (i==osa_ || i==osb_) ? 1 : 0;
  ExEnv::out0() << " ]\n";

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

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

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

OSSSCF::~OSSSCF()
{
  if (ndocc_) {
    delete[] ndocc_;
    ndocc_=0;
  }
}

void
OSSSCF::save_data_state(StateOut& s)
{
  SCF::save_data_state(s);
  
  cl_fock_.save_data_state(s);
  cl_fock_.result_noupdate().save(s);
  
  op_focka_.save_data_state(s);
  op_focka_.result_noupdate().save(s);
  
  op_fockb_.save_data_state(s);
  op_fockb_.result_noupdate().save(s);
  
  s.put(user_occupations_);
  s.put(tndocc_);
  s.put(nirrep_);
  s.put(ndocc_,nirrep_);
  s.put(osa_);
  s.put(osb_);
}

double
OSSSCF::occupation(int ir, int i)
{
  if (i < ndocc_[ir])
    return 2.0;
  else if ((ir==osa_ || ir==osb_) && (i == ndocc_[ir]))
    return 1.0;
  return 0.0;
}

double
OSSSCF::alpha_occupation(int ir, int i)
{
  if (i < ndocc_[ir] || (ir==osa_ && i==ndocc_[ir]))
    return 1.0;
  return 0.0;
}

double
OSSSCF::beta_occupation(int ir, int i)
{
  if (i < ndocc_[ir] || (ir==osb_ && i==ndocc_[ir]))
    return 1.0;
  return 0.0;
}

int
OSSSCF::n_fock_matrices() const
{
  return 3;
}

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

  if (n==0)
    return cl_fock_.result();
  else if (n==1)
    return op_focka_.result();
  else
    return op_fockb_.result();
}

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

void
OSSSCF::print(ostream&o) const
{
  int i;
  
  SCF::print(o);
  o << indent << "OSSSCF Parameters:\n" << incindent
    << indent << "ndocc = " << tndocc_ << endl
    << indent << "docc = [";
  for (i=0; i < nirrep_; i++)
    o << " " << ndocc_[i];
  o << " ]" << endl;

  o << indent << "socc = [";
  for (i=0; i < nirrep_; i++)
    o << " " << (i==osa_ || i==osb_) ? 1 : 0;
  o << " ]" << endl << decindent << endl;
}

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

void
OSSSCF::set_occupations(const RefDiagSCMatrix& ev)
{
  if (user_occupations_)
    return;
  
  int i,j;
  
  RefDiagSCMatrix evals;
  
  if (ev.null()) {
    initial_vector(0);
    evals = eigenvalues_.result_noupdate();
  }
  else
    evals = ev;

  // first convert evals to something we can deal with easily
  BlockedDiagSCMatrix *evalsb = require_dynamic_cast<BlockedDiagSCMatrix*>(evals,
                                                 "OSSSCF::set_occupations");
  
  double **vals = new double*[nirrep_];
  for (i=0; i < nirrep_; i++) {
    int nf=oso_dimension()->blocks()->size(i);
    if (nf) {
      vals[i] = new double[nf];
      evalsb->block(i)->convert(vals[i]);
    } else {
      vals[i] = 0;
    }
  }

  // now loop to find the tndocc_ lowest eigenvalues and populate those
  // MO's
  int *newdocc = new int[nirrep_];
  memset(newdocc,0,sizeof(int)*nirrep_);

  for (i=0; i < tndocc_; i++) {
    // find lowest eigenvalue
    int lir=0,ln=0;
    double lowest=999999999;

    for (int ir=0; ir < nirrep_; ir++) {
      int nf=oso_dimension()->blocks()->size(ir);
      if (!nf)
        continue;
      for (j=0; j < nf; j++) {
        if (vals[ir][j] < lowest) {
          lowest=vals[ir][j];
          lir=ir;
          ln=j;
        }
      }
    }
    vals[lir][ln]=999999999;
    newdocc[lir]++;
  }

  int osa=-1, osb=-1;
  
  for (i=0; i < 2; i++) {
    // find lowest eigenvalue
    int lir=0,ln=0;
    double lowest=999999999;

    for (int ir=0; ir < nirrep_; ir++) {
      int nf=oso_dimension()->blocks()->size(ir);
      if (!nf)
        continue;
      for (j=0; j < nf; j++) {
        if (vals[ir][j] < lowest) {
          lowest=vals[ir][j];
          lir=ir;
          ln=j;
        }
      }
    }
    vals[lir][ln]=999999999;

    if (!i) {
      osa=lir;
    } else {
      if (lir==osa) {
        i--;
        continue;
      }
      osb=lir;
    }
  }

  // get rid of vals
  for (i=0; i < nirrep_; i++)
    if (vals[i])
      delete[] vals[i];
  delete[] vals;

  if (!ndocc_) {
    ndocc_=newdocc;
    osa_=osa;
    osb_=osb;
  } else {
    // test to see if newocc is different from ndocc_
    for (i=0; i < nirrep_; i++) {
      if (ndocc_[i] != newdocc[i]) {
        ExEnv::err0() << indent << "OSSSCF::set_occupations:  WARNING!!!!\n"
             << incindent << indent
             << scprintf("occupations for irrep %d have changed\n", i+1)
             << indent
             << scprintf("ndocc was %d, changed to %d", ndocc_[i], newdocc[i])
             << endl << decindent;
      }
      if ((osa != osa_ && osa != osb_) || (osb != osb_ && osb != osa_)) {
        ExEnv::err0() << indent << "OSSSCF::set_occupations:  WARNING!!!!\n"
             << incindent << indent
             << "open shell occupations have changed"
             << endl << decindent;
        osa_=osa;
        osb_=osb;