coor.h

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

    void add_bonds(const Ref<SetIntCoor>& list, BitArrayLTri& bonds, Molecule& m);
    void add_bends(const Ref<SetIntCoor>& list, BitArrayLTri& bonds, Molecule& m);
    void add_tors(const Ref<SetIntCoor>& list, BitArrayLTri& bonds, Molecule& m);
    void add_out(const Ref<SetIntCoor>& list, BitArrayLTri& bonds, Molecule& m);
  public:
    /** Create an IntCoorGen given a Molecule and, optionally, extra bonds.
        IntCoorGen keeps a reference to extra and deletes it when the
        destructor is called. */
    IntCoorGen(const Ref<Molecule>&, int nextra=0, int *extra=0);
    /** The KeyVal constructor.
        <dl>

        <dt><tt>molecule</tt><dd> A Molecule object.  There is no default.

        <dt><tt>radius_scale_factor</tt><dd> If the distance between two
        atoms is less than the radius scale factor times the sum of the
        atoms' atomic radii, then a bond is placed between the two atoms
        for the purpose of finding internal coordinates.  The default is
        1.1.

        <dt><tt>linear_bend_threshold</tt><dd> A bend angle in degress
        greater than 180 minus this keyword's floating point value is
        considered a linear bend. The default is 1.0.

        <dt><tt>linear_tors_threshold</tt><dd> The angles formed by atoms
        a-b-c and b-c-d are checked for near linearity.  If an angle in
        degrees is greater than 180 minus this keyword's floating point
        value, then the torsion is classified as a linear torsion. The
        default is 1.0.

        <dt><tt>linear_bend</tt><dd> Generate BendSimpleCo objects to
        describe linear bends.  The default is false.

        <dt><tt>linear_lbend</tt><dd> Generate pairs of LinIPSimpleCo and
        LinIPSimpleCo objects to describe linear bends.  The default is
        true.

        <dt><tt>linear_tors</tt><dd> Generate TorsSimpleCo objects to
        described linear torsions.  The default is false.

        <dt><tt>linear_stors</tt><dd> Generate ScaledTorsSimpleCo objects
        to described linear torsions.  The default is true.

        <dt><tt>extra_bonds</tt><dd> This is a vector of atom numbers,
        where elements \f$2 (i-1) + 1\f$ and \f$2 i\f$ specify the atoms
        which are bound in extra bond \f$i\f$.  The extra_bonds keyword
        should only be needed for weakly interacting fragments, otherwise
        all the needed bonds will be found.

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

    ~IntCoorGen();

    /// Standard member.
    void save_data_state(StateOut&);

    /// This generates a set of internal coordinates.
    virtual void generate(const Ref<SetIntCoor>&);

    /// Print out information about this.
    virtual void print(std::ostream& out=ExEnv::out0()) const;
};


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


/** The MolecularCoor abstract class describes the coordinate system used
to describe a molecule.  It is used to convert a molecule's cartesian
coordinates to and from this coordinate system. */
class MolecularCoor: public SavableState
{
  protected:
    Ref<Molecule> molecule_;
    RefSCDimension dnatom3_; // the number of atoms x 3
    Ref<SCMatrixKit> matrixkit_; // used to construct matrices

    int debug_;
  public:
    MolecularCoor(Ref<Molecule>&);
    MolecularCoor(StateIn&);
    /** The KeyVal constructor.
        <dl>

        <dt><tt>molecule</tt><dd> A Molecule object.  There is no default.

        <dt><tt>debug</tt><dd> An integer which, if nonzero, will cause
        extra output.

        <dt><tt>matrixkit</tt><dd> A SCMatrixKit object.  It is usually
        unnecessary to give this keyword.

        <dt><tt>natom3</tt><dd> An SCDimension object for the dimension of
        the vector of cartesian coordinates.  It is usually unnecessary to
        give this keyword.

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

    virtual ~MolecularCoor();

    void save_data_state(StateOut&);

    /** Returns a smart reference to an SCDimension equal to the
        number of atoms in the molecule times 3. */
    RefSCDimension dim_natom3() { return dnatom3_; }

    /// Returns the molecule.
    Ref<Molecule> molecule() const { return molecule_; }

    /// Print the coordinate.
    virtual void print(std::ostream& =ExEnv::out0()) const = 0;
    virtual void print_simples(std::ostream& =ExEnv::out0()) const = 0;

    /** Returns a smart reference to an SCDimension equal to the number of
        coordinates (be they Cartesian, internal, or whatever) that are
        being optimized. */
    virtual RefSCDimension dim() = 0;
    
    /** Given a set of displaced internal coordinates, update the cartesian
        coordinates of the Molecule contained herein.  This function does
        not change the vector ``internal''. */
    int to_cartesian(const RefSCVector&internal);
    virtual int to_cartesian(const Ref<Molecule>&mol,
                             const RefSCVector&internal) = 0;

    /** Fill in the vector ``internal'' with the current internal
        coordinates.  Note that this member will update the values of the
        variable internal coordinates. */
    virtual int to_internal(RefSCVector&internal) = 0;

    /** Convert the internal coordinate gradients in ``internal'' to
        Cartesian coordinates and copy these Cartesian coordinate gradients
        to ``cartesian''. Only the variable internal coordinate gradients
        are transformed. */
    virtual int to_cartesian(RefSCVector&cartesian,RefSCVector&internal) = 0;

    /** Convert the Cartesian coordinate gradients in ``cartesian'' to
        internal coordinates and copy these internal coordinate gradients
        to ``internal''.  Only the variable internal coordinate gradients
        are calculated. */
    virtual int to_internal(RefSCVector&internal,RefSCVector&cartesian) = 0;

    /** Convert the internal coordinate Hessian ``internal'' to Cartesian
        coordinates and copy the result to ``cartesian''.  Only the variable
        internal coordinate force constants are transformed. */
    virtual int to_cartesian(RefSymmSCMatrix&cartesian,
                              RefSymmSCMatrix&internal) =0;

    /** Convert the Cartesian coordinate Hessian ``cartesian'' to internal
        coordinates and copy the result to ``internal''.  Only the variable
        internal coordinate force constants are calculated. */
    virtual int to_internal(RefSymmSCMatrix&internal,
                             RefSymmSCMatrix&cartesian) = 0;

    /** Calculate an approximate hessian and place the result in
        ``hessian''. */
    virtual void guess_hessian(RefSymmSCMatrix&hessian) = 0;

    /** Given an Hessian, return the inverse of that hessian.  For singular
        matrices this should return the generalized inverse. */
    virtual RefSymmSCMatrix inverse_hessian(RefSymmSCMatrix&) = 0;

    /// Returns the number of constrained coordinates.
    virtual int nconstrained();

    /** When this is called, MoleculeCoor may select a new internal
        coordinate system and return a transform to it.  The default action
        is to not change anything and return an IdentityTransform. */
    virtual Ref<NonlinearTransform> change_coordinates();

    Ref<SCMatrixKit> matrixkit() const { return matrixkit_; }
};


/** The IntMolecularCoor abstract class describes a molecule's coordinates
in terms of internal coordinates. */
class IntMolecularCoor: public MolecularCoor
{
  protected:
    Ref<IntCoorGen> generator_;

    void form_K_matrix(RefSCDimension& dredundant,
                       RefSCDimension& dfixed,
                       RefSCMatrix& K,
                       int*& is_totally_symmetric);

    RefSCDimension dim_; // corresponds to the number of variable coordinates
    RefSCDimension dvc_; // the number of variable + constant coordinates

    Ref<SetIntCoor> variable_; // the variable internal coordinates
    Ref<SetIntCoor> constant_; // the constant internal coordinates
    
    Ref<SetIntCoor> fixed_;
    Ref<SetIntCoor> watched_;
    Ref<IntCoor> followed_;

    // these are all of the basic coordinates
    Ref<SetIntCoor> bonds_;
    Ref<SetIntCoor> bends_;
    Ref<SetIntCoor> tors_;
    Ref<SetIntCoor> outs_;
    // these are provided by the user or generated coordinates that
    // could not be assigned to any of the above catagories
    Ref<SetIntCoor> extras_;

    Ref<SetIntCoor> all_;

    // Useful relationships
    // variable_->n() + constant_->n() = 3N-6(5)
    // symm_->n() + asymm_->n() = 3N-6(5)

    int update_bmat_;  // if 1 recompute the b matrix during to_cartesian
    int only_totally_symmetric_; // only coors with tot. symm comp. are varied
    double symmetry_tolerance_; // tol used to find coors with tot. sym. comp.
    double simple_tolerance_; // tol used to see if a simple is included
    double coordinate_tolerance_; // tol used to see if a coor is included
    double cartesian_tolerance_;  // tol used in intco->cart transformation
    double scale_bonds_; // scale factor for bonds
    double scale_bends_; // scale factor for bends
    double scale_tors_;  // scale factor for tors
    double scale_outs_;  // scale factor for outs

    int nextra_bonds_;
    int* extra_bonds_;

    int given_fixed_values_; // if true make molecule have given fixed values

    int decouple_bonds_;
    int decouple_bends_;

    int max_update_steps_;
    double max_update_disp_;

    /** This is called by the constructors of classes derived from
        IntMolecularCoor.  It initialized the lists of simple internal
        coordinates, and then calls the form_coordinates() member. */
    virtual void init();
    /** Allocates memory for the SetIntCoor's used to store the
        simple and internal coordinates. */
    virtual void new_coords();
    /// Reads the KeyVal input.
    virtual void read_keyval(const Ref<KeyVal>&);

    // control whether or not to print coordinates when they are formed
    int form_print_simples_;
    int form_print_variable_;
    int form_print_constant_;
    int form_print_molecule_;
  public:
    IntMolecularCoor(StateIn&);
    IntMolecularCoor(Ref<Molecule>&mol);
    /** The KeyVal constructor.
        <dl>

        <dt><tt>variable</tt><dd> Gives a SetIntCoor object that specifies
        the internal coordinates that can be varied. If this is not given,
        the variable coordinates will be generated.

        <dt><tt>followed</tt><dd> Gives a IntCoor object that specifies a
        coordinate to used as the first coordinate in the variable