solvent.cc

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

//
// solvent.cc
//
// Copyright (C) 1997 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/misc/timer.h>
#include <util/misc/formio.h>
#include <util/state/stateio.h>
#include <chemistry/qc/basis/petite.h>
#include <chemistry/qc/wfn/solvent.h>

#include <math/isosurf/volume.h>
#include <chemistry/qc/dft/integrator.h>
#include <chemistry/qc/dft/functional.h>

#include <iomanip>

using namespace std;
using namespace sc;

namespace sc {

//. The \clsnm{NElFunctional} computes the number of electrons.
//. It is primarily for testing the integrator.
class NElInShapeFunctional: public DenFunctional {
  private:
    Ref<Volume> vol_;
    double isoval_;
  public:
    NElInShapeFunctional(const Ref<Volume> &, double);
    ~NElInShapeFunctional();

    void point(const PointInputData&, PointOutputData&);
};

/////////////////////////////////////////////////////////////////////////////
// NElFunctional

static ClassDesc NElInShapeFunctional_cd(
  typeid(NElInShapeFunctional),"NElInShapeFunctional",1,"public DenFunctional",
  0, 0, 0);

NElInShapeFunctional::NElInShapeFunctional(const Ref<Volume>& vol,
                                           double isoval)
{
  vol_ = vol;
  isoval_ = isoval;
}

NElInShapeFunctional::~NElInShapeFunctional()
{
}

void
NElInShapeFunctional::point(const PointInputData &id,
                            PointOutputData &od)
{
  vol_->set_x(id.r);
  if (vol_->value() <= isoval_) {
      od.energy = id.a.rho + id.b.rho;
    }
  else {
      od.energy = 0.0;
    }
}

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

static ClassDesc BEMSolventH_cd(
  typeid(BEMSolventH),"BEMSolventH",1,"public AccumH",
  0, create<BEMSolventH>, create<BEMSolventH>);

BEMSolventH::BEMSolventH(const Ref<KeyVal>&keyval):
  AccumH(keyval)
{
  charge_positions_ = 0;
  normals_ = 0;
  efield_dot_normals_ = 0;
  charges_ = 0;
  charges_n_ = 0;
  solvent_ << keyval->describedclassvalue("solvent");
  gamma_ = keyval->doublevalue("gamma");
  if (keyval->error() != KeyVal::OK) {
      Ref<Units> npm = new Units("dyne/cm");
      gamma_ = 72.75 * npm->to_atomic_units();
    }
  // If onebody add a term to the one body hamiltonian, h.
  // Otherwise the energy contribution is scalar.
  onebody_ = keyval->booleanvalue("onebody");
  if (keyval->error() != KeyVal::OK) onebody_ = 1;
  // Normalize the charges if normalize_q is set.
  normalize_q_ = keyval->booleanvalue("normalize_q");
  if (keyval->error() != KeyVal::OK) normalize_q_ = 1;
  // Compute separately contributes to the energy from surfaces
  // charges induced by the nuclear and electronic charge densities.
  separate_surf_charges_ = keyval->booleanvalue("separate_surf_charges");
  if (keyval->error() != KeyVal::OK) separate_surf_charges_ = 0;
  // The Cammi-Tomasi Y term is set equal to the J term (as it formally is).
  y_equals_j_ = keyval->booleanvalue("y_equals_j");
  if (keyval->error() != KeyVal::OK) y_equals_j_ = 0;
  // As a test, integrate the number of electrons inside the surface.
  integrate_nelectron_ = keyval->booleanvalue("integrate_nelectron");
  if (keyval->error() != KeyVal::OK) integrate_nelectron_ = 0;
}

BEMSolventH::BEMSolventH(StateIn&s):
  SavableState(s),
  AccumH(s)
{
  charge_positions_ = 0;
  normals_ = 0;
  efield_dot_normals_ = 0;
  charges_ = 0;
  charges_n_ = 0;
  escalar_ = 0;

  wfn_ << SavableState::restore_state(s);
  //solvent_.restore_state(s);
  abort();
}

BEMSolventH::~BEMSolventH()
{
  // just in case
  done();
}

void
BEMSolventH::save_data_state(StateOut&s)
{
  AccumH::save_data_state(s);

  SavableState::save_state(wfn_.pointer(),s);
  //solvent_.save_state(s);
  abort();
}

void
BEMSolventH::init(const Ref<Wavefunction>& wfn)
{
  tim_enter("solvent");
  tim_enter("init");
  wfn_ = wfn;
  // just in case
  done();
  solvent_->init();
  charge_positions_ = solvent_->alloc_charge_positions();
  normals_ = solvent_->alloc_normals();
  efield_dot_normals_ = solvent_->alloc_efield_dot_normals();
  charges_ = solvent_->alloc_charges();
  charges_n_ = solvent_->alloc_charges();

  // get the positions of the charges
  solvent_->charge_positions(charge_positions_);

  // get the surface normals
  solvent_->normals(normals_);

  if (integrate_nelectron_) {
      Ref<DenIntegrator> integrator = new RadialAngularIntegrator();
      Ref<DenFunctional> functional
          = new NElInShapeFunctional(solvent_->surface()->volume_object(),
                                     solvent_->surface()->isovalue());
      integrator->init(wfn_);
      integrator->integrate(functional);
      integrator->done();
      ExEnv::out0() << indent
           << scprintf("N(e) in isosurf = %12.8f", integrator->value())
           << endl;
    }

  edisprep_ = solvent_->disprep();

  tim_exit("init");
  tim_exit("solvent");
}

// This adds J + X to h, where J and X are the matrices defined
// by Canni and Tomasi, J Comp Chem, 16(12), 1457, 1995.
// The resulting SCF free energy expression is
//    G = 1/2TrP[h' + F'] + Une + Unn + Vnn
//        -1/2(Uee+Uen+Une+Unn)
// which in the Canni-Tomasi notation is
//      = 1/2TrP[h+1/2(X+J+Y+G)] + Vnn + 1/2Unn
// which is identical to the Canni-Tomasi energy expression.
// My Fock matrix is
//       F' = h + J + X + G
// while the Canni-Tomasi Fock matrix is F' = h + 1/2(J+Y) + X + G.
// However, since J = Y formally, (assuming no numerical errors
// and all charge is enclosed, Canni-Tomasi use F' = h + J + X + G
// to get better numerical results.
//
//   If the y_equals_j option is true, the energy expression used
// here is G = 1/2TrP[h+1/2(X+2J+G)] + Vnn + 1/2Unn, however, THIS
// IS NOT RECOMMENDED.
void
BEMSolventH::accum(const RefSymmSCMatrix& h)
{
  tim_enter("solvent");
  tim_enter("accum");
  int i,j;

  //// compute the polarization charges

  // compute the e-field at each point and dot with normals
  tim_enter("efield");
  int ncharge = solvent_->ncharge();
  Ref<EfieldDotVectorData> efdn_dat = new EfieldDotVectorData;
  Ref<OneBodyInt> efdn = wfn_->integral()->efield_dot_vector(efdn_dat);
  Ref<SCElementOp> efdn_op = new OneBodyIntOp(efdn);
  RefSymmSCMatrix ao_density = wfn_->ao_density()->copy();
  RefSymmSCMatrix efdn_mat(ao_density->dim(), ao_density->kit());
  // for the scalar products, scale the density's off-diagonals by two
  ao_density->scale(2.0);
  ao_density->scale_diagonal(0.5);
  Ref<SCElementScalarProduct> sp = new SCElementScalarProduct;
  Ref<SCElementOp2> generic_sp(sp.pointer());
  for (i=0; i<ncharge; i++) {
      efdn_dat->set_position(charge_positions_[i]);
      efdn_dat->set_vector(normals_[i]);
      efdn->reinitialize();