Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

Parameters.h

Go to the documentation of this file.
00001 
00007 /*
00008    The Parameters class is used to read in and store all of the parameters
00009    from the parameter files.  These parameters are stored and used to assign
00010    constants to bonds and atoms as they are read in from the psf file
00011 */
00012 
00013 #ifndef PARAM_H
00014 
00015 #define PARAM_H
00016 
00017 #include "parm.h"
00018 
00019 #include "common.h"
00020 #include "structures.h"
00021 #include "strlib.h"
00022 #include "MStream.h"
00023 //****** BEGIN CHARMM/XPLOR type changes
00024 #include "SimParameters.h"
00025 //****** END CHARMM/XPLOR type changes
00026 
00027 #include "GromacsTopFile.h"
00028 
00029 #ifndef cbrt
00030   // cbrt() not in math.h on goneril
00031   #define cbrt(x)  pow((double)x,(double)(1.0/3.0))
00032 #endif
00033 
00034 class Communicate;
00035 class StringList;
00036 
00037 //  The class Parameters is used to store and query the parameters for
00038 //  bonds and atoms.  If the Parameter object resides on the master
00039 //  process, it is responsible for reading in all the parameters and
00040 //  then communicating them to the other processors.  To do this, first
00041 //  the routine read_paramter_file is called to read in each parameter
00042 //  file.  The information that is read in is stored in binary trees
00043 //  for vdw, bond, and angle information and in linked lists for 
00044 //  dihedrals and impropers.  Once all of the files have been read
00045 //  in, the routine done_reading_files is called.  At this point, all
00046 //  of the data that has been read in is copied to arrays.  This is
00047 //  so that each bond and atom can be assigned an index into these
00048 //  arrays to retreive the parameters in constant time.
00049 //
00050 //  Then the psf file is read in.  Each bond and atom is assigned an
00051 //  index into the parameter lists using the functions assign_*_index.
00052 //  Once the psf file has been read in, then the routine 
00053 //  done_reading_structre is called to tell the object that it no
00054 //  longer needs to store the binary trees and linked lists that were
00055 //  used to query the parameters based on atom type.  From this point
00056 //  on, only the indexes will be used.
00057 //
00058 //  The master node then uses send_parameters to send all of these
00059 //  parameters to the other nodes and the objects on all of the other
00060 //  nodes use receive_parameters to accept these parameters.
00061 //
00062 //  From this point on, all of the parameter data is static and the
00063 //  functions get_*_params are used to retrieve the parameter data
00064 //  that is desired.
00065 
00066 //  Define the number of Multiples that a Dihedral or Improper
00067 //  bond can have with the Charm22 parameter set
00068 #define MAX_MULTIPLICITY        6
00069 
00070 //  Number of characters maximum allowed for storing atom type names
00071 #define MAX_ATOMTYPE_CHARS      6
00072 
00073 //****** BEGIN CHARMM/XPLOR type changes
00074 //  Define the numbers associated with each possible parameter-file 
00075 //  type (format) NAMD can handle.
00076 #define paraXplor               0
00077 #define paraCharmm              1
00078 //****** END CHARMM/XPLOR type changes
00079 
00080 
00081 class BondValue {
00082 public:
00083         Real k;         //  Force constant for the bond
00084         Real x0;        //  Rest distance for the bond
00085 };
00086 
00087 class AngleValue {
00088 public:
00089         Real k;         //  Force constant for angle
00090         Real theta0;    //  Rest angle for angle
00091         Real k_ub;      //  Urey-Bradley force constant
00092         Real r_ub;      //  Urey-Bradley distance
00093   int normal; // Whether we use harmonic (0) or cos-based (1) angle terms
00094 };
00095 
00096 typedef struct four_body_consts
00097 {
00098         Real k;         //  Force constant
00099         int n;          //  Periodicity
00100         Real delta;     //  Phase shift
00101 } FourBodyConsts;
00102 
00103 class DihedralValue {
00104 public:
00105         int multiplicity;
00106         FourBodyConsts values[MAX_MULTIPLICITY];
00107 };
00108 
00109 class ImproperValue {
00110 public:
00111         int multiplicity;
00112         FourBodyConsts values[MAX_MULTIPLICITY];
00113 };
00114 
00115 struct CrosstermData { BigReal d00,d01,d10,d11; };
00116 
00117 class CrosstermValue {
00118 public:
00119         enum {dim=25};
00120         CrosstermData c[dim][dim];  // bicubic interpolation coefficients
00121 };
00122 
00123 class NonbondedExclValue {
00124 public:
00125         // need to put parameters here...
00126         // for now, copy bond
00127         Real k;         //  Force constant for the bond
00128         Real x0;        //  Rest distance for the bond
00129 };
00130 
00131 typedef struct vdw_val
00132 {
00133         Real sigma;     //  Sigma value
00134         Real epsilon;   //  Epsilon value
00135         Real sigma14;   //  Sigma value for 1-4 interactions
00136         Real epsilon14; //  Epsilon value for 1-4 interactions
00137 } VdwValue;
00138 
00139 //  IndexedVdwPair is used to form a binary search tree that is
00140 //  indexed by vwd_type index.  This is the tree that will be
00141 //  used to search during the actual simulation
00142 
00143 typedef struct indexed_vdw_pair
00144 {
00145         Index ind1;             //  Index for first atom type
00146         Index ind2;             //  Index for second atom type
00147         Real A;                 //  Parameter A for this pair
00148         Real A14;               //  Parameter A for 1-4 interactions
00149         Real B;                 //  Parameter B for this pair
00150         Real B14;               //  Parameter B for 1-4 interactions
00151         struct indexed_vdw_pair *right;  //  Right child
00152         struct indexed_vdw_pair *left;   //  Left child
00153 } IndexedVdwPair;
00154 
00155 
00156 //  IndexedTablePair is used to form a binary search tree that is
00157 //  indexed by table_type index.  This is the tree that will be
00158 //  used to search during the actual simulation
00159 
00160 typedef struct indexed_table_pair
00161 {
00162         Index ind1;             //  Index for first atom type
00163         Index ind2;             //  Index for second atom type
00164         int K;                  //  number of the table type for this pair
00165         struct indexed_table_pair *right;        //  Right child
00166         struct indexed_table_pair *left;         //  Left child
00167 } IndexedTablePair;
00168 
00169 //  Structures that are defined in Parameters.C
00170 struct bond_params;
00171 struct angle_params;
00172 struct improper_params;
00173 struct dihedral_params;
00174 struct crossterm_params;
00175 struct vdw_params;
00176 struct vdw_pair_params;
00177 struct table_pair_params;
00178 
00179 class Parameters
00180 {
00181 private:
00182         void initialize();                      //  zeros out pointers
00183 
00184         char *atomTypeNames;                    //  Names of atom types
00185         Bool AllFilesRead;                      //  Flag TRUE imples that all
00186                                                 //  of the parameter files
00187                                                 //  have been read in and
00188                                                 //  the arrays have been
00189                                                 //  created.
00190 //****** BEGIN CHARMM/XPLOR type changes
00191         int paramType;                          //  Type (format) of parameter-file
00192 //****** END CHARMM/XPLOR type changes
00193   bool cosAngles; // True if some angles may be cos-based
00194         struct bond_params *bondp;              //  Binary tree of bond params
00195         struct angle_params *anglep;            //  Binary tree of angle params
00196         struct improper_params *improperp;      //  Linked list of improper par.
00197         struct dihedral_params *dihedralp;      //  Linked list of dihedral par.
00198         struct crossterm_params *crosstermp;    //  Linked list of cross-term par.
00199         struct vdw_params *vdwp;                //  Binary tree of vdw params
00200         struct vdw_pair_params *vdw_pairp;      //  Binary tree of vdw pairs
00201         struct table_pair_params *table_pairp;  //  Binary tree of table pairs
00202 public:
00203         BondValue *bond_array;                  //  Array of bond params
00204         AngleValue *angle_array;                //  Array of angle params
00205         DihedralValue *dihedral_array;          //  Array of dihedral params
00206         ImproperValue *improper_array;          //  Array of improper params
00207         CrosstermValue *crossterm_array;        //  Array of crossterm params
00208         VdwValue *vdw_array;                    //  Array of vdw params
00209         int numenerentries;                     //  Number of entries for enertable
00210         int rowsize;
00211         int columnsize;
00212         BigReal* table_ener;                    //  Table for tabulated energies
00213         IndexedVdwPair *vdw_pair_tree;          //  Tree of vdw pair params
00214         IndexedTablePair *tab_pair_tree;                //  Tree of vdw pair params
00215   int tablenumtypes;
00216         int NumBondParams;                      //  Number of bond parameters
00217         int NumAngleParams;                     //  Number of angle parameters
00218   int NumCosAngles;       // Number of cosine-based angles
00219         int NumDihedralParams;                  //  Number of dihedral params
00220         int NumImproperParams;                  //  Number of improper params
00221         int NumCrosstermParams;                 //  Number of cross-term params
00222         int NumVdwParams;                       //  Number of vdw parameters
00223         int NumTableParams;                     //  Number of table parameters
00224   int NumVdwParamsAssigned;               //  Number actually assigned
00225         int NumVdwPairParams;                   //  Number of vdw_pair params
00226         int NumTablePairParams;                 //  Number of vdw_pair params
00227 private:
00228         ResizeArray<char *> error_msgs;         //  Avoids repeating warnings
00229 
00230         int *maxDihedralMults;                  //  Max multiplicity for
00231                                                 //  dihedral bonds
00232         int *maxImproperMults;                  //  Max multiplicity for
00233                                                 //  improper bonds
00234 
00235         void skip_stream_read(char *, FILE *);  // skip part of stream file
00236 
00237         void add_bond_param(char *);            //  Add a bond parameter
00238         struct bond_params *add_to_bond_tree(struct bond_params * , 
00239                                      struct bond_params *);
00240 
00241         void add_angle_param(char *);           //  Add an angle parameter
00242         struct angle_params *add_to_angle_tree(struct angle_params * , 
00243                                      struct angle_params *);
00244 
00245         void add_dihedral_param(char *, FILE *); //  Add a dihedral parameter
00246         void add_to_dihedral_list(struct dihedral_params *);
00247         void add_to_charmm_dihedral_list(struct dihedral_params *);
00248 
00249         void add_improper_param(char *, FILE *); //  Add an improper parameter
00250         void add_to_improper_list(struct improper_params *);
00251 
00252         void add_crossterm_param(char *, FILE *); //  Add an cross-term parameter
00253         void add_to_crossterm_list(struct crossterm_params *);
00254 
00255         void add_vdw_param(char *);             //  Add a vdw parameter
00256         struct vdw_params *add_to_vdw_tree(struct vdw_params *, 
00257                                      struct vdw_params *);
00258 
00259         void add_vdw_pair_param(char *);        //  Add a vdw pair parameter
00260         void add_table_pair_param(char *);      //  Add a table pair parameter
00261         void add_to_vdw_pair_list(struct vdw_pair_params *);
00262         void add_to_table_pair_list(struct table_pair_params *);
00263 
00264         void add_hb_pair_param(char *); //  Add a hydrogen bond pair parameter
00265 
00266         //  All of the traverse routines are used for debugging purposes
00267         //  to print out the appropriate list of parameters
00268         void traverse_vdw_pair_params(struct vdw_pair_params *);
00269         void traverse_vdw_params(struct vdw_params *);
00270         void traverse_dihedral_params(struct dihedral_params *);
00271         void traverse_improper_params(struct improper_params *);
00272         void traverse_crossterm_params(struct crossterm_params *);
00273         void traverse_angle_params(struct angle_params *);
00274         void traverse_bond_params(struct bond_params *);
00275 
00276         //  The index_* routines are used to index each of 
00277         //  the parameters and build the arrays that will be used
00278         //  for constant time access
00279         Index index_bonds(struct bond_params *, Index);
00280         Index index_angles(struct angle_params *, Index);
00281         Index index_vdw(struct vdw_params *, Index);
00282         void index_dihedrals();
00283         void index_impropers();
00284         void index_crossterms();
00285         
00286         void convert_vdw_pairs();
00287         void convert_table_pairs();
00288         IndexedVdwPair *add_to_indexed_vdw_pairs(IndexedVdwPair *, IndexedVdwPair *);
00289         IndexedTablePair *add_to_indexed_table_pairs(IndexedTablePair *, IndexedTablePair *);
00290         
00291         int vdw_pair_to_arrays(int *, int *, Real *, Real *, Real *, Real *, 
00292                                int, IndexedVdwPair *);
00293 
00294         int table_pair_to_arrays(int *, int *, int *, int, IndexedTablePair *);
00295 
00296         //  The free_* routines are used by the destructor to deallocate
00297         //  memory
00298         void free_bond_tree(struct bond_params *);
00299         void free_angle_tree(struct angle_params *);
00300         void free_dihedral_list(struct dihedral_params *);
00301         void free_improper_list(struct improper_params *);
00302         void free_crossterm_list(struct crossterm_params *);
00303         void free_vdw_tree(struct vdw_params *);
00304         void free_vdw_pair_tree(IndexedVdwPair *);
00305         void free_table_pair_tree(IndexedTablePair *);
00306         void free_vdw_pair_list();
00307 
00308   BigReal interp_lin(BigReal, BigReal, BigReal, BigReal, BigReal); // Perform a linear interpolation for energy table 
00309 
00310         /* does the actual initialization, once the variables have all
00311            been given default values.  See Parameters() below */
00312         void read_parm(const GromacsTopFile *, Bool min);
00313 public:
00314         //****** BEGIN CHARMM/XPLOR type changes
00316         Parameters();
00317         Parameters(SimParameters *, StringList *f);
00318         //****** END CHARMM/XPLOR type changes
00319         
00320         Parameters(Ambertoppar *, BigReal);
00321         void read_parm(Ambertoppar *, BigReal);
00322 
00323         /* initializes this to hold the set of parameters in the
00324            GROMACS topology file <gf>.  If the parameter <min> is on,
00325            this assumes that we are going to do minimization and
00326            therefore can't have atoms with zero VDW - it will add a
00327            small repulsive term to these. */
00328         Parameters(const GromacsTopFile *gf, Bool min);
00329         
00330         ~Parameters();                          //  Destructor
00331 
00332         // return a string for the Nth atom type.  This can only be
00333         // called after all the param files have been read and the type
00334         // names have been indexed.  The Nth atom type refers to the same
00335         // index of the Nth vdw parameter (i.e. there are NumVdwParams names).
00336         char *atom_type_name(Index a) {
00337           return (atomTypeNames + (a * (MAX_ATOMTYPE_CHARS + 1)));
00338         }
00339 
00340         //  Read a parameter file
00341         void read_parameter_file(char *);
00342 
00343         //****** BEGIN CHARMM/XPLOR type changes
00344         void read_charmm_parameter_file(char *);
00345         //****** END CHARMM/XPLOR type changes
00346 
00347         //  Signal the parameter object that all of
00348         //  the parameter files have been read in
00349         void done_reading_files();
00350 
00351         //  Signal the parameter object that the
00352         //  structure file has been read in
00353         void done_reading_structure();
00354     
00355         //  The assign_*_index routines are used to assign
00356         //  an index to atoms or bonds.  If an specific atom
00357         //  or bond type can't be found, then the program 
00358         //  terminates
00359     #ifdef MEM_OPT_VERSION
00360     void assign_vdw_index(char *, AtomCstInfo *);
00361     #else
00362         void assign_vdw_index(char *, Atom *);  //  Assign a vdw index to
00363     #endif
00364                                                 //  an atom
00365         void assign_bond_index(char *, char *, Bond *); 
00366                                                 //  Assign a bond index
00367                                                 //  to a bond
00368         void assign_angle_index(char *, char *, char *, Angle *);
00369                                                 //  Assign an angle index
00370                                                 //  to an angle
00371         void assign_dihedral_index(char *, char*, char*, char *, Dihedral *, int);
00372                                                 //  Assign a dihedral index
00373                                                 //  to a dihedral
00374         void assign_improper_index(char *, char*, char*, char *, Improper *, int);
00375                                                 //  Assign an improper index
00376                                                 //  to an improper
00377         void assign_crossterm_index(char *, char*, char*, char *, char *, char*, char*, char *, Crossterm *);
00378 
00379         //  send_parameters is used by the master process to
00380         //  communicate the paramters to all the other processors
00381         void send_Parameters(MOStream *);
00382 
00383         //  receive_parameters is used by all the child processes
00384         //  to receive the parameters from the master process
00385         void receive_Parameters(MIStream *);
00386 
00387         //  The get_*_params routines are the routines that really
00388         //  do all the work for this object.  Given an index, they
00389         //  access the parameters and return the relevant information
00390         void get_bond_params(Real *k, Real *x0, Index index)
00391         {
00392                 *k = bond_array[index].k;
00393                 *x0 = bond_array[index].x0;
00394         }
00395 
00396         void get_angle_params(Real *k, Real *theta0, Real *k_ub, Real *r_ub,
00397                               Index index)
00398         {
00399                 *k = angle_array[index].k;
00400                 *theta0 = angle_array[index].theta0;
00401                 *k_ub = angle_array[index].k_ub;
00402                 *r_ub = angle_array[index].r_ub;
00403         }
00404 
00405         int get_improper_multiplicity(Index index)
00406         {
00407                 return(improper_array[index].multiplicity);
00408         }
00409 
00410         int get_dihedral_multiplicity(Index index)
00411         {
00412                 return(dihedral_array[index].multiplicity);
00413         }
00414 
00415         void get_improper_params(Real *k, int *n, Real *delta, 
00416                                  Index index, int mult)
00417         {
00418                 if ( (mult<0) || (mult>MAX_MULTIPLICITY) )
00419                 {
00420                         NAMD_die("Bad mult index in Parameters::get_improper_params");
00421                 }
00422 
00423                 *k = improper_array[index].values[mult].k;
00424                 *n = improper_array[index].values[mult].n;
00425                 *delta = improper_array[index].values[mult].delta;
00426         }
00427 
00428         void get_dihedral_params(Real *k, int *n, Real *delta, 
00429                                  Index index, int mult)
00430         {
00431                 if ( (mult<0) || (mult>MAX_MULTIPLICITY) )
00432                 {
00433                         NAMD_die("Bad mult index in Parameters::get_dihedral_params");
00434                 }
00435 
00436                 *k = dihedral_array[index].values[mult].k;
00437                 *n = dihedral_array[index].values[mult].n;
00438                 *delta = dihedral_array[index].values[mult].delta;
00439         }
00440 
00441         void get_vdw_params(Real *sigma, Real *epsilon, Real *sigma14, 
00442                             Real *epsilon14, Index index)
00443   {
00444     if ( vdw_array ) {
00445       *sigma = vdw_array[index].sigma;
00446       *epsilon = vdw_array[index].epsilon;
00447       *sigma14 = vdw_array[index].sigma14;
00448       *epsilon14 = vdw_array[index].epsilon14;
00449     } else {
00450       // sigma and epsilon from A and B for each vdw type's interaction with itself
00451       Real A; Real B; Real A14; Real B14;
00452       if (get_vdw_pair_params(index, index, &A, &B, &A14, &B14) ) {
00453         if (A && B) {
00454           *sigma = sqrt(cbrt(A)) / sqrt(cbrt(B));
00455           *epsilon = B*B / (4*A);
00456         }
00457         else {
00458           *sigma = 0; *epsilon = 0;
00459         }
00460         if (A14 && B14) {
00461           *sigma14 = sqrt(cbrt(A14)) / sqrt(cbrt(B14));
00462           *epsilon14 = B14*B14 / (4*A14);
00463         }
00464         else {
00465           *sigma14 = 0; *epsilon14 = 0;
00466         }
00467       }
00468       else {NAMD_die ("Function get_vdw_params failed to derive Lennard-Jones sigma and epsilon from A and B values\n");}
00469     }
00470   }
00471 
00472         int get_vdw_pair_params(Index ind1, Index ind2, Real *, Real *, Real *, Real *);
00473                                                 //  Find a vwd_pair parameter
00474 
00475         int get_num_vdw_params(void) { return NumVdwParamsAssigned; }
00476 
00477   int get_table_pair_params(Index, Index, int*);
00478 
00479         //  The print_*_params are provided for debugging purposes
00480         void print_bond_params();               //  Print bonds params
00481         void print_angle_params();              //  Print angle params
00482         void print_dihedral_params();           //  Print dihedral params
00483         void print_improper_params();           //  Print improper params
00484         void print_crossterm_params();          //  Print cross-term params
00485         void print_vdw_params();                //  Print vdw params
00486         void print_vdw_pair_params();           //  Print vdw_pair params
00487         void print_param_summary();             //  Print a summary of params
00488         void read_ener_table(SimParameters*); // Read an energy table file
00489   int get_int_table_type(char*); // Return the integer type for a named table interaction
00490 
00491   int read_energy_type(FILE*, const int, BigReal*, const float, const float);
00492   int read_energy_type_cubspline(FILE*, const int, BigReal*, const float, const float);
00493   int read_energy_type_bothcubspline(FILE*, const int, BigReal*, const float, const float);
00494 };
00495 
00496 #endif
00497 

Generated on Mon Nov 23 04:59:22 2009 for NAMD by  doxygen 1.3.9.1