obwfn.cc

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

//
// obwfn.cc
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <[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 <util/state/stateio.h>
#include <util/misc/formio.h>

#include <math/symmetry/corrtab.h>
#include <math/scmat/local.h>
#include <math/scmat/blocked.h>

#include <chemistry/qc/basis/integral.h>
#include <chemistry/qc/basis/obint.h>
#include <chemistry/qc/basis/petite.h>
#include <chemistry/qc/wfn/obwfn.h>

using namespace std;
using namespace sc;

#define DEBUG 0

#ifndef DBL_EPSILON
#define DBL_EPSILON 1.0e-15
#endif

static ClassDesc OneBodyWavefunction_cd(
  typeid(OneBodyWavefunction),"OneBodyWavefunction",1,"public Wavefunction",
  0, 0, 0);

OneBodyWavefunction::OneBodyWavefunction(const Ref<KeyVal>&keyval):
  Wavefunction(keyval),
  density_(this),
  oso_eigenvectors_(this),
  eigenvalues_(this),
  nirrep_(0),
  nvecperirrep_(0),
  occupations_(0),
  alpha_occupations_(0),
  beta_occupations_(0)
{
  double acc = keyval->doublevalue("eigenvector_accuracy");
  if (keyval->error() != KeyVal::OK)
    acc = value_.desired_accuracy();
  
  oso_eigenvectors_.set_desired_accuracy(acc);
  eigenvalues_.set_desired_accuracy(acc);

  if (oso_eigenvectors_.desired_accuracy() < DBL_EPSILON) {
    oso_eigenvectors_.set_desired_accuracy(DBL_EPSILON);
    eigenvalues_.set_desired_accuracy(DBL_EPSILON);
  }
}

OneBodyWavefunction::OneBodyWavefunction(StateIn&s):
  SavableState(s),
  Wavefunction(s),
  density_(this),
  oso_eigenvectors_(this),
  eigenvalues_(this),
  nirrep_(0),
  nvecperirrep_(0),
  occupations_(0),
  alpha_occupations_(0),
  beta_occupations_(0)
{
  oso_eigenvectors_.result_noupdate() =
    basis_matrixkit()->matrix(oso_dimension(), oso_dimension());
  oso_eigenvectors_.restore_state(s);
  oso_eigenvectors_.result_noupdate().restore(s);

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

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

OneBodyWavefunction::~OneBodyWavefunction()
{
  if (nvecperirrep_) {
    delete[] nvecperirrep_;
    delete[] occupations_;
    delete[] alpha_occupations_;
    delete[] beta_occupations_;
  }
  nirrep_=0;
  nvecperirrep_=0;
  occupations_=0;
  alpha_occupations_=0;
  beta_occupations_=0;
}

void
OneBodyWavefunction::save_data_state(StateOut&s)
{
  Wavefunction::save_data_state(s);

  oso_eigenvectors_.save_data_state(s);
  oso_eigenvectors_.result_noupdate().save(s);

  eigenvalues_.save_data_state(s);
  eigenvalues_.result_noupdate().save(s);

  density_.save_data_state(s);
  density_.result_noupdate().save(s);
}

RefSCMatrix
OneBodyWavefunction::projected_eigenvectors(const Ref<OneBodyWavefunction>& owfn,
                                            int alp)
{
  //............................................................
  // first obtain the guess density matrix

  // The old density in the old SO basis
  RefSymmSCMatrix oldP_so;
  if (owfn->spin_unrestricted()) {
    if (alp) oldP_so = owfn->alpha_density();
    else oldP_so = owfn->beta_density();
  }
  else oldP_so = owfn->density();

  ExEnv::out0() << endl << indent
       << "Projecting the guess density.\n"
       << endl;
  ExEnv::out0() << incindent;

  // The old overlap
  RefSymmSCMatrix oldS = owfn->overlap();
  ExEnv::out0() << indent
       << "The number of electrons in the guess density = "
       << (oldP_so*oldS).trace() << endl;

  // Transform the old SO overlap into the orthogonal SO basis, oSO
  RefSCMatrix old_so_to_oso = owfn->so_to_orthog_so();
  RefSymmSCMatrix oldP_oso(owfn->oso_dimension(), owfn->basis_matrixkit());
  oldP_oso->assign(0.0);
  oldP_oso->accumulate_transform(old_so_to_oso, oldP_so);

  //............................................................
  // transform the guess density into the current basis

  // the transformation matrix is the new basis/old basis overlap
  integral()->set_basis(owfn->basis(), basis());
  RefSCMatrix old_to_new_ao(owfn->basis()->basisdim(), basis()->basisdim(),
                            basis()->matrixkit());
  Ref<SCElementOp> op = new OneBodyIntOp(integral()->overlap());
  old_to_new_ao.assign(0.0);
  old_to_new_ao.element_op(op);
  op = 0;
  integral()->set_basis(basis());

  // now must transform the transform into the SO basis
  Ref<PetiteList> pl = integral()->petite_list();
  Ref<PetiteList> oldpl = owfn->integral()->petite_list();
  RefSCMatrix blocked_old_to_new_ao(oldpl->AO_basisdim(), pl->AO_basisdim(),
                                    basis()->so_matrixkit());
  blocked_old_to_new_ao->convert(old_to_new_ao);
  RefSCMatrix old_to_new_so
    = oldpl->sotoao() * blocked_old_to_new_ao * pl->aotoso();

  // now must transform the transform into the orthogonal SO basis
  RefSCMatrix so_to_oso = so_to_orthog_so();
  RefSCMatrix old_to_new_oso = owfn->so_to_orthog_so_inverse().t()
                             * old_to_new_so
                             * so_to_oso.t();
  old_so_to_oso = 0;
  old_to_new_so = 0;
  
  // The old density transformed to the new orthogonal SO basis
  RefSymmSCMatrix newP_oso(oso_dimension(), basis_matrixkit());
  newP_oso->assign(0.0);
  newP_oso->accumulate_transform(old_to_new_oso.t(), oldP_oso);
  old_to_new_oso = 0;
  oldP_oso = 0;
  //newP_oso.print("projected orthoSO density");

  ExEnv::out0() << indent
       << "The number of electrons in the projected density = "
       << newP_oso.trace() << endl;

  //............................................................

  // reverse the sign of the density so the eigenvectors will
  // be ordered in the right way
  newP_oso.scale(-1.0);

  // use the guess density in the new basis to find the orbitals
  // (the density should be diagonal in the MO basis--this proceedure
  // will not give canonical orbitals, but they should at least give
  // a decent density)
  RefDiagSCMatrix newP_oso_vals(newP_oso.dim(), basis_matrixkit());
  RefSCMatrix newP_oso_vecs(newP_oso.dim(), newP_oso.dim(), basis_matrixkit());
  newP_oso.diagonalize(newP_oso_vals, newP_oso_vecs);
  //newP_oso_vals.print("eigenvalues of projected density");

  // Reordering of the vectors isn't needed because of the way
  // the density was scaled above.
  RefSCMatrix newvec_oso = newP_oso_vecs;

  if (debug_ >= 2) {
      newvec_oso.print("projected ortho SO vector");
      so_to_oso.print("SO to ortho SO transformation");
    }

  ExEnv::out0() << decindent;
  return newvec_oso;
}

// this is a hack for big basis sets where the core hamiltonian eigenvalues
// are total garbage.  Use the old wavefunction's occupied eigenvalues, and
// set all others to 99.

RefDiagSCMatrix
OneBodyWavefunction::projected_eigenvalues(const Ref<OneBodyWavefunction>& owfn,
                                           int alp)
{
  // get the old eigenvalues and the new core hamiltonian evals
  RefDiagSCMatrix oval;
  if (owfn->spin_unrestricted()) {
    if (alp)
      oval = owfn->alpha_eigenvalues();
    else
      oval = owfn->beta_eigenvalues();
  } else
    oval = owfn->eigenvalues();

  BlockedDiagSCMatrix *ovalp =
    require_dynamic_cast<BlockedDiagSCMatrix*>(oval.pointer(),
      "OneBodyWavefunction::projected_eigenvalues: oval"
      );
    
  // get the core hamiltonian eigenvalues
  RefDiagSCMatrix val;
  hcore_guess(val);
  BlockedDiagSCMatrix *valp =
    require_dynamic_cast<BlockedDiagSCMatrix*>(val.pointer(),
      "OneBodyWavefunction::projected_eigenvalues: val"
      );

  RefSCDimension oso = oso_dimension();
  RefSCDimension ooso = owfn->oso_dimension();

  for (int irrep=0; irrep < valp->nblocks(); irrep++) {
    // find out how many occupied orbitals there should be

    int nf = oso->blocks()->size(irrep);
    int nfo = ooso->blocks()->size(irrep);

    int nocc = 0;
    if (owfn->spin_unrestricted()) {
      if (alp)
        while (owfn->alpha_occupation(irrep,nocc) &&
               nocc < nfo) nocc++;
      else
        while (owfn->beta_occupation(irrep,nocc) &&
               nocc < nfo) nocc++;
    } else
      while (owfn->occupation(irrep,nocc) &&
             nocc < nfo) nocc++;

    if (!nf)
      continue;
    
    double *vals = new double[nf];
    valp->block(irrep)->convert(vals);

    int i;
    if (nfo) {
      double *ovals = new double[nfo];
      ovalp->block(irrep)->convert(ovals);
      for (i=0; i < nocc; i++) vals[i] = ovals[i];
      delete[] ovals;
    }
    
    for (i=nocc; i < nf; i++)
      vals[i] = 99.0;

    valp->block(irrep)->assign(vals);

    delete[] vals;