symmcoor.cc

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

//
// symmcoor.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.
//

#include <math.h>

#include <util/misc/formio.h>
#include <util/state/stateio.h>
#include <math/scmat/matrix.h>
#include <math/scmat/elemop.h>
#include <chemistry/molecule/localdef.h>
#include <chemistry/molecule/molecule.h>
#include <chemistry/molecule/coor.h>
#include <chemistry/molecule/simple.h>

#include <util/container/bitarray.h>

using namespace std;
using namespace sc;

#define VERBOSE 0

namespace sc {

///////////////////////////////////////////////////////////////////////////
// SymmCoorTransform

class SymmCoorTransform: public NonlinearTransform {
  private:
    Ref<Molecule> molecule_;
    RefSCDimension dnatom3_;
    Ref<SetIntCoor> oldintcoor_;
    Ref<SetIntCoor> newintcoor_;
    Ref<SCMatrixKit> matrixkit_;
    int transform_hessian_;
  public:
    SymmCoorTransform(const Ref<Molecule>& molecule,
                      const RefSCDimension& dnatom3,
                      const Ref<SCMatrixKit>& kit,
                      const Ref<SetIntCoor>& oldintcoor,
                      const Ref<SetIntCoor>& newintcoor,
                      int transform_hessian);
    void to_cartesian(const RefSCVector& new_internal);
    void transform_coordinates(const RefSCVector& x);
    void transform_hessian(const RefSymmSCMatrix& h);
};

SymmCoorTransform::SymmCoorTransform(const Ref<Molecule>& molecule,
                                     const RefSCDimension& dnatom3,
                                     const Ref<SCMatrixKit>& kit,
                                     const Ref<SetIntCoor>& oldintcoor,
                                     const Ref<SetIntCoor>& newintcoor,
                                     int transform_hessian)
{
  molecule_ = new Molecule(*molecule.pointer());
  dnatom3_ = dnatom3;
  matrixkit_ = kit;
  oldintcoor_ = oldintcoor;
  newintcoor_ = newintcoor;
  transform_hessian_ = transform_hessian;
}

void
SymmCoorTransform::to_cartesian(const RefSCVector& new_internal)
{
  Ref<SCMatrixKit> kit = new_internal.kit();

  // get a reference to Molecule for convenience
  Molecule& molecule = *(molecule_.pointer());

  RefSCDimension dim = new_internal.dim();

  // don't bother updating the bmatrix when the error is less than this
  const double update_tolerance = 1.0e-3;
  const double cartesian_tolerance = 1.0e-8;

  // compute the internal coordinate displacements
  RefSCVector old_internal(new_internal.dim(),kit);

  RefSCMatrix internal_to_cart_disp;
  double maxabs_cart_diff = 0.0;
  for (int step = 0; step < 100; step++) {
      // compute the old internal coordinates
      oldintcoor_->update_values(molecule_);
      oldintcoor_->values_to_vector(old_internal);

      // the displacements
      RefSCVector displacement = new_internal - old_internal;

      if (maxabs_cart_diff>update_tolerance
          || internal_to_cart_disp.null()) {

          int i;
          RefSCMatrix bmat(dim,dnatom3_,kit);

          // form the bmatrix
          oldintcoor_->bmat(molecule_,bmat);

          // Compute the singular value decomposition of B
          RefSCMatrix U(dim,dim,kit);
          RefSCMatrix V(dnatom3_,dnatom3_,kit);
          RefSCDimension min;
          if (dnatom3_.n()<dim.n()) min = dnatom3_;
          else min = dim;
          int nmin = min.n();
          RefDiagSCMatrix sigma(min,kit);
          bmat.svd(U,sigma,V);

          // compute the epsilon rank of B
          int rank = 0;
          for (i=0; i<nmin; i++) {
              if (fabs(sigma(i)) > 0.0001) rank++;
            }

          RefSCDimension drank = new SCDimension(rank);
          RefDiagSCMatrix sigma_i(drank,kit);
          for (i=0; i<rank; i++) {
              sigma_i(i) = 1.0/sigma(i);
            }
          RefSCMatrix Ur(dim, drank, kit);
          RefSCMatrix Vr(dnatom3_, drank, kit);
          Ur.assign_subblock(U,0, dim.n()-1, 0, drank.n()-1, 0, 0);
          Vr.assign_subblock(V,0, dnatom3_.n()-1, 0, drank.n()-1, 0, 0);
          internal_to_cart_disp = Vr * sigma_i * Ur.t();
        }

      // compute the cartesian displacements
      RefSCVector cartesian_displacement = internal_to_cart_disp*displacement;
      // update the geometry
      for (int i=0; i < dnatom3_.n(); i++) {
          molecule.r(i/3,i%3) += cartesian_displacement(i);
        }

      // check for convergence
      Ref<SCElementMaxAbs> maxabs = new SCElementMaxAbs();
      Ref<SCElementOp> op = maxabs.pointer();
      cartesian_displacement.element_op(op);
      maxabs_cart_diff = maxabs->result();
      if (maxabs_cart_diff < cartesian_tolerance) {
          oldintcoor_->update_values(molecule_);
          return;
        }
    }

  ExEnv::err0() << indent
       << "WARNING: SymmCoorTransform::to_cartesian(RefSCVector&):"
       << " too many iterations in geometry update\n";

  new_internal.print("SymmCoorTransform: desired internal coordinates");
  (new_internal
   - old_internal).print("SymmCoorTransform: difference of desired and actual coordinates");

  abort();
}

void
SymmCoorTransform::transform_coordinates(const RefSCVector& x)
{
  if (x.null()) return;

  Ref<SCMatrixKit> kit = x.kit();
  RefSCDimension dim = x.dim();

  // using the old coordinates update molecule
  to_cartesian(x);

  // compute the new coordinates
  newintcoor_->update_values(molecule_);
  newintcoor_->values_to_vector(x);

  // compute the linear transformation information

  // the old B matrix
  RefSCMatrix B(dim, dnatom3_, kit);
  oldintcoor_->bmat(molecule_, B);

  // get the B matrix for the new coordinates
  RefSCMatrix Bnew(dim, dnatom3_, kit);
  newintcoor_->update_values(molecule_);
  newintcoor_->bmat(molecule_, Bnew);

  // the transform from cartesian to new internal coordinates
  RefSymmSCMatrix bmbt(dim,kit);
  bmbt.assign(0.0);
  bmbt.accumulate_symmetric_product(Bnew);
  RefSCMatrix cart_to_new_internal = bmbt.gi() * Bnew;
  Bnew = 0;
  bmbt = 0;

  linear_transform_ = cart_to_new_internal * B.t();
#if VERBOSE
  linear_transform_.print("old internal to new");
#endif
}

void
SymmCoorTransform::transform_hessian(const RefSymmSCMatrix& h)
{
  if (transform_hessian_) {
      NonlinearTransform::transform_hessian(h);
    }
  else {
      ExEnv::err0() << indent
           << "WARNING: SymmCoorTransform::transform_hessian: "
           << "skipping hessian transform";
    }
}

///////////////////////////////////////////////////////////////////////////
// members of SymmMolecularCoor

static ClassDesc SymmMolecularCoor_cd(
  typeid(SymmMolecularCoor),"SymmMolecularCoor",1,"public IntMolecularCoor",
  0, create<SymmMolecularCoor>, create<SymmMolecularCoor>);

SymmMolecularCoor::SymmMolecularCoor(Ref<Molecule>&mol):
  IntMolecularCoor(mol)
{
  init();
}

SymmMolecularCoor::SymmMolecularCoor(const Ref<KeyVal>& keyval):
  IntMolecularCoor(keyval)
{
  init();

  int itmp;
  double dtmp;
  itmp = keyval->booleanvalue("change_coordinates");
  if (keyval->error() == KeyVal::OK) change_coordinates_ = itmp;
  itmp = keyval->booleanvalue("transform_hessian");
  if (keyval->error() == KeyVal::OK) transform_hessian_ = itmp;
  dtmp = keyval->doublevalue("max_kappa2");
  if (keyval->error() == KeyVal::OK) max_kappa2_ = dtmp;
}

SymmMolecularCoor::SymmMolecularCoor(StateIn& s):
  IntMolecularCoor(s)
{
  s.get(change_coordinates_);
  s.get(transform_hessian_);
  s.get(max_kappa2_);
}

SymmMolecularCoor::~SymmMolecularCoor()
{
}

void
SymmMolecularCoor::save_data_state(StateOut&s)