hamolmech.h

Go to the documentation of this file.

#ifndef HAMOLMECH_H
#define HAMOLMECH_H

#include "hastl.h"
#include "hastring.h"
#include "command.h"
#include "hacompmod.h"
#include "halinalg.h"
#include "vec3d.h"
#include "haatom.h"

#define AMBER_INT  1

class HaMolMechMod;
class AtomList;
class HaField3D;

enum SET_METHOD { NOT_SET = 0, SET_DEFAULT = 1, SET_FF_FIELD = 2, SET_RES_TEMPL = 3, SET_SPEC = 4 }; 

class MMBond
{
public:
        MMBond();
        MMBond(MatPoint* new_pt1, MatPoint* new_pt2);
        ~MMBond();

        bool operator==(const MMBond& rhs) const;
        bool operator< (const MMBond& rhs) const;

        double r0; // Equilibrium distance (in Ang)
        double fc; // force constant

        int set_type; 

public:
        MatPoint* pt1;
        MatPoint* pt2;

};

class MMNonBndParam
{
public:
        MMNonBndParam() {}
        virtual ~MMNonBndParam() {}

        double cn_6;
        double cn_12;

        bool operator == (const MMNonBndParam& rhs) const;
        bool operator <  (const MMNonBndParam& rhs) const;

};


class MMValAngle
{
public:
        MMValAngle();
        MMValAngle(MatPoint* new_pt1, MatPoint* new_pt2, 
                       MatPoint* new_pt3);
        ~MMValAngle();

public:
        MatPoint* pt1;
        MatPoint* pt2;
        MatPoint* pt3;

        double a0; 
        double fc; 

        int set_type; 

        bool operator == (const MMValAngle& rhs) const;
        bool operator <  (const MMValAngle& rhs) const;
};


class MMDihedral
{
public:
    MMDihedral();
        MMDihedral(MatPoint* new_pt1, MatPoint* new_pt2, 
                       MatPoint* new_pt3, MatPoint* new_pt4, bool improper_flag = false);
        virtual ~MMDihedral();

public:
        bool improper; 
        bool calc_14;  

        MatPoint* pt1;
        MatPoint* pt2;
        MatPoint* pt3;
        MatPoint* pt4;

        int set_type; 

// Torsional potential is expressed as:         
// E_tors = ( PK/IDIVF) * (1 + cos(PN*phi - PHASE))

        int GetNTerms() const { return pk.size(); } 
        void ClearParams() { pn.clear(); phase.clear(); pk.clear(); idivf.clear(); }  
        
        vector<double> pn;    
        vector<double> phase; 
        vector<double> pk;    
        vector<double> idivf; 

        int AddTerm( double pn_new, double phase_new, double pk_new, double idivf_new = 1.0 );

        bool operator == (const MMDihedral& rhs) const;
        bool operator <  (const MMDihedral& rhs) const;

};

typedef set<MatPoint*, less<MatPoint*> > SetMatPt;

// list of valid parameters for MM setup/execution:
const int MD_RUN = 0, MIN_RUN = 1, ENER_RUN = 2; // for run_type, IMIN in AMBER
const int NO_NMR_RESTRAIN = 0, READ_NMR_RESTRAIN=1, // for nmr_restrain
          NOESY_NMR_RESTRAIN = 2;                   // NMROPT in AMBER
const int CONJ_GRAD = 0, STEEPEST_DESCENT = 2; // for min_type - type of minimization, NTMIN in AMBER 
const int NO_PERIODICITY = 0, CONST_VOL = 1, CONST_PRES = 2; // for period_bcond, NTB in AMBER
const int NO_BORN = 0, PAIR_GEN_BORN = 1, VACUUM = 2, DIST_DEPEND_DIEL = 3; 
const int WRT_FORM = 0, WRT_BIN = 1; 
const int READ_X_FORM = 1, READ_X_BIN =2, READ_XV_BIN = 4, 
          READ_XV_FORM = 5, READ_XVBOX_BIN=6, READ_XVBOX_FORM = 7; // for init_read_coord =
const int NO_RESTART=0, DO_RESTART=1; // for restart_flag , IREST in AMBER
const int RESTRAINT_BIN=0, RESTRAINT_FORM =1; // for restraint_format , NTRX in AMBER 
const int WRITE_X_BIN = 0, WRITE_X_FORM = 1;  // for write_coord_format, NTXO in AMBER
const int CALC_ALL_INTER = 1, OMIT_BONDS_H = 2, OMIT_BONDS = 3 , // for omit_interactions
          OMIT_BONDS_VANG_H = 4, OMITS_BONDS_VANG = 5,           // NTF in AMBER
                  OMIT_BONDS_VANG_DIH_H = 6, OMIT_BONDS_VANG_DIH = 7; 
const int WATER_PAIRLIST_OXYGEN = 0, WATER_PAIRLIST_ATOMS = 86; // for water_pairlist_method, NTID in AMBER
const int DIST_DEP_DIEL = 0, CONST_DIEL = 1, GEN_BORN = 2; // for electr_model
const int NO_SOFT_REPULSION = 0, SOFT_REPULSION_NO_HBOND = 1, // for soft_repulsion
          SOFT_REPULSION_ALL = 2;           //  ISFTRP in AMBER
const int NO_POLAR = 0, ATOM_POLAR = 1, T3BODY_POLAR = 2; // for polar_method, IPOL in AMBER
const int MAXWELL_START_VEL = 3, READ_START_VEL = 4; // for start_vel_method, INIT in AMBER
const int CONST_ENE_MD = 0, CONST_TEMP_UNIFORM = 1,  // for temp_scale_method, NTT in AMBER
          CONST_TEMP_SOLUTE = 2, CONST_TEMP_THRESH = 3,
                  CONST_TEMP_QUICK_RESCALE = 4, CONST_TEMP_SOLUTE_SOLV = 5,
                  CONST_TEMP_MAXWELL = -1;
const int NO_CONST_P=0, ISOTROP_CONST_P=1, ANISOTROP_CONST_P=2; // for pressure_reg_method
const int ATOM_PRESS_SCALE = 0, MOL_PRESS_SCALE=1; // for pressure_scale_method, NPSCAL in AMBER
const int NO_SHAKE = 1, H_ATOM_SHAKE=2, ALL_ATOM_SHAKE=3; // for shake_constr, NTC in AMBER 
const int NORMAL_FAST_WATER = 0, REDEF_WATER_NAMES = 1, ONLY_SHAKE_FAST_WATER = 2,
          ONLY_SHAKE_FAST_WATER_REDEF_NAMES = 3, NO_FAST_WATER = 4; // for fast_water_method, JFASTW in AMBER
const int NMR_ABS_VAL_ERR = 1, NMR_SQ_SUM_ERR = 2, NMR_NOESY_PENALTY = 3; 

const int CALC_VDW_NO = 0, CALC_VDW_NORMAL = 1, CALC_VDW_NO_ATTRACT = 2; 
      
class CWinThread;     

const int SCRIPT_CONTINUE = 0;
const int SCRIPT_START = 1;
const int SCRIPT_STOP = 2;


//--------------------------------------------------

class HaMolMechMod : public HaCompMod
{
public:

        HaMolMechMod(HaMolSet* new_pmset=NULL);
        ~HaMolMechMod();

        int SetStdParams();
        int Initialize();        

        int SetStdValParams();       
        int SetStdVdWParams();       

        vector<SetMatPt> excluded_atom_list;    
        vector<SetMatPt> nonbond_contact_list;  

        bool BuildExcludedAtomList();    
        bool BuildNonBondContactList();  
        
    bool BuildGrpGrpNonBondList(AtomList* group1, AtomList* group2); 

        virtual int OnDelAtoms(AtomCollection& del_atoms); 
        bool DeletePoints( set< MatPoint*, less<MatPoint*> > *pt_set); 
        bool Print_info(ostream& sout, const int level);               

        vector<MatPoint*>    Points;          
        set<MMBond, less<MMBond> >       MBonds;          
        set<MMValAngle, less<MMValAngle> > ValAngles;       
        vector<MMDihedral>   Dihedrals;       
        vector<MMDihedral>   ImprDihedrals;   
        vector<MMBond>       BondConstrains;  

        MMBond*     GetMMBond(MatPoint* pt1, MatPoint* pt2); 
        MMValAngle* GetValAngle(MatPoint* pt1, MatPoint* pt2, MatPoint* pt3); 
        MMDihedral* GetDihedral(MatPoint* pt1, MatPoint* pt2, MatPoint* pt3,MatPoint* pt4); 
        MMDihedral* GetImprDihedral(MatPoint* pt1, MatPoint* pt2, MatPoint* pt3, MatPoint* pt4); 
        
        multimap<unsigned long, unsigned long, less<unsigned long> > res_impr_dih_map; 

        typedef vector<MatPoint*>  POINTS_TYPE;
        
        list<MMNonBndParam>    HBondParams;            



    static StrVecMap  symb_ppar_map;     
        static StrVecMap  bond_param_map;    
        static StrVecMap  vang_param_map;    
    static StrVecMap  dih_param_map;     
    static StrVecMap  impdih_param_map;  

        static HaVec_double  FindPointParamFromSymbol( const char* ats1);                  
        static HaVec_double  FindBondParamFromSymbol ( const char* ats1,const char* ats2);
        static HaVec_double  FindValAngleParamFromSymbol( const char* ats1,const char* ats2,const char* ats3);
        static HaVec_double  FindDihedralParamFromSymbol( const char* ats1,const char* ats2,const char* ats3,const char* ats4, bool improper_flag = false);

    static int init_ff_maps;  

        MMDihedral* AddImprDihedral(MatPoint* pt1, MatPoint* pt2, MatPoint* pt3, MatPoint* pt4);
        
        int SetMMBond(MatPoint* pt1,MatPoint* pt2,double r0,double fc,int set_type = NOT_SET); 
        int SetValAngle(MatPoint* pt1,MatPoint* pt2,MatPoint* pt3,double a0,double fc,int set_type = NOT_SET); 
        
        int GetResImprAngles(HaResidue* pres, list<MMDihedral*>& res_idih_list); 



        int Set14interDihFlags(); 





        int cur_nstep;  
        double time;       
        double temp;       
        double press;      

    double cur_energy;    

        double kin_ene;    
        double pot_ene;    
        double bond_ene;   
        double vang_ene;   
        double dihed_ene;  
        
        double vdw_ene;       
        double vdw_ene_14;    
        double vdw_ene_nb;    

        double electr_ene;    
        double electr_ene_14; 
        double electr_ene_nb; 
        double polar_ene;     

        double hbond_ene;     

        double constraint_ene; 

        double density;        
        
        double rms_ene;        
        double grad_ene_max;   

        int module_to_init_flag; 
        int build_nb_contact_list_flag; 
        int param_to_init_flag;  
        int init_charges_flag; 
        int InitParamSet();   



        HaString amber_run_file;  
        HaString amber_inp_file;  
        HaString amber_top_file;  
        HaString amber_crd_file;  

        HaString amber_out_file;        
        HaString amber_rst_file;        

        HaString amber_trj_coord_file;   
        HaString amber_trj_vel_file;     
        HaString amber_trj_ene_file;     
        HaString amber_constr_crd_file; 

        FILE* trj_coord_fp;   
        FILE* trj_vel_fp;  
        FILE* trj_ene_fp;  

        int SaveAmberRunFile();  
        int SaveAmberInpFile();  
        int SaveAmberTopFile();  
        int SaveAmberCrdFile();  
        int SaveAmberAllInpFiles();  



        int run_internal_flag;    
        int RunSanderInt();        

        int LoadAmberRestartFile(); 
        int LoadAmberMDInfoFile();  

        int RunAmber();              
        int RunAmberThread();         
        HaString sander_exe_fname;     

        int OpenSanderFilesToRead();  
        int ReadTrajPoint();          

        HaString traj_script;         

        long amber_proc_id;
        
        int StopCalc(); 

        CWinThread* run_amber_thread; 
        
        int delay_time;           
        int stop_calc_flag;       
        int skip_pt_init;         
        int skip_pt_between;      
        int play_back_flag;       

        int CalculateEnergy(); 

        bool CalcNonBondPt(const MatPoint* pt1, const MatPoint* pt2, double& vdw_at_ene, double& el_at_ene);  


    int time_limit;   
        int run_type;     
        int nmr_restrain; 

        int amber_version;  



    int init_read_coord;   
        int restart_flag;      
        int restraint_format;  



        int write_coord_format;   
        int print_freq;           
    int nbstp_wrt_rstrt;      
    int nbstp_wrt_coord;      
    int nbstp_wrt_vel;        
        int nbstp_wrt_ener;       
        int traj_wrt_format;      
        int limit_wrt_atoms;      

        int wrap_coord;           



    int omit_interactions;   
        int period_bcond;        
        int electr_model;        

        double scale_14_electr;   
        double scale_14_vdw;      
        double nonb_cutoff_dist;  
                                  
//      int diel_fun_type;        //!< the type of dielectric constant, IDIEL in AMBER, derived from electr_model
        double diel_const;        

        int nonb_list_flag;        
    int nonb_list_update_freq; 
        int water_pairlist_method; 
        int neutral_end_hydr_flag; 

        int calc_vdw_flag;        
        int calc_electr_flag;     
    
        int soft_repulsion;          
        double soft_repulsion_const; 


        int polar_method;    


    int num_3body_inter; 
    int num_ions;        

// \name 3-Body interaction parameters:
    HaVec_int atom2_3body;    
    HaVec_int atom3_3body;    
        HaVec_double pexp_3body;  
        HaVec_double beta_3body;  
        HaVec_double gamma_3body; 



        double atom_restr_const;        

        int SetMoveOnly( const HaString& new_mov_atom_list_name); 
        int SetMoveAll();                                         

        AtomList* GetMovingAtoms(); 

        HaString moving_atoms;  

        int SetRestrainedAtoms( const char* restr_atom_list_name); 
        AtomList* GetRestrAtoms(); 
        HaString restrained_atoms; 

        bool SetHBondConstraints(double force_const); 
        int AddHarmConstraint(HaAtom* atom1,HaAtom* atom2, double eq_dist, double force_const); 
        bool ClearConstraints(); 



        int min_type;       

    int max_comp_cycles;        
        int num_step_steep_descent; 

        double init_min_step;       
        double max_min_step;        
        double min_cnvrg_criterium; 



        int num_md_run;         
        int length_md_run;      
        int num_dfreedom_sub;   

    int remove_init_motion_flag; 


        int remove_motion_freq;      


        int start_vel_method;   
        double start_time;      
        double md_time_step;    



        double ref_temp;       
        double init_temp;      
    int random_seed;       
    double scale_init_vel; 
        int temp_scale_method; 


    int last_solute_atom;           
    double temp_deviation;          
        double temp_relax_time_solute;  
    double temp_relax_time_solvent; 
    double vel_limit;               



    double peaks_rate;        
        double peaks_coupling;    
        double peaks_init_pot;    



        int pressure_reg_method;    
        double ref_pressure;        
        double compressibility;     
        double pressure_relax_time; 
        int pressure_scale_method;  



        int shake_constr; 
        double shake_tol; 




        int solute_solvent_image_flag; 


        int remove_solute_nb_cut_flag; 

    int fast_water_method; 


        HaString wat_res_name; 
        HaString wat_ox_name;  
        HaString wat_h1_name;  
        HaString wat_h2_name;  



        int water_cap_flag; 
        int cap_atom_num;   
        double cap_fconst;  




        int noesy_vol_freq; 


        int nmr_penalty_min_method; 

        int nmr_max_num_submol; 




        int pmesh_ewald_flag;   
        double pme_box_x;       
        double pme_box_y;       // 
        double pme_box_z;       
        double pme_box_alpha;   
        double pme_box_beta;
        double pme_box_gamma;   

        int pme_grid_nx;   
        int pme_grid_ny;
        int pme_grid_nz;   

    int pme_spline_order;      
    int pme_neutralize_system; 
    int pme_verbose;           
        int exact_ewald_flag;      
        double pme_dsum_tol;       

protected:

};


#if defined(AMBER_INT)
extern "C" {
        void sandermn_(doublereal* x,integer* ix, integer* ih,
                            integer* ipairs, doublereal* r_stack, integer* i_stack); 

typedef struct{
        char mdin[80],   mdout[80],  inpcrd[80], parm[80],  restrt[80], 
             refc[80],   mdvel[80],  mden[80],   mdcrd[80], mdinfo[80], 
                     nmr[80],    mincor[80], vecs[80],   radii[80], freqe[80], owrite[80],
             rstdip[80], mddip[80],  inpdip[80];         
          } files_type;
 
          extern files_type files_;

};
#endif

#endif // end if !defined(HAMOLMECH_H) 

---