coor.h

//
// coor.h
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <cljanss@limitpt.com>
// 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
 
#include <iostream>
#include <vector>
 
#include <math/scmat/matrix.h>
#include <math/optimize/transform.h>
#include <chemistry/molecule/molecule.h>
 
namespace sc {
 
 
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&);
    IntCoor(const char* label = 0);
    IntCoor(const Ref<KeyVal>&);
    
    virtual ~IntCoor();
    void save_data_state(StateOut&);
 
    virtual const char* label() const;
    virtual double value() const;
    virtual void set_value(double);
    virtual double preferred_value() const;
    virtual const char* ctype() const = 0;
    virtual void print(std::ostream & o=ExEnv::out0()) const;
    virtual void print_details(const Ref<Molecule> &, std::ostream& =ExEnv::out0()) const;
    virtual double force_constant(Ref<Molecule>&) = 0;
    virtual void update_value(const Ref<Molecule>&) = 0;
    virtual void bmat(const Ref<Molecule>&,RefSCVector&bmat,double coef=1.0) = 0;
    virtual int equivalent(Ref<IntCoor>&) = 0;
};
 
class SumIntCoor: public IntCoor {
  private:
    std::vector<double> coef_;
    std::vector<Ref<IntCoor> > coor_;
  public:
    SumIntCoor(StateIn&);
    SumIntCoor(const char *);
    SumIntCoor(const Ref<KeyVal>&);
 
    ~SumIntCoor();
    void save_data_state(StateOut&);
 
    int n();
    void add(Ref<IntCoor>&,double coef);
    void normalize();
 
    // IntCoor overrides
    double preferred_value() const;
    const char* ctype() const;
    void print_details(const Ref<Molecule> &, std::ostream& =ExEnv::out0()) const;
    double force_constant(Ref<Molecule>&);
    void update_value(const Ref<Molecule>&);
    void bmat(const Ref<Molecule>&,RefSCVector&bmat,double coef = 1.0);
    int equivalent(Ref<IntCoor>&);
};
 
class SetIntCoor: public SavableState {
  private:
    std::vector<Ref<IntCoor> > coor_;
  public:
    SetIntCoor();
    SetIntCoor(StateIn&);
    SetIntCoor(const Ref<KeyVal>&);
 
    virtual ~SetIntCoor();
    void save_data_state(StateOut&);
 
    void add(const Ref<IntCoor>&);
    void add(const Ref<SetIntCoor>&);
    void pop();
    void clear();
    int n() const;
    Ref<IntCoor> coor(int i) const;
    virtual void fd_bmat(const Ref<Molecule>&,RefSCMatrix&);
    virtual void bmat(const Ref<Molecule>&, RefSCMatrix&);
    virtual void guess_hessian(Ref<Molecule>&,RefSymmSCMatrix&);
    virtual void print_details(const Ref<Molecule> &,std::ostream& =ExEnv::out0()) const;
    virtual void update_values(const Ref<Molecule>&);
    virtual void values_to_vector(const RefSCVector&);
};
 
 
// ////////////////////////////////////////////////////////////////////////
 
class BitArrayLTri;
 
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);
 
    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);
 
    void connect_subgraphs(const Molecule& mol,
                           BitArrayLTri& adjacency_matrix);
 
  public:
    IntCoorGen(const Ref<Molecule>&, int nextra=0, int *extra=0);
    IntCoorGen(const Ref<KeyVal>&);
    IntCoorGen(StateIn&);
 
    ~IntCoorGen();
 
    void save_data_state(StateOut&);
 
    virtual void generate(const Ref<SetIntCoor>&);
 
    virtual void print(std::ostream& out=ExEnv::out0()) const;
 
    static BitArrayLTri adjacency_matrix(const Molecule& mol,
                                         double radius_scaling_factor = 1.1);
    static std::vector<std::set<int> > find_disconnected_subgraphs(const BitArrayLTri& adjmat);
};
 
 
// ////////////////////////////////////////////////////////////////////////
 
 
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&);
    MolecularCoor(const Ref<KeyVal>&);
 
    virtual ~MolecularCoor();
 
    void save_data_state(StateOut&);
 
    RefSCDimension dim_natom3() { return dnatom3_; }
 
    Ref<Molecule> molecule() const { return molecule_; }
 
    virtual void print(std::ostream& =ExEnv::out0()) const = 0;
    virtual void print_simples(std::ostream& =ExEnv::out0()) const = 0;
 
    virtual RefSCDimension dim() = 0;
    
    int to_cartesian(const RefSCVector&internal);
    virtual int to_cartesian(const Ref<Molecule>&mol,
                             const RefSCVector&internal) = 0;
 
    virtual int to_internal(RefSCVector&internal) = 0;
 
    virtual int to_cartesian(RefSCVector&cartesian,RefSCVector&internal) = 0;
 
    virtual int to_internal(RefSCVector&internal,RefSCVector&cartesian) = 0;
 
    virtual int to_cartesian(RefSymmSCMatrix&cartesian,
                              RefSymmSCMatrix&internal) =0;
 
    virtual int to_internal(RefSymmSCMatrix&internal,
                             RefSymmSCMatrix&cartesian) = 0;
 
    virtual void guess_hessian(RefSymmSCMatrix&hessian) = 0;
 
    virtual RefSymmSCMatrix inverse_hessian(RefSymmSCMatrix&) = 0;
 
    virtual int nconstrained();
 
    virtual Ref<NonlinearTransform> change_coordinates();
 
    Ref<SCMatrixKit> matrixkit() const { return matrixkit_; }
};
 
 
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_;
 
    virtual void init();
    virtual void new_coords();
    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);
    IntMolecularCoor(const Ref<KeyVal>&);
 
    virtual ~IntMolecularCoor();
    void save_data_state(StateOut&);
  
    virtual void form_coordinates(int keep_variable=0) =0;
    
    virtual int all_to_cartesian(const Ref<Molecule> &,RefSCVector&internal);
    virtual int all_to_internal(const Ref<Molecule> &,RefSCVector&internal);
 
    virtual RefSCDimension dim();
    virtual int to_cartesian(const Ref<Molecule> &,const RefSCVector&internal);
    virtual int to_internal(RefSCVector&internal);
    virtual int to_cartesian(RefSCVector&cartesian,RefSCVector&internal);
    virtual int to_internal(RefSCVector&internal,RefSCVector&cartesian);
    virtual int to_cartesian(RefSymmSCMatrix&cart,RefSymmSCMatrix&internal);
    virtual int to_internal(RefSymmSCMatrix&internal,RefSymmSCMatrix&cart);
    virtual void print(std::ostream& =ExEnv::out0()) const;
    virtual void print_simples(std::ostream& =ExEnv::out0()) const;
    virtual void print_variable(std::ostream& =ExEnv::out0()) const;
    virtual void print_constant(std::ostream& =ExEnv::out0()) const;
    int nconstrained();
};
 
// ///////////////////////////////////////////////////////////////////////
 
class SymmMolecularCoor: public IntMolecularCoor
{
  protected:
    // true if coordinates should be changed during optimization
    int change_coordinates_;
    // true if hessian should be transformed too (should always be true)
    int transform_hessian_;
    // max value for the condition number if coordinates can be changed
    double max_kappa2_;
 
    void init();
  public:
    SymmMolecularCoor(Ref<Molecule>&mol);
    SymmMolecularCoor(StateIn&);
    SymmMolecularCoor(const Ref<KeyVal>&);
 
    virtual ~SymmMolecularCoor();
    void save_data_state(StateOut&);
 
    void form_coordinates(int keep_variable=0);
 
    void guess_hessian(RefSymmSCMatrix&hessian);
    RefSymmSCMatrix inverse_hessian(RefSymmSCMatrix&);
 
    Ref<NonlinearTransform> change_coordinates();
 
    void print(std::ostream& =ExEnv::out0()) const;
};
 
// ///////////////////////////////////////////////////////////////////////
 
class RedundMolecularCoor: public IntMolecularCoor
{
 
  public:
    RedundMolecularCoor(Ref<Molecule>&mol);
    RedundMolecularCoor(StateIn&);
    RedundMolecularCoor(const Ref<KeyVal>&);
 
    virtual ~RedundMolecularCoor();
    void save_data_state(StateOut&);
 
    void form_coordinates(int keep_variable=0);
    void guess_hessian(RefSymmSCMatrix&hessian);
    RefSymmSCMatrix inverse_hessian(RefSymmSCMatrix&);
};
 
// ///////////////////////////////////////////////////////////////////////
 
class CartMolecularCoor: public MolecularCoor
{
  private:
  protected:
    RefSCDimension dim_; // the number of atoms x 3
 
    virtual void init();
  public:
    CartMolecularCoor(Ref<Molecule>&mol);
    CartMolecularCoor(StateIn&);
    CartMolecularCoor(const Ref<KeyVal>&);
 
    virtual ~CartMolecularCoor();
  
    void save_data_state(StateOut&);
 
    virtual RefSCDimension dim();
    virtual int to_cartesian(const Ref<Molecule>&,const RefSCVector&internal);
    virtual int to_internal(RefSCVector&internal);
    virtual int to_cartesian(RefSCVector&cartesian,RefSCVector&internal);
    virtual int to_internal(RefSCVector&internal,RefSCVector&cartesian);
    virtual int to_cartesian(RefSymmSCMatrix&cart,RefSymmSCMatrix&internal);
    virtual int to_internal(RefSymmSCMatrix&internal,RefSymmSCMatrix&cart);
    virtual void print(std::ostream& =ExEnv::out0()) const;
    virtual void print_simples(std::ostream& =ExEnv::out0()) const;
    void guess_hessian(RefSymmSCMatrix&hessian);
    RefSymmSCMatrix inverse_hessian(RefSymmSCMatrix&);
};
 
// end of addtogroup ChemistryMolecule
 
}
 
#endif
 
// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: