3f7b9b77680b555a2cc5e7f27704560874526b62
[namd.git] / src / Parameters.h
1 /**
2 ***  Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
3 ***  The Board of Trustees of the University of Illinois.
4 ***  All rights reserved.
5 **/
6
7 /*
8    The Parameters class is used to read in and store all of the parameters
9    from the parameter files.  These parameters are stored and used to assign
10    constants to bonds and atoms as they are read in from the psf file
11 */
12
13 #ifndef PARAM_H
14
15 #define PARAM_H
16
17 #include "parm.h"
18
19 #include "common.h"
20 #include "structures.h"
21 #include "strlib.h"
22 #include "MStream.h"
23 //****** BEGIN CHARMM/XPLOR type changes
24 #include "SimParameters.h"
25 //****** END CHARMM/XPLOR type changes
26
27 #include "GromacsTopFile.h"
28
29 #ifndef cbrt
30   // cbrt() not in math.h on goneril
31   #define cbrt(x)  pow((double)x,(double)(1.0/3.0))
32 #endif
33
34 class Communicate;
35 class StringList;
36
37 //  The class Parameters is used to store and query the parameters for
38 //  bonds and atoms.  If the Parameter object resides on the master
39 //  process, it is responsible for reading in all the parameters and
40 //  then communicating them to the other processors.  To do this, first
41 //  the routine read_paramter_file is called to read in each parameter
42 //  file.  The information that is read in is stored in binary trees
43 //  for vdw, bond, and angle information and in linked lists for 
44 //  dihedrals and impropers.  Once all of the files have been read
45 //  in, the routine done_reading_files is called.  At this point, all
46 //  of the data that has been read in is copied to arrays.  This is
47 //  so that each bond and atom can be assigned an index into these
48 //  arrays to retreive the parameters in constant time.
49 //
50 //  Then the psf file is read in.  Each bond and atom is assigned an
51 //  index into the parameter lists using the functions assign_*_index.
52 //  Once the psf file has been read in, then the routine 
53 //  done_reading_structre is called to tell the object that it no
54 //  longer needs to store the binary trees and linked lists that were
55 //  used to query the parameters based on atom type.  From this point
56 //  on, only the indexes will be used.
57 //
58 //  The master node then uses send_parameters to send all of these
59 //  parameters to the other nodes and the objects on all of the other
60 //  nodes use receive_parameters to accept these parameters.
61 //
62 //  From this point on, all of the parameter data is static and the
63 //  functions get_*_params are used to retrieve the parameter data
64 //  that is desired.
65
66 //  Define the number of Multiples that a Dihedral or Improper
67 //  bond can have with the Charm22 parameter set
68 #define MAX_MULTIPLICITY        6
69
70 //  Number of characters maximum allowed for storing atom type names
71 #define MAX_ATOMTYPE_CHARS      6
72
73 //****** BEGIN CHARMM/XPLOR type changes
74 //  Define the numbers associated with each possible parameter-file 
75 //  type (format) NAMD can handle.
76 #define paraXplor               0
77 #define paraCharmm              1
78 //****** END CHARMM/XPLOR type changes
79
80
81 class BondValue {
82 public:
83         Real k;         //  Force constant for the bond
84         Real x0;        //  Rest distance for the bond
85 };
86
87 class AngleValue {
88 public:
89         Real k;         //  Force constant for angle
90         Real theta0;    //  Rest angle for angle
91         Real k_ub;      //  Urey-Bradley force constant
92         Real r_ub;      //  Urey-Bradley distance
93   int normal; // Whether we use harmonic (1) or cos-based (0) angle terms
94 };
95
96 typedef struct four_body_consts
97 {
98         Real k;         //  Force constant
99         int n;          //  Periodicity
100         Real delta;     //  Phase shift
101 } FourBodyConsts;
102
103 class DihedralValue {
104 public:
105         int multiplicity;
106         FourBodyConsts values[MAX_MULTIPLICITY];
107 };
108
109 class ImproperValue {
110 public:
111         int multiplicity;
112         FourBodyConsts values[MAX_MULTIPLICITY];
113 };
114
115 struct CrosstermData { BigReal d00,d01,d10,d11; };
116
117 class CrosstermValue {
118 public:
119         enum {dim=25};
120         CrosstermData c[dim][dim];  // bicubic interpolation coefficients
121 };
122
123 // JLai
124 class GromacsPairValue {
125 public:
126     BigReal pairC6;
127     BigReal pairC12;
128 /*    BigReal A;
129     BigReal B;
130     BigReal A14;
131     BigReal B14;
132     BigReal sqrt26;
133     BigReal expo;*/
134 };
135 // End of JLai
136
137
138 class NonbondedExclValue {
139 public:
140         // need to put parameters here...
141         // for now, copy bond
142         Real k;         //  Force constant for the bond
143         Real x0;        //  Rest distance for the bond
144 };
145
146 typedef struct vdw_val
147 {
148         Real sigma;     //  Sigma value
149         Real epsilon;   //  Epsilon value
150         Real sigma14;   //  Sigma value for 1-4 interactions
151         Real epsilon14; //  Epsilon value for 1-4 interactions
152 } VdwValue;
153
154 //  IndexedVdwPair is used to form a binary search tree that is
155 //  indexed by vwd_type index.  This is the tree that will be
156 //  used to search during the actual simulation
157
158 typedef struct indexed_vdw_pair
159 {
160         Index ind1;             //  Index for first atom type
161         Index ind2;             //  Index for second atom type
162         Real A;                 //  Parameter A for this pair
163         Real A14;               //  Parameter A for 1-4 interactions
164         Real B;                 //  Parameter B for this pair
165         Real B14;               //  Parameter B for 1-4 interactions
166         struct indexed_vdw_pair *right;  //  Right child
167         struct indexed_vdw_pair *left;   //  Left child
168 } IndexedVdwPair;
169
170 typedef struct indexed_nbthole_pair
171 {
172         Index ind1;             //  Index for first atom type
173         Index ind2;             //  Index for second atom type
174         Real alphai;           //  Parameter alpha for this pair
175         Real alphaj;           //  Parameter alpha for this pair
176         Real tholeij;          //  Parameter thole for this pair
177         struct indexed_nbthole_pair *right;  //  Right child
178         struct indexed_nbthole_pair *left;   //  Left child
179 } IndexedNbtholePair;
180
181 typedef struct nbthole_pair_value
182 {
183         Index ind1;             //  Index for first atom type
184         Index ind2;             //  Index for second atom type
185         Real alphai;           //  Parameter alpha for this pair
186         Real alphaj;           //  Parameter alpha for this pair
187         Real tholeij;          //  Parameter thole for this pair
188 } NbtholePairValue;
189
190 //  IndexedTablePair is used to form a binary search tree that is
191 //  indexed by table_type index.  This is the tree that will be
192 //  used to search during the actual simulation
193
194 typedef struct indexed_table_pair
195 {
196         Index ind1;             //  Index for first atom type
197         Index ind2;             //  Index for second atom type
198         int K;                  //  number of the table type for this pair
199         struct indexed_table_pair *right;        //  Right child
200         struct indexed_table_pair *left;         //  Left child
201 } IndexedTablePair;
202
203 //  Structures that are defined in Parameters.C
204 struct bond_params;
205 struct angle_params;
206 struct improper_params;
207 struct dihedral_params;
208 struct crossterm_params;
209 struct vdw_params;
210 struct vdw_pair_params;
211 struct nbthole_pair_params;
212 struct table_pair_params;
213
214 class Parameters
215 {
216 private:
217         void initialize();                      //  zeros out pointers
218
219         char *atomTypeNames;                    //  Names of atom types
220         Bool AllFilesRead;                      //  Flag TRUE imples that all
221                                                 //  of the parameter files
222                                                 //  have been read in and
223                                                 //  the arrays have been
224                                                 //  created.
225 //****** BEGIN CHARMM/XPLOR type changes
226         int paramType;                          //  Type (format) of parameter-file
227 //****** END CHARMM/XPLOR type changes
228   bool cosAngles; // True if some angles may be cos-based
229         struct bond_params *bondp;              //  Binary tree of bond params
230         struct angle_params *anglep;            //  Binary tree of angle params
231         struct improper_params *improperp;      //  Linked list of improper par.
232         struct dihedral_params *dihedralp;      //  Linked list of dihedral par.
233         struct crossterm_params *crosstermp;    //  Linked list of cross-term par.
234         struct vdw_params *vdwp;                //  Binary tree of vdw params
235         struct vdw_pair_params *vdw_pairp;      //  Binary tree of vdw pairs
236         struct nbthole_pair_params *nbthole_pairp;      //  Binary tree of nbthole pairs
237         struct table_pair_params *table_pairp;  //  Binary tree of table pairs
238 public:
239         BondValue *bond_array;                  //  Array of bond params
240         AngleValue *angle_array;                //  Array of angle params
241         DihedralValue *dihedral_array;          //  Array of dihedral params
242         ImproperValue *improper_array;          //  Array of improper params
243         CrosstermValue *crossterm_array;        //  Array of crossterm params
244         // JLai
245         GromacsPairValue *gromacsPair_array;    //  Array of gromacsPair params
246         // End of JLai
247         VdwValue *vdw_array;                    //  Array of vdw params
248         NbtholePairValue *nbthole_array;        //  Array of nbthole params
249         int numenerentries;                     //  Number of entries for enertable
250         int rowsize;
251         int columnsize;
252         BigReal* table_ener;                    //  Table for tabulated energies
253         IndexedVdwPair *vdw_pair_tree;          //  Tree of vdw pair params
254         IndexedNbtholePair *nbthole_pair_tree;      //  Tree of nbthole pair params
255         IndexedTablePair *tab_pair_tree;                //  Tree of vdw pair params
256   int tablenumtypes;
257         int NumBondParams;                      //  Number of bond parameters
258         int NumAngleParams;                     //  Number of angle parameters
259   int NumCosAngles;       // Number of cosine-based angles
260         int NumDihedralParams;                  //  Number of dihedral params
261         int NumImproperParams;                  //  Number of improper params
262         int NumCrosstermParams;                 //  Number of cross-term params
263         // JLai
264         int NumGromacsPairParams;               //  Number of gromacsPair params
265         // End of JLai
266         int NumVdwParams;                       //  Number of vdw parameters
267         int NumTableParams;                     //  Number of table parameters
268   int NumVdwParamsAssigned;               //  Number actually assigned
269         int NumVdwPairParams;                   //  Number of vdw_pair params
270         int NumNbtholePairParams;                   //  Number of nbthole_pair params
271         int NumTablePairParams;                 //  Number of vdw_pair params
272 private:
273         ResizeArray<char *> error_msgs;         //  Avoids repeating warnings
274
275         int *maxDihedralMults;                  //  Max multiplicity for
276                                                 //  dihedral bonds
277         int *maxImproperMults;                  //  Max multiplicity for
278                                                 //  improper bonds
279
280         void skip_stream_read(char *, FILE *);  // skip part of stream file
281
282         void add_bond_param(char *);            //  Add a bond parameter
283         struct bond_params *add_to_bond_tree(struct bond_params * , 
284                                      struct bond_params *);
285
286         void add_angle_param(char *);           //  Add an angle parameter
287         struct angle_params *add_to_angle_tree(struct angle_params * , 
288                                      struct angle_params *);
289
290         void add_dihedral_param(char *, FILE *); //  Add a dihedral parameter
291         void add_to_dihedral_list(struct dihedral_params *);
292         void add_to_charmm_dihedral_list(struct dihedral_params *);
293
294         void add_improper_param(char *, FILE *); //  Add an improper parameter
295         void add_to_improper_list(struct improper_params *);
296
297         void add_crossterm_param(char *, FILE *); //  Add an cross-term parameter
298         void add_to_crossterm_list(struct crossterm_params *);
299
300         void add_vdw_param(char *);             //  Add a vdw parameter
301         struct vdw_params *add_to_vdw_tree(struct vdw_params *, 
302                                      struct vdw_params *);
303
304         void add_vdw_pair_param(char *);        //  Add a vdw pair parameter
305         void add_nbthole_pair_param(char *);        //  Add a nbthole pair parameter        
306         void add_table_pair_param(char *);      //  Add a table pair parameter
307         void add_to_vdw_pair_list(struct vdw_pair_params *);
308         void add_to_nbthole_pair_list(struct nbthole_pair_params *);
309         void add_to_table_pair_list(struct table_pair_params *);
310
311         void add_hb_pair_param(char *); //  Add a hydrogen bond pair parameter
312
313         //  All of the traverse routines are used for debugging purposes
314         //  to print out the appropriate list of parameters
315         void traverse_vdw_pair_params(struct vdw_pair_params *);
316         void traverse_nbthole_pair_params(struct nbthole_pair_params *);
317         void traverse_vdw_params(struct vdw_params *);
318         void traverse_dihedral_params(struct dihedral_params *);
319         void traverse_improper_params(struct improper_params *);
320         void traverse_crossterm_params(struct crossterm_params *);
321         void traverse_angle_params(struct angle_params *);
322         void traverse_bond_params(struct bond_params *);
323
324         //  The index_* routines are used to index each of 
325         //  the parameters and build the arrays that will be used
326         //  for constant time access
327         Index index_bonds(struct bond_params *, Index);
328         Index index_angles(struct angle_params *, Index);
329         Index index_vdw(struct vdw_params *, Index);
330         void index_dihedrals();
331         void index_impropers();
332         void index_crossterms();
333         
334         void convert_vdw_pairs();
335         void convert_nbthole_pairs();
336         void convert_table_pairs();
337         IndexedVdwPair *add_to_indexed_vdw_pairs(IndexedVdwPair *, IndexedVdwPair *);
338         IndexedNbtholePair *add_to_indexed_nbthole_pairs(IndexedNbtholePair *, IndexedNbtholePair *);
339         IndexedTablePair *add_to_indexed_table_pairs(IndexedTablePair *, IndexedTablePair *);
340         
341         int vdw_pair_to_arrays(int *, int *, Real *, Real *, Real *, Real *, 
342                                int, IndexedVdwPair *);
343
344         int nbthole_pair_to_arrays(int *, int *, Real *, Real *, Real *, int, IndexedNbtholePair *);
345         
346         int table_pair_to_arrays(int *, int *, int *, int, IndexedTablePair *);
347
348         //  The free_* routines are used by the destructor to deallocate
349         //  memory
350         void free_bond_tree(struct bond_params *);
351         void free_angle_tree(struct angle_params *);
352         void free_dihedral_list(struct dihedral_params *);
353         void free_improper_list(struct improper_params *);
354         void free_crossterm_list(struct crossterm_params *);
355         void free_vdw_tree(struct vdw_params *);
356         void free_vdw_pair_tree(IndexedVdwPair *);
357         void free_nbthole_pair_tree(IndexedNbtholePair *);
358         void free_table_pair_tree(IndexedTablePair *);
359         void free_vdw_pair_list();
360         void free_nbthole_pair_list(); 
361
362   BigReal interp_lin(BigReal, BigReal, BigReal, BigReal, BigReal); // Perform a linear interpolation for energy table 
363
364         /* does the actual initialization, once the variables have all
365            been given default values.  See Parameters() below */
366         void read_parm(const GromacsTopFile *, Bool min);
367 public:
368         //****** BEGIN CHARMM/XPLOR type changes
369         //// added SimParameters to argument list
370         Parameters();
371         Parameters(SimParameters *, StringList *f);
372         //****** END CHARMM/XPLOR type changes
373         
374         Parameters(Ambertoppar *, BigReal);
375         void read_parm(Ambertoppar *, BigReal);
376
377         /* initializes this to hold the set of parameters in the
378            GROMACS topology file <gf>.  If the parameter <min> is on,
379            this assumes that we are going to do minimization and
380            therefore can't have atoms with zero VDW - it will add a
381            small repulsive term to these. */
382         Parameters(const GromacsTopFile *gf, Bool min);
383         
384         ~Parameters();                          //  Destructor
385
386         // return a string for the Nth atom type.  This can only be
387         // called after all the param files have been read and the type
388         // names have been indexed.  The Nth atom type refers to the same
389         // index of the Nth vdw parameter (i.e. there are NumVdwParams names).
390         char *atom_type_name(Index a) {
391           return (atomTypeNames + (a * (MAX_ATOMTYPE_CHARS + 1)));
392         }
393
394         //  Read a parameter file
395         void read_parameter_file(char *);
396
397         //****** BEGIN CHARMM/XPLOR type changes
398         void read_charmm_parameter_file(char *);
399         //****** END CHARMM/XPLOR type changes
400
401         //  Signal the parameter object that all of
402         //  the parameter files have been read in
403         void done_reading_files();
404
405         //  Signal the parameter object that the
406         //  structure file has been read in
407         void done_reading_structure();
408     
409         //  The assign_*_index routines are used to assign
410         //  an index to atoms or bonds.  If an specific atom
411         //  or bond type can't be found, then the program 
412         //  terminates
413     #ifdef MEM_OPT_VERSION
414     void assign_vdw_index(char *, AtomCstInfo *);
415     #else
416         void assign_vdw_index(char *, Atom *);  //  Assign a vdw index to
417     #endif
418                                                 //  an atom
419         void assign_bond_index(char *, char *, Bond *); 
420                                                 //  Assign a bond index
421                                                 //  to a bond
422         void assign_angle_index(char *, char *, char *, Angle *, int);
423                                                 //  Assign an angle index
424                                                 //  to an angle
425         void assign_dihedral_index(char *, char*, char*, char *, Dihedral *, int, int);
426                                                 //  Assign a dihedral index
427                                                 //  to a dihedral
428         void assign_improper_index(char *, char*, char*, char *, Improper *, int);
429                                                 //  Assign an improper index
430                                                 //  to an improper
431         void assign_crossterm_index(char *, char*, char*, char *, char *, char*, char*, char *, Crossterm *);
432
433         //  send_parameters is used by the master process to
434         //  communicate the paramters to all the other processors
435         void send_Parameters(MOStream *);
436
437         //  receive_parameters is used by all the child processes
438         //  to receive the parameters from the master process
439         void receive_Parameters(MIStream *);
440
441         //  The get_*_params routines are the routines that really
442         //  do all the work for this object.  Given an index, they
443         //  access the parameters and return the relevant information
444         void get_bond_params(Real *k, Real *x0, Index index)
445         {
446                 *k = bond_array[index].k;
447                 *x0 = bond_array[index].x0;
448         }
449
450         void get_angle_params(Real *k, Real *theta0, Real *k_ub, Real *r_ub,
451                               Index index)
452         {
453                 *k = angle_array[index].k;
454                 *theta0 = angle_array[index].theta0;
455                 *k_ub = angle_array[index].k_ub;
456                 *r_ub = angle_array[index].r_ub;
457         }
458
459         int get_improper_multiplicity(Index index)
460         {
461                 return(improper_array[index].multiplicity);
462         }
463
464         int get_dihedral_multiplicity(Index index)
465         {
466                 return(dihedral_array[index].multiplicity);
467         }
468
469         void get_improper_params(Real *k, int *n, Real *delta, 
470                                  Index index, int mult)
471         {
472                 if ( (mult<0) || (mult>MAX_MULTIPLICITY) )
473                 {
474                         NAMD_die("Bad mult index in Parameters::get_improper_params");
475                 }
476
477                 *k = improper_array[index].values[mult].k;
478                 *n = improper_array[index].values[mult].n;
479                 *delta = improper_array[index].values[mult].delta;
480         }
481
482         void get_dihedral_params(Real *k, int *n, Real *delta, 
483                                  Index index, int mult)
484         {
485                 if ( (mult<0) || (mult>MAX_MULTIPLICITY) )
486                 {
487                         NAMD_die("Bad mult index in Parameters::get_dihedral_params");
488                 }
489
490                 *k = dihedral_array[index].values[mult].k;
491                 *n = dihedral_array[index].values[mult].n;
492                 *delta = dihedral_array[index].values[mult].delta;
493         }
494
495         void get_vdw_params(Real *sigma, Real *epsilon, Real *sigma14, 
496                             Real *epsilon14, Index index)
497   {
498     if ( vdw_array ) {
499       *sigma = vdw_array[index].sigma;
500       *epsilon = vdw_array[index].epsilon;
501       *sigma14 = vdw_array[index].sigma14;
502       *epsilon14 = vdw_array[index].epsilon14;
503     } else {
504       // sigma and epsilon from A and B for each vdw type's interaction with itself
505       Real A; Real B; Real A14; Real B14;
506       if (get_vdw_pair_params(index, index, &A, &B, &A14, &B14) ) {
507         if (A && B) {
508           *sigma = sqrt(cbrt(A)) / sqrt(cbrt(B));
509           *epsilon = B*B / (4*A);
510         }
511         else {
512           *sigma = 0; *epsilon = 0;
513         }
514         if (A14 && B14) {
515           *sigma14 = sqrt(cbrt(A14)) / sqrt(cbrt(B14));
516           *epsilon14 = B14*B14 / (4*A14);
517         }
518         else {
519           *sigma14 = 0; *epsilon14 = 0;
520         }
521       }
522       else {NAMD_die ("Function get_vdw_params failed to derive Lennard-Jones sigma and epsilon from A and B values\n");}
523     }
524   }
525
526         int get_vdw_pair_params(Index ind1, Index ind2, Real *, Real *, Real *, Real *);
527                                                 //  Find a vwd_pair parameter
528
529         int get_num_vdw_params(void) { return NumVdwParamsAssigned; }
530
531   int get_table_pair_params(Index, Index, int*);
532
533         //  The print_*_params are provided for debugging purposes
534         void print_bond_params();               //  Print bonds params
535         void print_angle_params();              //  Print angle params
536         void print_dihedral_params();           //  Print dihedral params
537         void print_improper_params();           //  Print improper params
538         void print_crossterm_params();          //  Print cross-term params
539         void print_vdw_params();                //  Print vdw params
540         void print_vdw_pair_params();           //  Print vdw_pair params
541         void print_nbthole_pair_params();           //  Print nbthole_pair params
542         void print_param_summary();             //  Print a summary of params
543         void read_ener_table(SimParameters*); // Read an energy table file
544   int get_int_table_type(char*); // Return the integer type for a named table interaction
545
546   int read_energy_type(FILE*, const int, BigReal*, const float, const float);
547   int read_energy_type_cubspline(FILE*, const int, BigReal*, const float, const float);
548   int read_energy_type_bothcubspline(FILE*, const int, BigReal*, const float, const float);
549 };
550
551 #endif
552