imcoor.cc

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

//
// imcoor.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 DEFAULT_SIMPLE_TOLERANCE 1.0e-3

///////////////////////////////////////////////////////////////////////////
// members of IntMolecularCoor

static ClassDesc IntMolecularCoor_cd(
  typeid(IntMolecularCoor),"IntMolecularCoor",6,"public MolecularCoor",
  0, 0, 0);

IntMolecularCoor::IntMolecularCoor(Ref<Molecule>&mol):
  MolecularCoor(mol),
  update_bmat_(0),
  only_totally_symmetric_(1),
  symmetry_tolerance_(1.0e-5),
  simple_tolerance_(DEFAULT_SIMPLE_TOLERANCE),
  coordinate_tolerance_(1.0e-7),
  cartesian_tolerance_(1.0e-12),
  scale_bonds_(1.0),
  scale_bends_(1.0),
  scale_tors_(1.0),
  scale_outs_(1.0),
  given_fixed_values_(0),
  decouple_bonds_(0),
  decouple_bends_(0),
  max_update_steps_(100),
  max_update_disp_(0.5),
  form_print_simples_(0),
  form_print_variable_(0),
  form_print_constant_(0),
  form_print_molecule_(0)
{
  new_coords();
  generator_ = new IntCoorGen(mol);
}

IntMolecularCoor::IntMolecularCoor(const Ref<KeyVal>& keyval):
  MolecularCoor(keyval),
  update_bmat_(0),
  only_totally_symmetric_(1),
  symmetry_tolerance_(1.0e-5),
  simple_tolerance_(DEFAULT_SIMPLE_TOLERANCE),
  coordinate_tolerance_(1.0e-7),
  cartesian_tolerance_(1.0e-12),
  scale_bonds_(1.0),
  scale_bends_(1.0),
  scale_tors_(1.0),
  scale_outs_(1.0),
  decouple_bonds_(0),
  decouple_bends_(0)
{
  // intialize the coordinate sets
  new_coords();

  // actually read the keyval info
  read_keyval(keyval);
}

IntMolecularCoor::IntMolecularCoor(StateIn& s):
  MolecularCoor(s)
{
  generator_ << SavableState::restore_state(s);

  if (s.version(::class_desc<IntMolecularCoor>()) >= 3) {
      s.get(decouple_bonds_);
      s.get(decouple_bends_);
    }
  else {
      decouple_bonds_ = 0;
      decouple_bends_ = 0;
    }
  
  if (s.version(::class_desc<IntMolecularCoor>()) >= 2) {
    s.get(max_update_steps_);
    s.get(max_update_disp_);
    s.get(given_fixed_values_);
  } else {
    max_update_steps_ = 100;
    max_update_disp_ = 0.5;
    given_fixed_values_ = 0;
  }
  
  if (s.version(::class_desc<IntMolecularCoor>()) >= 4) {
    s.get(form_print_simples_);
    s.get(form_print_variable_);
    s.get(form_print_constant_);
  } else {
    form_print_simples_ = 0;
    form_print_variable_ = 0;
    form_print_constant_ = 0;
  }
  
  if (s.version(::class_desc<IntMolecularCoor>()) >= 5) {
    s.get(form_print_molecule_);
  } else {
    form_print_molecule_ = 0;
  }

  dim_ << SavableState::restore_state(s);
  dvc_ << SavableState::restore_state(s);

  all_ << SavableState::restore_state(s);

  variable_ << SavableState::restore_state(s);
  constant_ << SavableState::restore_state(s);

  fixed_ << SavableState::restore_state(s);
  followed_ << SavableState::restore_state(s);

  if (s.version(::class_desc<IntMolecularCoor>()) >= 6)
      watched_ << SavableState::restore_state(s);

  bonds_ << SavableState::restore_state(s);
  bends_ << SavableState::restore_state(s);
  tors_ << SavableState::restore_state(s);
  outs_ << SavableState::restore_state(s);
  extras_ << SavableState::restore_state(s);

  s.get(update_bmat_);
  s.get(only_totally_symmetric_);
  s.get(scale_bonds_);
  s.get(scale_bends_);
  s.get(scale_tors_);
  s.get(scale_outs_);
  s.get(simple_tolerance_);
  s.get(symmetry_tolerance_);
  s.get(coordinate_tolerance_);
  s.get(cartesian_tolerance_);
}

void
IntMolecularCoor::new_coords()
{
  // intialize the coordinate sets
  all_ = new SetIntCoor; // all redundant coors
  variable_ = new SetIntCoor; // internal coors to be varied
  constant_ = new SetIntCoor; // internal coors to be fixed
  bonds_ = new SetIntCoor;
  bends_ = new SetIntCoor;
  tors_ = new SetIntCoor;
  outs_ = new SetIntCoor;
  extras_ = new SetIntCoor;
  fixed_ = new SetIntCoor;
  followed_ = 0;
  watched_ = 0;
}

void
IntMolecularCoor::read_keyval(const Ref<KeyVal>& keyval)
{
  variable_ << keyval->describedclassvalue("variable");
  if (variable_.null()) variable_ = new SetIntCoor;
  fixed_ << keyval->describedclassvalue("fixed");
  if (fixed_.null()) fixed_ = new SetIntCoor;
  followed_ << keyval->describedclassvalue("followed");
  watched_ << keyval->describedclassvalue("watched");

  decouple_bonds_ = keyval->booleanvalue("decouple_bonds");
  decouple_bends_ = keyval->booleanvalue("decouple_bends");

  given_fixed_values_ = keyval->booleanvalue("have_fixed_values");

  max_update_steps_ = keyval->intvalue("max_update_steps");
  if (keyval->error() != KeyVal::OK) max_update_steps_ = 100;

  max_update_disp_ = keyval->doublevalue("max_update_disp");
  if (keyval->error() != KeyVal::OK) max_update_disp_ = 0.5;

  generator_ << keyval->describedclassvalue("generator");

  if (generator_.null()) {
      // the extra_bonds list is given as a vector of atom numbers
      // (atom numbering starts at 1)
      int nextra_bonds = keyval->count("extra_bonds");
      nextra_bonds /= 2;
      int *extra_bonds;
      if (nextra_bonds) {
          extra_bonds = new int[nextra_bonds*2];
          for (int i=0; i<nextra_bonds*2; i++) {
              extra_bonds[i] = keyval->intvalue("extra_bonds",i);
              if (keyval->error() != KeyVal::OK) {
                  ExEnv::err0() << indent << "IntMolecularCoor:: keyval CTOR: "
                       << "problem reading \"extra_bonds:" << i << "\"\n";
                  abort();
                }
            }
        }
      else {
          extra_bonds = 0;
        }
      generator_ = new IntCoorGen(molecule_, nextra_bonds, extra_bonds);
    }
          

  update_bmat_ = keyval->booleanvalue("update_bmat");

  only_totally_symmetric_ = keyval->booleanvalue("only_totally_symmetric");
  if (keyval->error() != KeyVal::OK) only_totally_symmetric_ = 1;

  double tmp;
  tmp = keyval->doublevalue("scale_bonds");
  if (keyval->error() == KeyVal::OK) scale_bonds_ = tmp;
  tmp = keyval->doublevalue("scale_bends");
  if (keyval->error() == KeyVal::OK) scale_bends_ = tmp;
  tmp = keyval->doublevalue("scale_tors");
  if (keyval->error() == KeyVal::OK) scale_tors_ = tmp;
  tmp = keyval->doublevalue("scale_outs");
  if (keyval->error() == KeyVal::OK) scale_outs_ = tmp;
  tmp = keyval->doublevalue("symmetry_tolerance");
  if (keyval->error() == KeyVal::OK) symmetry_tolerance_ = tmp;
  tmp = keyval->doublevalue("simple_tolerance");
  if (keyval->error() == KeyVal::OK) simple_tolerance_ = tmp;
  tmp = keyval->doublevalue("coordinate_tolerance");
  if (keyval->error() == KeyVal::OK) coordinate_tolerance_ = tmp;
  tmp = keyval->doublevalue("cartesian_tolerance");
  if (keyval->error() == KeyVal::OK) cartesian_tolerance_ = tmp;

  form_print_simples_ = keyval->booleanvalue("form:print_simple");
  if (keyval->error() != KeyVal::OK) form_print_simples_ = 0;
  form_print_variable_ = keyval->booleanvalue("form:print_variable");
  if (keyval->error() != KeyVal::OK) form_print_variable_ = 0;
  form_print_constant_ = keyval->booleanvalue("form:print_constant");
  if (keyval->error() != KeyVal::OK) form_print_constant_ = 0;
  form_print_molecule_ = keyval->booleanvalue("form:print_molecule");
  if (keyval->error() != KeyVal::OK) form_print_molecule_ = 0;
}

void
IntMolecularCoor::init()
{
  Ref<SetIntCoor> redundant = new SetIntCoor;
  generator_->generate(redundant);

  // sort out the simple coordinates by type
  int i;
  for (i=0; i<redundant->n(); i++) {
      Ref<IntCoor> coor = redundant->coor(i);
      if (coor->class_desc()
          == ::class_desc<StreSimpleCo>()) {
          bonds_->add(coor);
        }
      else if (coor->class_desc() == ::class_desc<BendSimpleCo>()
               || coor->class_desc() == ::class_desc<LinIPSimpleCo>()
               || coor->class_desc() == ::class_desc<LinOPSimpleCo>()) {
          bends_->add(coor);
        }
      else if (coor->class_desc()
               == ::class_desc<TorsSimpleCo>()
               || coor->class_desc()
               == ::class_desc<ScaledTorsSimpleCo>()) {
          tors_->add(coor);
        }
      else if (coor->class_desc()
               == ::class_desc<OutSimpleCo>()) {
          outs_->add(coor);
        }
      else {
          extras_->add(coor);
        }
    }

  all_->add(bonds_);
  all_->add(bends_);
  all_->add(tors_);
  all_->add(outs_);
  all_->add(extras_);

  // don't let form_coordinates create new variables coordinates
  // if they were given by the user
  int keep_variable = (variable_->n() != 0);

  if (given_fixed_values_) {
      // save the given coordinate values
      RefSCDimension original_dfixed
          = new SCDimension(fixed_->n(),"Nfix");
      RefSCVector given_fixed_coords(original_dfixed,matrixkit_);
      for (i=0; i<original_dfixed.n(); i++) {
          given_fixed_coords(i) = fixed_->coor(i)->value();
        }

      // find the current fixed coordinates
      RefSCVector current_fixed_coords(original_dfixed,matrixkit_);
      fixed_->update_values(molecule_);
      for (i=0; i<original_dfixed.n(); i++) {
          current_fixed_coords(i) = fixed_->coor(i)->value();
        }

      // the difference between current fixed and desired fixed
      RefSCVector diff_fixed_coords = given_fixed_coords-current_fixed_coords;

      // break up the displacement into several manageable steps
      double maxabs = diff_fixed_coords.maxabs();
      int nstep = int(maxabs/max_update_disp_) + 1;
      diff_fixed_coords.scale(1.0/nstep);
      ExEnv::out0() << indent << "IntMolecularCoor: "
           << "displacing fixed coordinates to the requested values in "
           << nstep << " steps\n";
      for (int istep=1; istep<=nstep; istep++) {
          form_coordinates(keep_variable);

          dim_ = new SCDimension(variable_->n(), "Nvar");
          dvc_ = new SCDimension(variable_->n()+constant_->n(),
                             "Nvar+Nconst");

          RefSCVector new_internal_coordinates(dvc_,matrixkit_);
          for (i=0; i<variable_->n(); i++) {
              new_internal_coordinates(i) = variable_->coor(i)->value();
            }
          int j;
          for (j=0; j<original_dfixed.n(); j++,i++) {
              new_internal_coordinates(i)
                  = current_fixed_coords(j)+istep*double(diff_fixed_coords(j));
            }
          for (; j<constant_->n(); i++,j++) {
              new_internal_coordinates(i) = constant_->coor(j)->value();
            }

          all_to_cartesian(molecule_, new_internal_coordinates);
        }

      // make sure that the coordinates have exactly the
      // original values to avoid round-off error
      for (i=0; i<original_dfixed.n(); i++) {
          fixed_->coor(i)->set_value(given_fixed_coords(i));
        }
    }

  form_coordinates(keep_variable);

  dim_ = new SCDimension(variable_->n(), "Nvar");
  dvc_ = new SCDimension(variable_->n()+constant_->n(),
                         "Nvar+Nconst");

#if 0 // this will always think the rank has changed with redundant coordinates
    {
      const double epsilon = 0.001;

      // compute the condition number
      RefSCMatrix B(dim_, dnatom3_,matrixkit_);
      variable_->bmat(molecule_, B);