coor.h

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

//
// coor.h
//
// 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.
//

#ifndef _chemistry_molecule_coor_h
#define _chemistry_molecule_coor_h

#ifdef __GNUC__
#pragma interface
#endif

#include <iostream>
#include <vector>

#include <math/scmat/matrix.h>
#include <math/optimize/transform.h>
#include <chemistry/molecule/molecule.h>

namespace sc {

/** The IntCoor abstract class describes an internal coordinate of a
molecule. */
class IntCoor: public SavableState {
  protected:
    // conversion factors from radians, bohr to the preferred units
    static double bohr_conv;
    static double radian_conv;
    char *label_;
    double value_;
  public:
    IntCoor(StateIn&);
    IntCoor(const IntCoor&);
    /** This constructor takes a string containing a label for the
        internal coordinate.  The string is copied. */
    IntCoor(const char* label = 0);
    /** The KeyVal constructor.
        <dl>

        <dt><tt>label</tt><dd> A label for the coordinate using only to
        identify the coordinate to the user in printouts.  The default is
        no label.

        <dt><tt>value</tt><dd> A value for the coordinate.  In the way that
        coordinates are usually used, the default is to compute a value
        from the cartesian coordinates in a Molecule object.

        <dt><tt>unit</tt><dd> The unit in which the value is given.  This
        can be bohr, anstrom, radian, and degree.  The default is bohr for
        lengths and radian for angles.

        </dl> */
    IntCoor(const Ref<KeyVal>&);
    
    virtual ~IntCoor();
    void save_data_state(StateOut&);

    /// Returns the string containing the label for the internal coordinate.
    virtual const char* label() const;
    /// Returns the value of the coordinate in atomic units or radians.
    virtual double value() const;
    /// Sets the value of the coordinate in atomic units or radians.
    virtual void set_value(double);
    /// Returns the value of the coordinate in more familiar units.
    virtual double preferred_value() const;
    /// Returns a string representation of the type of coordinate this is.
    virtual const char* ctype() const = 0;
    /// Print information about the coordinate.
    virtual void print(std::ostream & o=ExEnv::out0()) const;
    virtual void print_details(const Ref<Molecule> &, std::ostream& =ExEnv::out0()) const;
    /** Returns the value of the force constant associated with this
        coordinate. */
    virtual double force_constant(Ref<Molecule>&) = 0;
    /// Recalculate the value of the coordinate.
    virtual void update_value(const Ref<Molecule>&) = 0;
    /// Fill in a row the the B matrix.
    virtual void bmat(const Ref<Molecule>&,RefSCVector&bmat,double coef=1.0) = 0;
    /** Test to see if this internal coordinate is equivalent to that one.
        The definition of equivalence is left up to the individual
        coordinates. */
    virtual int equivalent(Ref<IntCoor>&) = 0;
};

/** SumIntCoor is used to construct linear combinations of internal
coordinates.

The following is a sample ParsedKeyVal input for a SumIntCoor object:
<pre>
  sumintcoor<SumIntCoor>: (
    coor: [
      <StreSimpleCo>:( atoms = [ 1 2 ] )
      <StreSimpleCo>:( atoms = [ 2 3 ] )
      ]
    coef = [ 1.0 1.0 ]
    )
</pre>
*/
class SumIntCoor: public IntCoor {
  private:
    std::vector<double> coef_;
    std::vector<Ref<IntCoor> > coor_;
  public:
    SumIntCoor(StateIn&);
    /** This constructor takes a string containing a label for this
        coordinate. */
    SumIntCoor(const char *);
    /** The KeyVal constructor.
        <dl>

        <dt><tt>coor</tt><dd> A vector of IntCoor objects that define the
        summed coordinates.

        <dt><tt>coef</tt><dd> A vector of floating point numbers that gives
        the coefficients of the summed coordinates.

        </dl> */
    SumIntCoor(const Ref<KeyVal>&);

    ~SumIntCoor();
    void save_data_state(StateOut&);

    /// Returns the number of coordinates in this linear combination.
    int n();
    /** Add a coordinate to the linear combination.  coef is the
        coefficient for the added coordinate. */
    void add(Ref<IntCoor>&,double coef);
    /// This function normalizes all the coefficients.
    void normalize();

    // IntCoor overrides
    /// Returns the value of the coordinate in a.u. and radians.
    double preferred_value() const;
    /// Always returns ``SUM''.
    const char* ctype() const;
    /// Print the individual coordinates in the sum with their coefficients.
    void print_details(const Ref<Molecule> &, std::ostream& =ExEnv::out0()) const;
    /// Returns the weighted sum of the individual force constants.
    double force_constant(Ref<Molecule>&);
    /// Recalculate the value of the coordinate.
    void update_value(const Ref<Molecule>&);
    /// Fill in a row the the B matrix.
    void bmat(const Ref<Molecule>&,RefSCVector&bmat,double coef = 1.0);
    /// Always returns 0.
    int equivalent(Ref<IntCoor>&);
};

/** The SetIntCoor class describes a set of internal coordinates.
It can automatically generate these coordinates using a integral coordinate
generator object (see the IntCoorGen class) or the internal
coordinates can be explicity given.

The following is a sample ParsedKeyVal input for
a SetIntCoor object.
<pre>
  setintcoor<SetIntCoor>: [
    <SumIntCoor>: (
      coor: [
        <StreSimpleCo>:( atoms = [ 1 2 ] )
        <StreSimpleCo>:( atoms = [ 2 3 ] )
        ]
      coef = [ 1.0 1.0 ]
      )
    <BendSimpleCo>:( atoms = [ 1 2 3 ] )
  ]
</pre>
*/
class SetIntCoor: public SavableState {
  private:
    std::vector<Ref<IntCoor> > coor_;
  public:
    SetIntCoor();
    SetIntCoor(StateIn&);
    /** The KeyVal constructor.
        <dl>

        <dt><tt>generator</tt><dd> A IntCoorGen object that will be used to
        generate the internal coordinates.

        <dt><tt>i</tt><dd> A sequence of integer keywords, all \f$i\f$ for
        \f$0 \leq i < n\f$, can be assigned to IntCoor objects.

        </dl> */
    SetIntCoor(const Ref<KeyVal>&);

    virtual ~SetIntCoor();
    void save_data_state(StateOut&);

    /// Adds an internal coordinate to the set.
    void add(const Ref<IntCoor>&);
    /// Adds all the elements of another set to this one.
    void add(const Ref<SetIntCoor>&);
    /// Removes the last coordinate from this set.
    void pop();
    /// Removes all coordinates from the set.
    void clear();
    /// Returns the number of coordinates in the set.
    int n() const;
    /// Returns a reference to the i'th coordinate in the set.
    Ref<IntCoor> coor(int i) const;
    /// Compute the B matrix by finite displacements.
    virtual void fd_bmat(const Ref<Molecule>&,RefSCMatrix&);
    /// Compute the B matrix the old-fashioned way.
    virtual void bmat(const Ref<Molecule>&, RefSCMatrix&);
    /** Create an approximate Hessian for this set of coordinates.  This
        Hessian is a symmetric matrix whose i'th diagonal is the force
        constant for the i'th coordinate in the set. */
    virtual void guess_hessian(Ref<Molecule>&,RefSymmSCMatrix&);
    /// Print the coordinates in the set.
    virtual void print_details(const Ref<Molecule> &,std::ostream& =ExEnv::out0()) const;
    /// Recalculate the values of the internal coordinates in the set.
    virtual void update_values(const Ref<Molecule>&);
    /// Copy the values of the internal coordinates to a vector.
    virtual void values_to_vector(const RefSCVector&);
};


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

class BitArrayLTri;

/** IntCoorGen generates a set of simple internal coordinates
    for a molecule. */
class IntCoorGen: public SavableState
{
  protected:
    Ref<Molecule> molecule_;
    
    int linear_bends_;
    int linear_lbends_;
    int linear_tors_;
    int linear_stors_;
    int nextra_bonds_;
    int *extra_bonds_;
    double linear_bend_thres_;
    double linear_tors_thres_;
    double radius_scale_factor_;

    void init_constants();

    double cos_ijk(Molecule& m, int i, int j, int k);
    int hterminal(Molecule& m, BitArrayLTri& bonds, int i);
    int nearest_contact(int i, Molecule& m);