Fix REST2 config parameter name lookup
[namd.git] / src / Molecule.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    This class is used to store all of the structural       
9    information for a simulation.  It reads in this information 
10    from a .psf file, cross checks and obtains some information
11    from the Parameters object that is passed in, and then     
12    stores all this information for later use. 
13 */
14
15
16 #ifndef MOLECULE_H
17
18 #define MOLECULE_H
19
20 #include "parm.h"
21
22 #include "common.h"
23 #include "NamdTypes.h"
24 #include "structures.h"
25 #include "ConfigList.h"
26 #include "Vector.h"
27 #include "UniqueSet.h"
28 #include "Hydrogen.h"
29 #include "GromacsTopFile.h"
30 /* BEGIN gf */
31 #include "GridForceGrid.h"
32 #include "Tensor.h"
33 /* END gf */
34
35 // Go -- JLai
36 #define MAX_GO_CHAINS          10
37 #define MAX_RESTRICTIONS       10
38 // End of Go
39
40 #include "molfile_plugin.h"
41
42 #include <vector>
43
44 class SimParameters;
45 class Parameters;
46 class ConfigList;
47 class PDB;
48 class MIStream;
49 class MOStream;
50
51 class ExclElem;
52 class BondElem;
53 class AngleElem;
54 class DihedralElem;
55 class ImproperElem;
56 class TholeElem;  // Drude model
57 class AnisoElem;  // Drude model
58 class CrosstermElem;
59 // JLai
60 class GromacsPairElem;
61 // End of JLai
62 class ResidueLookupElem;
63
64 struct OutputAtomRecord;
65
66 template<class Type> class ObjectArena;
67
68 class ExclusionCheck {
69 public:
70   int32 min,max;
71   char *flags;
72
73   ExclusionCheck(){
74       min=0;
75       max=-1;
76       flags = NULL;
77   }
78 private:
79   ExclusionCheck(const ExclusionCheck& chk){
80   }
81   ExclusionCheck &operator=(const ExclusionCheck& chk){
82     return *this;
83   }
84 };
85 #define EXCHCK_FULL 1
86 #define EXCHCK_MOD 2
87
88 class ResidueLookupElem
89 {
90 public:
91   char mySegid[11];
92   ResidueLookupElem *next;      // stored as a linked list
93   int firstResid;               // valid resid is >= firstResid
94   int lastResid;                // valid resid is <= lastResid
95   ResizeArray<int> atomIndex;   // 0-based index for first atom in residue
96
97   ResidueLookupElem(void) { next = 0; firstResid = -1; lastResid = -1; }
98   ~ResidueLookupElem(void) { delete next; }
99   int lookup(const char *segid, int resid, int *begin, int *end) const;
100   ResidueLookupElem* append(const char *segid, int resid, int aid);
101 };
102
103 // Ported by JLai -- JE
104 typedef struct go_val
105 {
106   Real epsilon;      // Epsilon
107   int exp_a;         // First exponent for attractive L-J term
108   int exp_b;         // Second exponent for attractive L-J term
109   int exp_rep;       // Exponent for repulsive L-J term
110   Real sigmaRep;     // Sigma for repulsive term
111   Real epsilonRep;   // Epsilon for replusive term
112   Real cutoff;       // Cutoff distance for Go calculation
113   int  restrictions[MAX_RESTRICTIONS];  //  List of residue ID differences to be excluded from Go calculation
114 } GoValue;
115
116 typedef struct go_pair
117 {
118   int goIndxA;       // indexA
119   int goIndxB;       // indexB
120   double A;          // double A for the LJ pair
121   double B;          // double B for the LJ pair
122 } GoPair;
123 // End of port - JL
124
125 #define QMLSSMODEDIST 1
126 #define QMLSSMODECOM 2
127
128 #define QMFormatORCA 1
129 #define QMFormatMOPAC 2
130 #define QMFormatUSR 3
131
132 #define QMCHRGNONE 0
133 #define QMCHRGMULLIKEN 1
134 #define QMCHRGCHELPG 2
135
136 #define QMSCHEMECS 1
137 #define QMSCHEMERCD 2
138 #define QMSCHEMEZ1 3
139 #define QMSCHEMEZ2 4
140 #define QMSCHEMEZ3 5
141
142 //only used for compressing the molecule information
143 typedef struct seg_resid
144 {
145     char segname[11];
146     int resid;
147 }AtomSegResInfo;
148
149 // List maintaining the global atom indicies sorted by helix groups.
150 class Molecule
151 {
152  private:
153 typedef struct constraint_params
154 {
155    Real k;    //  Force constant
156    Vector refPos; //  Reference position for restraint
157 } ConstraintParams;
158
159
160
161 /* BEGIN gf */
162 typedef struct gridfrc_params
163 {
164     Real k;     // force multiplier
165     Charge q;   // charge
166 } GridforceParams;
167 /* END gf */
168
169
170 typedef struct stir_params
171 {
172   Real startTheta;
173   Vector refPos;  //  Reference position for stirring
174 } StirParams;
175
176 typedef struct movdrag_params
177 {
178   Vector v;            //  Linear velocity (A/step)
179 } MovDragParams;
180
181
182 typedef struct rotdrag_params
183 {
184    Real v;              //  Angular velocity (deg/step)
185    Vector a;            //  Rotation axis
186    Vector p;            //  Rotation pivot point
187 } RotDragParams;
188
189 typedef struct constorque_params
190 {
191    Real v;              //  "Torque" value (Kcal/(mol*A^2))
192    Vector a;            //  Rotation axis
193    Vector p;            //  Rotation pivot point
194 } ConsTorqueParams;
195
196 #ifdef MEM_OPT_VERSION
197 //indicate a range of atoms from aid1 to aid2
198 struct AtomSet{
199     AtomID aid1;
200     AtomID aid2;
201
202     int operator < (const AtomSet &a) const{
203                 return aid1 < a.aid1;
204         }
205 };
206 typedef ResizeArray<AtomSet> AtomSetList;
207
208   void load_atom_set(StringList *setfile, const char *setname,
209         int *numAtomsInSet, Molecule::AtomSetList **atomsSet) const;
210
211 #endif
212
213 friend class ExclElem;
214 friend class BondElem;
215 friend class AngleElem;
216 friend class DihedralElem;
217 friend class ImproperElem;
218 friend class TholeElem;  // Drude model
219 friend class AnisoElem;  // Drude model
220 friend class CrosstermElem;
221 // JLai
222 friend class GromacsPairElem;
223 // End of JLai
224 friend class WorkDistrib;
225
226 private:
227
228 #ifndef MEM_OPT_VERSION
229         Atom *atoms;    //  Array of atom structures
230         ObjectArena<char> *nameArena;
231         AtomNameInfo *atomNames;//  Array of atom name info.  Only maintained on node 0 for VMD interface
232         //replaced by atom signatures
233         Bond *bonds;    //  Array of bond structures
234         Angle *angles;    //  Array of angle structures
235         Dihedral *dihedrals;  //  Array of dihedral structures
236         Improper *impropers;  //  Array of improper structures                          
237         Crossterm *crossterms;  //  Array of cross-term structures
238         GromacsPair *gromacsPair; //  Array of gromacs-pair structures
239
240         //These will be replaced by exclusion signatures
241         Exclusion *exclusions;  //  Array of exclusion structures
242         UniqueSet<Exclusion> exclusionSet;  //  Used for building
243
244         int32 *cluster;   //  first atom of connected cluster
245
246         ObjectArena<int32> *tmpArena;
247         int32 **bondsWithAtom;  //  List of bonds involving each atom
248         ObjectArena<int32> *arena;
249
250         //function is replaced by atom signatures
251         int32 **bondsByAtom;  //  List of bonds owned by each atom
252         int32 **anglesByAtom;     //  List of angles owned by each atom
253         int32 **dihedralsByAtom;  //  List of dihedrals owned by each atom
254         int32 **impropersByAtom;  //  List of impropers owned by each atom
255         int32 **crosstermsByAtom;  //  List of crossterms owned by each atom
256
257         int32 **exclusionsByAtom; //  List of exclusions owned by each atom
258         int32 **fullExclusionsByAtom; //  List of atoms excluded for each atom
259         int32 **modExclusionsByAtom; //  List of atoms modified for each atom
260 // JLai
261         int32 **gromacsPairByAtom; // gromacsPair exclusion list which by definition should not have any exclusions (still not sure if it should be empty or zeroed out)
262 // End of JLai
263
264         ObjectArena<char> *exclArena;
265         ExclusionCheck *all_exclusions;
266
267         // DRUDE
268         int32 **tholesByAtom;  // list of Thole correction terms owned by each atom
269         int32 **anisosByAtom;  // list of anisotropic terms owned by each atom
270         // DRUDE
271
272 #else
273         //Indexing to constant pools to save space
274         AtomCstInfo *atoms;
275         Index *eachAtomMass; //after initialization, this could be freed (possibly)
276         Index *eachAtomCharge; //after initialization, this could be freed (possibly)
277         AtomNameIdx *atomNames;
278         ObjectArena<char> *nameArena; //the space for names to be allocated  
279
280         //A generally true assumption: all atoms are arranged in the order of clusters. In other words,
281         //all atoms for a cluster must appear before/after any atoms in other clusters
282         //The first atom in the cluster (which has the lowest atom id) stores the cluster size
283         //other atoms in the cluster stores -1
284         int32 *clusterSigs;
285         int32 numClusters;
286         
287         AtomSigID *eachAtomSig;
288         ExclSigID *eachAtomExclSig;
289
290     AtomSetList *fixedAtomsSet;    
291     AtomSetList *constrainedAtomsSet;    
292 #endif
293  
294   ResidueLookupElem *resLookup; // find residues by name
295
296   AtomSegResInfo *atomSegResids; //only used for generating compressed molecule info
297
298   Bond *donors;         //  Array of hydrogen bond donor structures
299   Bond *acceptors;  //  Array of hydrogen bond acceptor
300
301   // DRUDE
302   DrudeConst *drudeConsts;  // supplement Atom data (length of Atom array)
303   Thole *tholes;            // Thole (screened Coulomb) interactions
304   Aniso *anisos;            // anisotropic terms
305   Lphost *lphosts;          // lone pair hosts
306   int32 *lphostIndexes;     // index for each LP into lphosts array
307   // DRUDE
308
309   int32 *consIndexes; //  Constraint indexes for each atom
310   ConstraintParams *consParams;
311
312 /* BEGIN gf */
313   int32 **gridfrcIndexes;
314   GridforceParams **gridfrcParams;
315   GridforceGrid **gridfrcGrid;
316 /* END gf */
317
318         //  Parameters for each atom constrained
319   int32 *stirIndexes; //Stirring indexes for each atoms
320   StirParams *stirParams;
321                           // Paramters for each atoms stirred
322   int32 *movDragIndexes;  //  Moving drag indexes for each atom
323   MovDragParams *movDragParams;
324                                 //  Parameters for each atom moving-dragged
325   int32 *rotDragIndexes;  //  Rotating drag indexes for each atom
326   RotDragParams *rotDragParams;
327                                 //  Parameters for each atom rotation-dragged
328
329   Real *langevinParams;   //  b values for langevin dynamics
330   int32 *fixedAtomFlags;  //  1 for fixed, -1 for fixed group, else 0
331   int32 *exPressureAtomFlags; // 1 for excluded, -1 for excluded group.
332
333   //In the memory optimized version: it will be NULL if the general
334   //true assumption mentioned above holds. If not, its size is numClusters.
335   //In the ordinary version: its size is numAtoms, and indicates the size
336   //of connected cluster or 0.
337   int32 *clusterSize;                            
338
339         Real *rigidBondLengths;  //  if H, length to parent or 0. or
340         //  if not H, length between children or 0.
341 //fepb
342         unsigned char *fepAtomFlags; 
343 //fepe
344
345 //soluteScaling
346         unsigned char *ssAtomFlags;
347         unsigned int num_ss;
348 //soluteScaling
349
350   //occupancy and bfactor data from plugin-based IO implementation of loading structures
351   float *occupancy;
352   float *bfactor;
353
354
355   // added during the trasition from 1x to 2
356   SimParameters *simParams;
357   Parameters *params;
358
359 private:
360         void initialize(SimParameters *, Parameters *param);
361         // Sets most data to zero
362
363   //LCPO
364   int *lcpoParamType;
365
366 #ifndef MEM_OPT_VERSION
367         void read_psf_file(char *, Parameters *);
368         //  Read in a .psf file given
369         //  the filename and the parameter
370         //  object to use          
371   
372   void read_atoms(FILE *, Parameters *);
373         //  Read in atom info from .psf
374   void read_bonds(FILE *, Parameters *);
375         //  Read in bond info from .psf
376   void read_angles(FILE *, Parameters *);
377         //  Read in angle info from .psf
378   void read_dihedrals(FILE *, Parameters *);
379         //  Read in dihedral info from .psf
380   void read_impropers(FILE *, Parameters *);
381         //  Read in improper info from .psf
382   void read_crossterms(FILE *, Parameters *);
383         //  Read in cross-term info from .psf
384   void read_donors(FILE *);
385         //  Read in hydrogen bond donors from .psf
386   void read_acceptors(FILE *);
387         //  Read in hydrogen bond acceptors from .psf
388   void read_exclusions(FILE *);
389         //  Read in exclusion info from .psf
390   // JLai August 16th, 2012 Modifications
391   void read_exclusions(int*, int*, int);
392   /* Read in exclusion array and sort entries */
393   static bool goPairCompare (GoPair, GoPair);
394   // End of JLai August 16th, 2012 Modifications
395
396
397   // DRUDE: PSF reading
398   void read_lphosts(FILE *);
399         //  Read in lone pair hosts from Drude PSF
400   void read_anisos(FILE *);
401         //  Read in anisotropic terms from Drude PSF
402   // DRUDE
403
404   //LCPO
405   //input type is Charmm/Amber/other
406   //0 - Charmm/Xplor
407   //1 - Amber TODO
408   //2 - Plugin TODO
409   //3 - Gromacs TODO
410   void assignLCPOTypes(int inputType);
411
412   //pluginIO-based loading atoms' structure
413   void plgLoadAtomBasics(molfile_atom_t *atomarray);
414   void plgLoadBonds(int *from, int *to); //atom index is 1-based in the parameters
415   void plgLoadAngles(int *plgAngles);
416   void plgLoadDihedrals(int *plgDihedrals);
417   void plgLoadImpropers(int *plgImpropers);
418   void plgLoadCrossterms(int *plgCterms);
419
420         //  List of all exclusions, including
421         //  explicit exclusions and those calculated
422         //  from the bonded structure based on the
423         //  exclusion policy.  Also includes 1-4's.
424   void build_lists_by_atom();
425         //  Build the list of structures by atom
426
427   void build12excl(void);
428   void build13excl(void);
429   void build14excl(int);
430
431   // DRUDE: extend exclusions for Drude and LP
432   void build_inherited_excl(int);
433   // DRUDE
434   void stripFepExcl(void);
435
436   void build_exclusions();
437   // analyze the atoms, and determine which are oxygen, hb donors, etc.
438   // this is called after a molecule is sent our (or received in)
439   void build_atom_status(void);
440
441 #else
442   //the method to load the signatures of atoms etc. (i.e. reading the file in 
443   //text fomrat of the compressed psf file)
444   void read_mol_signatures(char *fname, Parameters *params, ConfigList *cfgList=0);     
445   void load_one_inputatom(int aid, OutputAtomRecord *one, InputAtom *fAtom);
446   void build_excl_check_signatures();  
447 #endif
448
449         void read_parm(const GromacsTopFile *);  
450           
451 public:
452
453   // for solute scaling
454   int ss_num_vdw_params;
455   int *ss_vdw_type;
456   int *ss_index;
457
458   // DRUDE
459   int is_drude_psf;      // flag for reading Drude PSF
460   int is_lonepairs_psf;  // flag for reading lone pairs from PSF
461   // DRUDE
462
463   // data for TIP4P
464   Real r_om;
465   Real r_ohc;
466
467   // data for tail corrections
468   BigReal tail_corr_ener;
469   BigReal tail_corr_virial;
470   BigReal tail_corr_dUdl_1, tail_corr_dUdl_2;
471   BigReal tail_corr_virial_1, tail_corr_virial_2;
472   void compute_LJcorrection();
473   BigReal getEnergyTailCorr(const BigReal);
474   BigReal getVirialTailCorr(const BigReal);
475
476   int const * getLcpoParamType() {
477     return lcpoParamType;
478   }
479
480   BigReal GetAtomAlpha(int i) const {
481     return drudeConsts[i].alpha;
482   }
483
484 #ifdef MEM_OPT_VERSION
485   AtomCstInfo *getAtoms() const { return atoms; }
486   AtomNameIdx *getAtomNames() const { return atomNames; }
487 #else
488   Atom *getAtoms () const { return atoms; }
489   AtomNameInfo *getAtomNames() const { return atomNames; }
490 #endif
491
492   AtomSegResInfo *getAtomSegResInfo() const { return atomSegResids; }
493   
494   // return number of fixed atoms, guarded by SimParameters
495   int num_fixed_atoms() const {
496     // local variables prefixed by s_
497     int s_NumFixedAtoms = (simParams->fixedAtomsOn ? numFixedAtoms : 0);
498     return s_NumFixedAtoms;  // value is "turned on" SimParameters
499   }
500
501   int num_fixed_groups() const {
502     // local variables prefixed by s_
503     int s_NumFixedAtoms = num_fixed_atoms();
504     int s_NumFixedGroups = (s_NumFixedAtoms ? numFixedGroups : 0);
505     return s_NumFixedGroups;  // value is "turned on" SimParameters
506   }
507
508   int64_t num_group_deg_freedom() const {
509     // local variables prefixed by s_
510     int64_t s_NumGroupDegFreedom = 3 * (int64_t) numHydrogenGroups;
511     int64_t s_NumFixedAtoms = num_fixed_atoms();
512     int s_NumFixedGroups = num_fixed_groups();
513     if (s_NumFixedGroups) s_NumGroupDegFreedom -= 3 * s_NumFixedGroups;
514     if ( ! (s_NumFixedAtoms || numConstraints
515           || simParams->comMove || simParams->langevinOn) ) {
516       s_NumGroupDegFreedom -= 3;
517     }
518     return s_NumGroupDegFreedom;
519   }
520
521   int64_t num_deg_freedom(int isInitialReport = 0) const {
522     // local variables prefixed by s_
523     int64_t s_NumDegFreedom = 3 * (int64_t) numAtoms;
524     int64_t s_NumFixedAtoms = num_fixed_atoms();
525     if (s_NumFixedAtoms) s_NumDegFreedom -= 3 * s_NumFixedAtoms;
526     if (numLonepairs) s_NumDegFreedom -= 3 * numLonepairs;
527     if ( ! (s_NumFixedAtoms || numConstraints
528           || simParams->comMove || simParams->langevinOn) ) {
529       s_NumDegFreedom -= 3;
530     }
531     if ( ! isInitialReport && simParams->pairInteractionOn) {
532       //
533       // DJH: a kludge?  We want to initially report system degrees of freedom
534       //
535       // this doesn't attempt to deal with fixed atoms or constraints
536       s_NumDegFreedom = 3 * numFepInitial;
537     }
538     int s_NumFixedRigidBonds = 
539       (simParams->fixedAtomsOn ? numFixedRigidBonds : 0);
540     if (simParams->watmodel == WAT_TIP4) {
541       // numLonepairs is subtracted here because all lonepairs have a rigid bond
542       // to oxygen, but all of the LP degrees of freedom are dealt with above
543       s_NumDegFreedom -= (numRigidBonds - s_NumFixedRigidBonds - numLonepairs);
544     }
545     else {
546       // Non-polarized systems don't have LPs.
547       // For Drude model, bonds that attach LPs are not counted as rigid;
548       // LPs have already been subtracted from degrees of freedom.
549       s_NumDegFreedom -= (numRigidBonds - s_NumFixedRigidBonds);
550     }
551     return s_NumDegFreedom;
552   }
553
554   int numAtoms;   //  Number of atoms                   
555
556   int numRealBonds;   //  Number of bonds for exclusion determination
557   int numBonds;   //  Number of bonds calculated, including extras
558   int numAngles;    //  Number of angles
559   int numDihedrals; //  Number of dihedrals
560   int suspiciousAlchBonds;    //  angles dropped due to crossing FEP partitions
561   int alchDroppedAngles;    //  angles dropped due to crossing FEP partitions
562   int alchDroppedDihedrals; //  dihedrals dropped due to crossing FEP partitions
563   int alchDroppedImpropers; //  impropers dropped due to crossing FEP partitions
564   int numImpropers; //  Number of impropers
565   int numCrossterms; //  Number of cross-terms
566   int numDonors;          //  Number of hydrogen bond donors
567   int numAcceptors; //  Number of hydrogen bond acceptors
568   int numExclusions;  //  Number of exclusions
569   int64 numTotalExclusions; //  Real Total Number of Exclusions // hack
570
571   // DRUDE
572   int numLonepairs; // Number of lone pairs
573   int numDrudeAtoms;  // Number of Drude particles
574   int numTholes;  // Number of Thole terms
575   int numAnisos;  // Number of anisotropic terms
576   int numLphosts;  ///< Number of lone pair host records in PSF
577   // DRUDE
578   
579   int numConstraints; //  Number of atoms constrained
580 /* BEGIN gf */
581   int numGridforceGrids;//  Number of gridforce grids
582   int *numGridforces;   //  Number of atoms in gridforce file (array, one per grid)
583 /* END gf */
584   int numMovDrag;         //  Number of atoms moving-dragged
585   int numRotDrag;         //  Number of atoms rotating-dragged
586   int numConsTorque;  //  Number of atoms "constant"-torqued
587   int numFixedAtoms;  //  Number of fixed atoms
588   int numStirredAtoms;  //  Number of stirred atoms
589   int numExPressureAtoms; //  Number of atoms excluded from pressure
590   int numHydrogenGroups;  //  Number of hydrogen groups
591   int maxHydrogenGroupSize;  //  Max atoms per hydrogen group
592   int numMigrationGroups;  //  Number of migration groups
593   int maxMigrationGroupSize;  //  Max atoms per migration group
594   int numFixedGroups; //  Number of totally fixed hydrogen groups
595   int numRigidBonds;  //  Number of rigid bonds
596   int numFixedRigidBonds; //  Number of rigid bonds between fixed atoms
597 //fepb
598   int numFepInitial;  // no. of fep atoms with initial flag
599   int numFepFinal;  // no. of fep atoms with final flag
600 //fepe
601
602   int numConsForce; //  Number of atoms that have constant force applied
603   int32 *consForceIndexes;//  Constant force indexes for each atom
604   Vector *consForce;  //  Constant force array
605
606   int32 *consTorqueIndexes; //  "Constant" torque indexes for each atom
607   ConsTorqueParams *consTorqueParams;
608                                 //  Parameters for each atom "constant"-torqued
609
610   // The following are needed for error checking because we
611   // eliminate bonds, etc. which involve only fixed atoms
612   int numCalcBonds; //  Number of bonds requiring calculation
613   int numCalcAngles;  //  Number of angles requiring calculation
614   int numCalcDihedrals; //  Number of dihedrals requiring calculation
615   int numCalcImpropers; //  Number of impropers requiring calculation
616   int numCalcCrossterms; //  Number of cross-terms requiring calculation
617   int64 numCalcExclusions;  //  Number of exclusions requiring calculation
618   int64 numCalcFullExclusions;  //  Number of full exclusions requiring calculation
619
620   // DRUDE
621   int numCalcTholes;  // Number of Thole correction terms requiring calculation
622   int numCalcAnisos;  // Number of anisotropic terms requiring calculation
623   // DRUDE
624
625   //  Number of dihedrals with multiple periodicity
626   int numMultipleDihedrals; 
627   //  Number of impropers with multiple periodicity
628   int numMultipleImpropers; 
629   // indexes of "atoms" sorted by hydrogen groups
630   HydrogenGroup hydrogenGroup;
631
632   // Ported by JLai -- JE - Go
633   int numGoAtoms;         //  Number of atoms subject to Go forces -- ported by JLai/ Original by JE
634   int32 *atomChainTypes;  //  Go chain type for each atom; from 1 to MAX_GO_CHAINS
635   int32 *goSigmaIndices;  //  Indices into goSigmas
636   int32 *goResidIndices;  //  Indices into goSigmas
637   Real  *goSigmas;        //  Sigma values for Go forces L-J type formula
638   bool *goWithinCutoff;   //  Whether the reference atom-atom distance is within the Go cutoff
639   Real  *goCoordinates;   //  Coordinates (x,y,z) for Go atoms in the native structure
640   int *goResids;          //  Residue ID from PDB
641   PDB *goPDB;             //  Pointer to PDB object to use
642   // NAMD-Go2 calculation code
643   int goNumLJPair;        //  Integer storing the total number of explicit pairs (LJ)
644   int *goIndxLJA;         //  Pointer to the array of atom indices for LJ atom A
645   int *goIndxLJB;         //  Pointer to the array of atom indices for LJ atom B
646   double *goSigmaPairA;  //  Pointer to the array of A LJ parameters
647   double *goSigmaPairB;  //  Pointer to the array of B LJ parameters
648   int *pointerToGoBeg;    //  Array of pointers to array
649   int *pointerToGoEnd;    //  Array of pointers to array
650   // Gromacs LJ Pair list calculation code
651   int numPair;            //  Integer storing the total number of explicit pairs (LJ + Gaussian)
652   int numLJPair;          //  Integer storing the number of explicit LJ pairs
653   int numCalcLJPair;         //  Number of explicit LJ pairs requiring calculation
654   int *pointerToLJBeg;       //  Array of pointers to array 
655   int *pointerToLJEnd;       //  Array of pointers to array B
656   int *indxLJA;           //  Pointer to the array of atom indices for LJ atom A
657   int *indxLJB;           //  Pointer to the array of atom indices for LJ atom B
658   Real *pairC6;           //  Pointer to the array of C6 LJ parameters
659   Real *pairC12;          //  Pointer to the array of C12 LJ parameters
660   int *gromacsPair_type;   //  
661   // Gromacs Gauss Pair list calculation code
662   int *pointerToGaussBeg;    //  Array of pointers to array B
663   int *pointerToGaussEnd;    //  Array of pointers to array B
664   int numGaussPair;       //  Integer storing the number of explicit Gaussian pairs  
665   int *indxGaussA;        //  Pointer to the array of atom indices for Gaussian atom A 
666   int *indxGaussB;        //  Pointer to the array of atom indices for Gaussian atom B 
667   Real *gA;               //  Pointer to the array of force constants to the Gaussian potential
668   Real *gMu1;             //  Pointer to the array of mu (shifts Gaussian)
669   Real *giSigma1;          //  Pointer to the array of sigma (controls spread of Gaussian)
670   Real *gMu2;             //  Pointer to the array of mu (shifts Gaussian 2)
671   Real *giSigma2;          //  Pointer to the array of sigma (controls spread of Gaussian 2)
672   Real *gRepulsive;       //  Pointer to the a LJ-12 repulsive parameter that adds to the Gaussian
673
674   // GO ENERGY CALCULATION CODE
675   BigReal energyNative;    // Holds the energy value of the native structure
676   BigReal energyNonnative; // Holds the energy value of the nonnative structure
677   // GO ENERGY CALCULATION CODE
678   // End of port - JL
679
680   
681 private:
682   /////////////////////////
683   // Begin QM
684   /////////////////////////
685     
686   int qmDroppedBonds, qmDroppedAngles, qmDroppedDihedrals;
687   int qmDroppedImpropers, qmDroppedCrossterms;
688   ResizeArray<Bond> qmExtraBonds;
689   
690   Bool qmReplaceAll;              // Indicates if all forces should be replaced.
691   int qmNumQMAtoms;           // Number of QM atoms, from all groups.
692   
693   // Array of abbreviated element names for all atoms.
694   char **qmElementArray;
695   // Array of elements for dummy atoms.
696   char **qmDummyElement;
697   
698   // Number of QM atoms in each QM group
699   int *qmGrpSizes;
700   
701   // QM group for each atom of the system. 0 means MM atom.
702   Real *qmAtomGroup;
703   // Number of different QM regions
704   int qmNumGrps ;
705   
706   // QM Atom charges.
707   Real *qmAtmChrg ;
708   // global indices of all QM atoms.
709   int *qmAtmIndx ;
710   
711   // Number of QM-MM bonds.
712   int qmNumBonds ;
713   // IDs of each group, that is, the value in the qmColumn, per group.
714   Real *qmGrpID ;
715   // The total charge of each QM group.
716   Real *qmGrpChrg;
717   // The multiplicity of QM groups
718   Real *qmGrpMult;
719   // Number of QM-MM bonds in each group.
720   int *qmGrpNumBonds ;
721   // Sequential list of bonds between a MM atom and a QM atom.
722   // (with the bonded atoms ins this order: MM first, QM second).
723   int **qmMMBond ;
724   // List of bond indexes for each QM group.
725   int **qmGrpBonds ;
726   // List of atom indexes for MM atoms in QM-MM bonds, per group.
727   // This is used to check if a point charge is actually an MM atom
728   // that need to be ignored, and will be substituted by a dummy atom (for example).
729   int **qmMMBondedIndx ;
730   // List of indices for MM atoms which will receive fractions of charges from 
731   // the MM atom of a QM-MM bond. This is used in the Charge Shift scheme.
732   int **qmMMChargeTarget;
733   // Number of target MM atoms per QM-MM bond. Targets will receive the partial
734   // charge of the MM atom in a QM-MM bond. (in Charge shift scheme)
735   int *qmMMNumTargs;
736   // List of distances from thr QM atom where the dummy atom will be placed.
737   BigReal *qmDummyBondVal;
738   // Indicate if point charges will be used in QM calculations.
739   Bool qmNoPC;
740   
741   /////////////////////////
742   // These variables are used in case mechanichal embedding is requested (NoPC = TRUE)
743   // AND there are QM-MM bonds in the system.
744   
745   // Indicates the number of items in the arrays for Mechanical embedding with QM-MM bonds.
746   int qmMeNumBonds;
747   // Indicates the indices of MM atoms which participate in QM-MM bonds.
748   int *qmMeMMindx;
749   // Indicates the QM group to which the QM-MM bonds belongs.
750   Real *qmMeQMGrp;
751   /////////////////////////
752   
753   // Indicates frequency of selection of point charges.
754   int qmPCFreq;
755   // Custom PC array;
756   int *qmCustomPCIdxs;
757   // Number of Indexes associated with each QM Group.
758   int *qmCustPCSizes;
759   // Total number of custom PC charges.
760   int qmTotCustPCs;
761   // Number of solvent molecules that will be selected and updated by NAMD, per group.
762   int *qmLSSSize;
763   // Number of atoms per solvent molecule. We need all molecules to have the same size.
764   int qmLSSResidueSize;
765   // IDs of all atoms belonging to QM solvent molecules, from all groups.
766   int *qmLSSIdxs;
767   // MAss of each atom in LSS solvent molecules.
768   Mass *qmLSSMass;
769   // Atom IDs of QM atoms which will be used to select solvent molecules by distance.
770   int *qmLSSRefIDs;
771   // Mass of each QM atom used as reference for solvent selection.
772   Mass *qmLSSRefMass ;
773   // Number of atom IDs/Mass per group.
774   int *qmLSSRefSize;
775   int qmLSSFreq;
776   int qmLSSTotalNumAtms;
777   // This will map each classical solvent atom with a global index of solvent residues,
778   // so we don't depend on segment names and repeated residue indices.
779   std::map<int,int> qmClassicSolv ;
780   
781   /////////////////////////
782   // Conditional SMD
783     
784     void read_qm_csdm_file(std::map<Real,int> &qmGrpIDMap) ;
785     
786     Bool qmcSMD;
787         // Total number of cSMD instances.
788         int cSMDnumInst;
789     // Num instances per QM group
790     int* cSMDindxLen;
791     // Index of Conditional SMD guides per QM group
792     int** cSMDindex;
793     // Atom indices for Origin and Target of cSMD
794     int** cSMDpairs;
795     // Spring constants for cSMD
796     Real* cSMDKs;
797     // Speed of movement of guide particles for cSMD.
798     Real* cSMDVels;
799     // Distance cutoff for guide particles for cSMD.
800     Real* cSMDcoffs;
801   
802   /////////////////////////
803   // end QM
804   /////////////////////////
805   
806 public:
807   // QM
808   void set_qm_replaceAll(Bool newReplaceAll) { qmReplaceAll = newReplaceAll; } ;
809   
810   const Real * get_qmAtomGroup() const {return qmAtomGroup; } ;
811   Real get_qmAtomGroup(int indx) const {return qmAtomGroup[indx]; } ;
812   
813   Real *get_qmAtmChrg() {return qmAtmChrg; } ;
814   const int *get_qmAtmIndx() {return qmAtmIndx; } ;
815   
816   int get_numQMAtoms() {return qmNumQMAtoms; } ;
817   const char *const * get_qmElements() {return qmElementArray;} ;
818   int get_qmNumGrps() const {return qmNumGrps; } ;
819   const int * get_qmGrpSizes() {return qmGrpSizes; } ;
820   Real * get_qmGrpID() { return qmGrpID; } ;
821   Real *get_qmGrpChrg() {return qmGrpChrg; } ;
822   Real *get_qmGrpMult() {return qmGrpMult; } ;
823   
824   int get_qmNumBonds() { return qmNumBonds; } ;
825   const BigReal * get_qmDummyBondVal() { return qmDummyBondVal; } ;
826   const int * get_qmGrpNumBonds() { return qmGrpNumBonds; } ;
827   const int *const * get_qmMMBond() { return qmMMBond; } ;
828   const int *const * get_qmGrpBonds() { return qmGrpBonds; } ;
829   const int *const * get_qmMMBondedIndx() { return qmMMBondedIndx; } ;
830   const char *const * get_qmDummyElement() { return qmDummyElement; } ;
831   
832   const int *const * get_qmMMChargeTarget() { return qmMMChargeTarget; } ;
833   const int * get_qmMMNumTargs() { return qmMMNumTargs; } ;
834   
835   int* get_qmLSSSize() { return qmLSSSize;} ;
836   int* get_qmLSSIdxs() { return qmLSSIdxs;} ;
837   Mass *get_qmLSSMass() { return qmLSSMass; } ;
838   int get_qmLSSFreq() { return qmLSSFreq; } ;
839   int get_qmLSSResSize() { return qmLSSResidueSize; } ;
840   int *get_qmLSSRefIDs() { return qmLSSRefIDs ; } ;
841   int *get_qmLSSRefSize() { return qmLSSRefSize ; } ;
842   Mass *get_qmLSSRefMass() {return qmLSSRefMass; } ;
843   std::map<int,int> &get_qmMMSolv() { return qmClassicSolv;} ;
844   
845   Bool get_qmReplaceAll() {return qmReplaceAll; } ;
846   
847   Bool get_noPC() { return qmNoPC; } ;
848   int get_qmMeNumBonds() { return qmMeNumBonds; } ;
849   int *get_qmMeMMindx() { return qmMeMMindx; } ;
850   Real *get_qmMeQMGrp() { return qmMeQMGrp; } ;
851   
852   int get_qmPCFreq() { return qmPCFreq; } ;
853   
854   int get_qmTotCustPCs() { return qmTotCustPCs; } ;
855   int *get_qmCustPCSizes()  { return qmCustPCSizes; } ;
856   int *get_qmCustomPCIdxs() { return qmCustomPCIdxs; } ;
857   
858   Bool get_qmcSMD() { return qmcSMD;} ;
859   int get_cSMDnumInst() { return cSMDnumInst;} ;
860   int const * get_cSMDindxLen() { return cSMDindxLen;} ;
861   int const * const * get_cSMDindex() {return cSMDindex;} ;
862   int const * const * get_cSMDpairs() {return cSMDpairs;} ;
863   const Real* get_cSMDKs() {return cSMDKs;} ;
864   const Real* get_cSMDVels() {return cSMDVels;} ;
865   const Real* get_cSMDcoffs() {return cSMDcoffs;} ;
866   
867   void prepare_qm(const char *pdbFileName, Parameters *params, ConfigList *cfgList) ;
868   void delete_qm_bonded() ;
869   // end QM
870   
871   
872   
873   Molecule(SimParameters *, Parameters *param);
874   Molecule(SimParameters *, Parameters *param, char *filename, ConfigList *cfgList=NULL);  
875   Molecule(SimParameters *simParams, Parameters *param, molfile_plugin_t *pIOHdl, void *pIOFileHdl, int natoms);
876   
877   Molecule(SimParameters *, Parameters *, Ambertoppar *);
878   void read_parm(Ambertoppar *);
879
880   Molecule(SimParameters *, Parameters *, const GromacsTopFile *);
881
882   ~Molecule();    //  Destructor
883
884   void send_Molecule(MOStream *);
885         //  send the molecular structure 
886         //  from the master to the clients
887
888   void receive_Molecule(MIStream *);
889         //  receive the molecular structure
890         //  from the master on a client
891   
892   void build_constraint_params(StringList *, StringList *, StringList *,
893              PDB *, char *);
894         //  Build the set of harmonic constraint 
895         // parameters
896
897 /* BEGIN gf */
898   void build_gridforce_params(StringList *, StringList *, StringList *, StringList *, PDB *, char *);
899         //  Build the set of gridForce-style force pars
900 /* END gf */
901
902   void build_movdrag_params(StringList *, StringList *, StringList *, 
903           PDB *, char *);
904         //  Build the set of moving drag pars
905
906   void build_rotdrag_params(StringList *, StringList *, StringList *,
907           StringList *, StringList *, StringList *,
908           PDB *, char *);
909         //  Build the set of rotating drag pars
910
911   void build_constorque_params(StringList *, StringList *, StringList *,
912              StringList *, StringList *, StringList *,
913              PDB *, char *);
914         //  Build the set of "constant" torque pars
915
916
917   void build_constant_forces(char *);
918         //  Build the set of constant forces
919
920   void build_langevin_params(BigReal coupling, BigReal drudeCoupling,
921       Bool doHydrogen);
922   void build_langevin_params(StringList *, StringList *, PDB *, char *);
923         //  Build the set of langevin dynamics parameters
924
925 #ifdef MEM_OPT_VERSION
926   void load_fixed_atoms(StringList *fixedFile);
927   void load_constrained_atoms(StringList *constrainedFile);
928 #endif
929
930   void build_fixed_atoms(StringList *, StringList *, PDB *, char *);
931         //  Determine which atoms are fixed (if any)
932
933   void build_stirred_atoms(StringList *, StringList *, PDB *, char *);
934         //  Determine which atoms are stirred (if any)
935
936   void build_extra_bonds(Parameters *parameters, StringList *file);
937
938 //fepb
939         void build_fep_flags(StringList *, StringList *, PDB *, char *, const char *);
940                                // selection of the mutant atoms
941         void delete_alch_bonded(void);
942 //fepe
943
944   /**
945    * Build the flags needed for solute scaling.
946    *
947    * A PDB file is read, indicating which atoms are to be scaled,
948    * and an array is maintained marking which are to be included.
949    * Each marked atom then has its corresponding van der Waals
950    * type number reassigned to enable extending the LJTable with
951    * scaled interaction values.
952    */
953   void build_ss_flags(
954       const StringList *ssfile,
955       /**< config "soluteScalingFile = my.pdb" for PDB filename */
956       const StringList *sscol,
957       /**< config "soluteScalingCol = O" for column of PDB ATOM records */
958       PDB *initial_pdb,         ///< the initial PDB file
959       const char *cwd           ///< current working directory
960       );
961
962   void build_exPressure_atoms(StringList *, StringList *, PDB *, char *);
963         //  Determine which atoms are excluded from
964                                 //  pressure (if any)
965
966   // Ported by JLai -- Original JE - Go -- Change the unsigned int to ints
967   void print_go_sigmas(); //  Print out Go sigma parameters
968   void build_go_sigmas(StringList *, char *);
969         //  Determine which atoms have Go forces applied
970         //  calculate sigmas from distances between Go atom pairs
971   void build_go_sigmas2(StringList *, char *);
972         //  Determine which atoms have Go forces applied
973         //  calculate sigmas from distances between Go atom pairs
974   void build_go_arrays(StringList *, char *);
975         //  Determine which atoms have Go forces applied
976   BigReal get_gro_force(BigReal, BigReal, BigReal, int, int) const;
977   BigReal get_gro_force2(BigReal, BigReal, BigReal, int, int, BigReal *, BigReal *) const;
978   BigReal get_go_force(BigReal, int, int, BigReal *, BigReal *) const;
979         //  Calculate the go force between a pair of atoms -- Modified to 
980         //  output Go energies
981   BigReal get_go_force_new(BigReal, int, int, BigReal *, BigReal *) const;
982         //  Calculate the go force between a pair of atoms
983   BigReal get_go_force2(BigReal, BigReal, BigReal, int, int, BigReal *, BigReal *) const;
984         //  Calculate the go force between a pair of atoms
985   Bool atoms_1to4(unsigned int, unsigned int);
986 // End of port -- JL  
987
988   void reloadCharges(float charge[], int n);
989
990         Bool is_lp(int);     // return true if atom is a lone pair
991         Bool is_drude(int);     // return true if atom is a Drude particle
992         Bool is_hydrogen(int);     // return true if atom is hydrogen
993         Bool is_oxygen(int);       // return true if atom is oxygen
994   Bool is_hydrogenGroupParent(int); // return true if atom is group parent
995   Bool is_water(int);        // return true if atom is part of water 
996   int  get_groupSize(int);     // return # atoms in (hydrogen) group
997         int get_mother_atom(int);  // return mother atom of a hydrogen
998
999   #ifdef MEM_OPT_VERSION
1000   //the way to get the cluster size if the atom ids of the cluster are
1001   //contiguous. The input parameter is the atom id that leads the cluster.
1002   int get_cluster_size_con(int aid) const { return clusterSigs[aid]; }  
1003   //the way to get the cluster size if the atoms ids of the cluster are
1004   //not contiguous. The input parameter is the cluster index.
1005   int get_cluster_size_uncon(int cIdx) const { return clusterSize[cIdx]; }
1006   int get_cluster_idx(int aid) const { return clusterSigs[aid]; }
1007   int get_num_clusters() const { return numClusters; }
1008   #else
1009   int get_cluster(int anum) const { return cluster[anum]; }
1010   int get_clusterSize(int anum) const { return clusterSize[anum]; }
1011   #endif
1012
1013 #ifndef MEM_OPT_VERSION
1014   const float *getOccupancyData() { return (const float *)occupancy; }
1015   void setOccupancyData(molfile_atom_t *atomarray);
1016   void freeOccupancyData() { delete [] occupancy; occupancy=NULL; }
1017
1018   const float *getBFactorData() { return (const float *)bfactor; }
1019   void setBFactorData(molfile_atom_t *atomarray);
1020   void freeBFactorData() { delete [] bfactor; bfactor=NULL; }
1021 #endif
1022
1023   //  Get the mass of an atom
1024   Real atommass(int anum) const
1025   {
1026     #ifdef MEM_OPT_VERSION
1027     return atomMassPool[eachAtomMass[anum]];
1028     #else
1029     return(atoms[anum].mass);
1030     #endif
1031   }
1032
1033   //  Get the charge of an atom
1034   Real atomcharge(int anum) const
1035   {
1036     #ifdef MEM_OPT_VERSION
1037     return atomChargePool[eachAtomCharge[anum]];
1038     #else
1039     return(atoms[anum].charge);
1040     #endif
1041   }
1042   
1043   //  Get the vdw type of an atom
1044   Index atomvdwtype(int anum) const
1045   {      
1046       return(atoms[anum].vdw_type);
1047   }
1048
1049   #ifndef MEM_OPT_VERSION
1050   //  Retrieve a bond structure
1051   Bond *get_bond(int bnum) const {return (&(bonds[bnum]));}
1052
1053   //  Retrieve an angle structure
1054   Angle *get_angle(int anum) const {return (&(angles[anum]));}
1055
1056   //  Retrieve an improper strutcure
1057   Improper *get_improper(int inum) const {return (&(impropers[inum]));}
1058
1059   //  Retrieve a dihedral structure
1060   Dihedral *get_dihedral(int dnum) const {return (&(dihedrals[dnum]));}
1061
1062   //  Retrieve a cross-term strutcure
1063   Crossterm *get_crossterm(int inum) const {return (&(crossterms[inum]));}
1064   #endif
1065
1066   // DRUDE: retrieve lphost structure
1067   Lphost *get_lphost(int atomid) const {
1068     // don't call unless simParams->drudeOn == TRUE
1069     // otherwise lphostIndexes array doesn't exist!
1070     int index = lphostIndexes[atomid];
1071     return (index != -1 ? &(lphosts[index]) : NULL);
1072   }
1073   // DRUDE
1074
1075   #ifndef MEM_OPT_VERSION
1076   Bond *getAllBonds() const {return bonds;}
1077   Angle *getAllAngles() const {return angles;}
1078   Improper *getAllImpropers() const {return impropers;}
1079   Dihedral *getAllDihedrals() const {return dihedrals;}
1080   Crossterm *getAllCrossterms() const {return crossterms;}
1081   #endif
1082
1083   // DRUDE: retrieve entire lphosts array
1084   Lphost *getAllLphosts() const { return lphosts; }
1085   // DRUDE
1086
1087   //  Retrieve a hydrogen bond donor structure
1088   Bond *get_donor(int dnum) const {return (&(donors[dnum]));}  
1089
1090   //  Retrieve a hydrogen bond acceptor structure
1091   Bond *get_acceptor(int dnum) const {return (&(acceptors[dnum]));} 
1092
1093   Bond *getAllDonors() const {return donors;}
1094   Bond *getAllAcceptors() const {return acceptors;}
1095
1096   //  Retrieve an exclusion structure
1097   #ifndef MEM_OPT_VERSION
1098   Exclusion *get_exclusion(int ex) const {return (&(exclusions[ex]));}
1099   #endif
1100
1101   //  Retrieve an atom type
1102   const char *get_atomtype(int anum) const
1103   {
1104     if (atomNames == NULL)
1105     {
1106       NAMD_die("Tried to find atom type on node other than node 0");
1107     }
1108
1109     #ifdef MEM_OPT_VERSION    
1110     return atomTypePool[atomNames[anum].atomtypeIdx];
1111     #else
1112     return(atomNames[anum].atomtype);
1113     #endif
1114   }
1115
1116   //  Lookup atom id from segment, residue, and name
1117   int get_atom_from_name(const char *segid, int resid, const char *aname) const;
1118
1119   //  Lookup number of atoms in residue from segment and residue
1120   int get_residue_size(const char *segid, int resid) const;
1121
1122   //  Lookup atom id from segment, residue, and index in residue
1123   int get_atom_from_index_in_residue(const char *segid, int resid, int index) const;
1124
1125   
1126   //  The following routines are used to get the list of bonds
1127   //  for a given atom.  This is used when creating the bond lists
1128   //  for the force objects  
1129
1130   #ifndef MEM_OPT_VERSION
1131   int32 *get_bonds_for_atom(int anum)
1132       { return bondsByAtom[anum]; } 
1133   int32 *get_angles_for_atom(int anum) 
1134       { return anglesByAtom[anum]; }
1135   int32 *get_dihedrals_for_atom(int anum) 
1136       { return dihedralsByAtom[anum]; }
1137   int32 *get_impropers_for_atom(int anum) 
1138       { return impropersByAtom[anum]; }  
1139   int32 *get_crossterms_for_atom(int anum) 
1140       { return crosstermsByAtom[anum]; }  
1141   int32 *get_exclusions_for_atom(int anum)
1142       { return exclusionsByAtom[anum]; }
1143   const int32 *get_full_exclusions_for_atom(int anum) const
1144       { return fullExclusionsByAtom[anum]; }
1145   const int32 *get_mod_exclusions_for_atom(int anum) const
1146       { return modExclusionsByAtom[anum]; }
1147   #endif
1148   
1149   //  Check for exclusions, either explicit or bonded.
1150         //  Returns 1 for full, 2 for 1-4 exclusions.
1151   #ifdef MEM_OPT_VERSION
1152   int checkExclByIdx(int idx1, int atom1, int atom2) const;
1153   const ExclusionCheck *get_excl_check_for_idx(int idx) const{      
1154       return &exclChkSigPool[idx];
1155   }
1156   #else
1157   int checkexcl(int atom1, int atom2) const;
1158
1159   const ExclusionCheck *get_excl_check_for_atom(int anum) const{      
1160       return &all_exclusions[anum];             
1161   }
1162   #endif
1163
1164 /* BEGIN gf */
1165   // Return true or false based on whether or not the atom
1166   // is subject to grid force
1167   Bool is_atom_gridforced(int atomnum, int gridnum) const
1168   {
1169       if (numGridforceGrids)
1170       {
1171           return(gridfrcIndexes[gridnum][atomnum] != -1);
1172       }
1173       else
1174       {
1175           return(FALSE);
1176       }
1177   }
1178 /* END gf */
1179
1180 #ifndef MEM_OPT_VERSION
1181   //  Return true or false based on whether the specified atom
1182   //  is constrained or not.
1183   Bool is_atom_constrained(int atomnum) const
1184   {
1185     if (numConstraints)
1186     {
1187       //  Check the index to see if it is constrained
1188       return(consIndexes[atomnum] != -1);
1189     }
1190     else
1191     {
1192       //  No constraints at all, so just return FALSE
1193       return(FALSE);
1194     }
1195   }
1196 #endif
1197
1198   //  Return true or false based on whether the specified atom
1199   //  is moving-dragged or not.
1200   Bool is_atom_movdragged(int atomnum) const
1201   {
1202     if (numMovDrag)
1203     {
1204       //  Check the index to see if it is constrained
1205       return(movDragIndexes[atomnum] != -1);
1206     }
1207     else
1208     {
1209       //  No constraints at all, so just return FALSE
1210       return(FALSE);
1211     }
1212   }
1213
1214   //  Return true or false based on whether the specified atom
1215   //  is rotating-dragged or not.
1216   Bool is_atom_rotdragged(int atomnum) const
1217   {
1218     if (numRotDrag)
1219     {
1220       //  Check the index to see if it is constrained
1221       return(rotDragIndexes[atomnum] != -1);
1222     }
1223     else
1224     {
1225       //  No constraints at all, so just return FALSE
1226       return(FALSE);
1227     }
1228   }
1229
1230   //  Return true or false based on whether the specified atom
1231   //  is "constant"-torqued or not.
1232   Bool is_atom_constorqued(int atomnum) const
1233   {
1234     if (numConsTorque)
1235     {
1236       //  Check the index to see if it is constrained
1237       return(consTorqueIndexes[atomnum] != -1);
1238     }
1239     else
1240     {
1241       //  No constraints at all, so just return FALSE
1242       return(FALSE);
1243     }
1244   }
1245
1246 #ifndef MEM_OPT_VERSION
1247   //  Get the harmonic constraints for a specific atom
1248   void get_cons_params(Real &k, Vector &refPos, int atomnum) const
1249   {
1250     k = consParams[consIndexes[atomnum]].k;
1251     refPos = consParams[consIndexes[atomnum]].refPos;
1252   }
1253 #endif
1254
1255 /* BEGIN gf */
1256   void get_gridfrc_params(Real &k, Charge &q, int atomnum, int gridnum) const
1257   {
1258       k = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].k;
1259       q = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].q;
1260   }
1261   
1262   GridforceGrid* get_gridfrc_grid(int gridnum) const
1263   {
1264       GridforceGrid *result = NULL;
1265       if (gridnum >= 0 && gridnum < numGridforceGrids) {
1266           result = gridfrcGrid[gridnum];
1267       }
1268       return result;
1269   }
1270   
1271   int set_gridfrc_grid(int gridnum, GridforceGrid *grid)
1272   {
1273       if (grid && gridnum >= 0 && gridnum < numGridforceGrids) {
1274           gridfrcGrid[gridnum] = grid;
1275           return 0;
1276       } else {
1277           return -1;
1278       }
1279   }
1280 /* END gf */
1281
1282   Real langevin_param(int atomnum) const
1283   {
1284     return(langevinParams ? langevinParams[atomnum] : 0.);
1285   }
1286
1287   //  Get the stirring constraints for a specific atom
1288   void get_stir_refPos(Vector &refPos, int atomnum) const
1289   {
1290     refPos = stirParams[stirIndexes[atomnum]].refPos;
1291   }
1292
1293
1294   void put_stir_startTheta(Real theta, int atomnum) const
1295   {
1296     stirParams[stirIndexes[atomnum]].startTheta = theta;
1297   }
1298
1299
1300   Real get_stir_startTheta(int atomnum) const
1301   {
1302     return stirParams[stirIndexes[atomnum]].startTheta;
1303   }
1304  
1305
1306   //  Get the moving drag factor for a specific atom
1307   void get_movdrag_params(Vector &v, int atomnum) const
1308   {
1309     v = movDragParams[movDragIndexes[atomnum]].v;
1310   }
1311
1312   //  Get the rotating drag pars for a specific atom
1313   void get_rotdrag_params(BigReal &v, Vector &a, Vector &p, 
1314         int atomnum) const
1315   {
1316     v = rotDragParams[rotDragIndexes[atomnum]].v;
1317     a = rotDragParams[rotDragIndexes[atomnum]].a;
1318     p = rotDragParams[rotDragIndexes[atomnum]].p;
1319   }
1320
1321   //  Get the "constant" torque pars for a specific atom
1322   void get_constorque_params(BigReal &v, Vector &a, Vector &p, 
1323         int atomnum) const
1324   {
1325     v = consTorqueParams[consTorqueIndexes[atomnum]].v;
1326     a = consTorqueParams[consTorqueIndexes[atomnum]].a;
1327     p = consTorqueParams[consTorqueIndexes[atomnum]].p;
1328   }
1329
1330 //fepb
1331   unsigned char get_fep_type(int anum) const
1332   {
1333     return(fepAtomFlags[anum]);
1334   }
1335
1336 //soluteScaling
1337         unsigned char get_ss_type(int anum) const
1338         {
1339                 return(ssAtomFlags[anum]);
1340         }
1341 //soluteScaling
1342
1343   /* BKR - Get the FEP type (i.e. 0, 1, or 2) of a bonded _interaction_ based 
1344      on the atom indices of the atoms involved (internally converted to FEP 
1345      types) and the order of the bonded interaction (i.e. 2, 3, or 4). This 
1346      automatically accounts for whether or not purely alchemical interactions
1347      are being scaled (e.g. bonds between two atoms of fep_type 1).
1348
1349      The logic here is admittedly a bit opaque. When the fep_type's are back
1350      mapped to -1,0,1, we can use the sum of all of the types to determine the 
1351      type of interaction for all bonded term types:
1352
1353      0  - only non-alchemical atoms involved
1354      >0 - _atleast_ one appearing atom
1355      <0 - _atleast_ one disappearing atom
1356
1357      If the magnitude of the sum is equal to the interaction order (i.e. 2, 3, 
1358      or 4), then it is a _purely_ alchemical interaction and it may be 
1359      desireable to retain it for sampling purposes (note that this adds a 
1360      constant to the free energy that will almost always differ at the 
1361      endpoints).  In order to avoid unexpected instability these interactions 
1362      are retained by default, but can be scaled along with everything else by 
1363      setting "alchBondDecouple on".
1364
1365      NB: All pure alchemical interactions beyond order 2 are ALWAYS discarded
1366      by alchemify. This is good, because higher order interactions would break
1367      the above logic. For example, if we had an angle term between atoms of 
1368      types (1,1,-1) the sum would be 1, but this term should receive no scaling
1369      because it involves groups -1 and 1 but not 0.
1370   */
1371   int get_fep_bonded_type(const int *atomID, unsigned int order) const
1372   {
1373     int typeSum = 0;
1374     for ( int i=0; i < order; ++i ) {
1375       typeSum += (fepAtomFlags[atomID[i]] == 2 ? -1 : fepAtomFlags[atomID[i]]);
1376     }
1377     // Increase the cutoff if scaling purely alchemical bonds. 
1378     // This really only applies when order = 2.
1379     if ( simParams->alchBondDecouple ) order++;
1380
1381     if ( typeSum == 0 ) return 0; // Most interactions get caught here.
1382     else if ( 0 < typeSum && typeSum < order ) return 1;
1383     else if ( 0 > typeSum && typeSum > -order ) return 2;
1384
1385     if ( simParams->alchBondDecouple ) {
1386       // Alchemify should always keep this from bombing, but just in case...
1387       NAMD_die("Unexpected alchemical bonded interaction!");
1388     }
1389   }
1390 //fepe
1391
1392 #ifndef MEM_OPT_VERSION
1393   Bool is_atom_fixed(int atomnum) const
1394   {
1395     return (numFixedAtoms && fixedAtomFlags[atomnum]);
1396   }
1397 #else
1398   //Since binary search is more expensive than direct array access,
1399   //and this function is usually called for consecutive atoms in this context,
1400   //the *listIdx returns the index to the range of atoms [aid1, aid2]
1401   //that are fixed. If the atom aid is fixed, then aid1=<aid<=aid2;
1402   //If the atom aid is not fixed, then aid1 indicates the smallest fixed atom
1403   //id that is larger than aid; so the listIdx could be equal the size of
1404   //fixedAtomsSet. --Chao Mei
1405   Bool is_atom_in_set(AtomSetList *localAtomsSet, int aid, int *listIdx) const;
1406   inline Bool is_atom_fixed(int aid, int *listIdx=NULL) const {
1407     return is_atom_in_set(fixedAtomsSet,aid,listIdx);
1408   }
1409   inline Bool is_atom_constrained(int aid, int *listIdx=NULL) const {
1410     return is_atom_in_set(constrainedAtomsSet,aid,listIdx);
1411   }
1412 #endif
1413         
1414   Bool is_atom_stirred(int atomnum) const
1415   {
1416     if (numStirredAtoms)
1417     {
1418       //  Check the index to see if it is constrained
1419       return(stirIndexes[atomnum] != -1);
1420     }
1421     else
1422     {
1423       //  No constraints at all, so just return FALSE
1424       return(FALSE);
1425     }
1426   }
1427   
1428
1429   Bool is_group_fixed(int atomnum) const
1430   {
1431     return (numFixedAtoms && (fixedAtomFlags[atomnum] == -1));
1432   }
1433   Bool is_atom_exPressure(int atomnum) const
1434   {
1435     return (numExPressureAtoms && exPressureAtomFlags[atomnum]);
1436   }
1437   // 0 if not rigid or length to parent, for parent refers to H-H length
1438   // < 0 implies MOLLY but not SHAKE, > 1 implies both if MOLLY is on
1439   Real rigid_bond_length(int atomnum) const
1440   {
1441     return(rigidBondLengths[atomnum]);
1442   }
1443
1444   void print_atoms(Parameters *); 
1445         //  Print out list of atoms
1446   void print_bonds(Parameters *); 
1447         //  Print out list of bonds
1448   void print_exclusions();//  Print out list of exclusions
1449
1450 public:  
1451   int isOccupancyValid, isBFactorValid;
1452
1453 #ifdef MEM_OPT_VERSION
1454   //read the per-atom file for the memory optimized version where the file 
1455   //name already exists the range of atoms to be read are 
1456   //[fromAtomID, toAtomID].
1457   void read_binary_atom_info(int fromAtomID, int toAtomID, InputAtomList &inAtoms);
1458
1459   int64 getNumCalcExclusions(){return numCalcExclusions;}
1460   void setNumCalcExclusions(int64 x){numCalcExclusions= x;}
1461
1462   Index getEachAtomMass(int i){return eachAtomMass[i];}
1463   Index getEachAtomCharge(int i){return eachAtomCharge[i];}
1464
1465   ExclSigID getAtomExclSigId(int aid) const {
1466       return eachAtomExclSig[aid];
1467   }
1468
1469   Real *getAtomMassPool(){return atomMassPool;}
1470   Real *getAtomChargePool(){return atomChargePool;}
1471   AtomCstInfo *getAtoms(){return atoms;}  
1472
1473   int atomSigPoolSize;
1474   AtomSignature *atomSigPool;
1475
1476   /* All the following are temporary variables for reading the compressed psf file */
1477   //declarations for atoms' constant information  
1478   int segNamePoolSize; //Its value is usually less than 5
1479   char **segNamePool; //This seems not to be important, but it only occupied very little space.
1480
1481   int resNamePoolSize;
1482   char **resNamePool;
1483
1484   int atomNamePoolSize;
1485   char **atomNamePool;
1486
1487   int atomTypePoolSize;
1488   char **atomTypePool;
1489
1490   int chargePoolSize;
1491   Real *atomChargePool;
1492
1493   int massPoolSize;
1494   Real *atomMassPool;
1495
1496   AtomSigID getAtomSigId(int aid) {
1497       return eachAtomSig[aid]; 
1498   }
1499
1500   //Indicates the size of both exclSigPool and exclChkSigPool
1501   int exclSigPoolSize;
1502   //this will be deleted after build_lists_by_atom
1503   ExclusionSignature *exclSigPool;
1504   //This is the final data structure we want to store  
1505   ExclusionCheck *exclChkSigPool;
1506
1507   void addNewExclSigPool(const std::vector<ExclusionSignature>&);  
1508
1509   void delEachAtomSigs(){    
1510       //for NAMD-smp version, only one Molecule object is held
1511       //on each node, therefore, only one deletion operation should
1512       //be taken on a node, otherwise, there possibly would be some
1513       //wierd memory problems. The same reason applies to other deletion
1514       //operations inside the Molecule object.   
1515       if(CmiMyRank()) return;
1516
1517       delete [] eachAtomSig;
1518       delete [] eachAtomExclSig;
1519       eachAtomSig = NULL;
1520       eachAtomExclSig = NULL;
1521   }
1522
1523   void delChargeSpace(){
1524       if(CmiMyRank()) return;
1525
1526       delete [] atomChargePool;
1527       delete [] eachAtomCharge;
1528       atomChargePool = NULL;
1529       eachAtomCharge = NULL;
1530   }
1531   
1532   void delMassSpace(){
1533       if(CmiMyRank()) return;
1534
1535       delete [] atomMassPool;
1536       delete [] eachAtomMass;
1537       atomMassPool = NULL;
1538       eachAtomMass = NULL;
1539   }
1540   
1541   void delClusterSigs() {
1542       if(CmiMyRank()) return;      
1543
1544       delete [] clusterSigs;
1545       clusterSigs = NULL;
1546   }
1547
1548   void delAtomNames(){
1549       if(CmiMyRank()) return;
1550       delete [] atomNamePool;
1551       delete [] atomNames;
1552       atomNamePool = NULL;
1553       atomNames = NULL;
1554   }
1555
1556   void delFixedAtoms(){
1557       if(CmiMyRank()) return;
1558       delete fixedAtomsSet;
1559       fixedAtomsSet = NULL;
1560   }
1561
1562 private:
1563   Index insert_new_mass(Real newMass);
1564
1565 #endif
1566
1567 // Go stuff
1568 public:
1569
1570 GoValue go_array[MAX_GO_CHAINS*MAX_GO_CHAINS];    //  Array of Go params -- JLai
1571 int go_indices[MAX_GO_CHAINS+1];        //  Indices from chainIDS to go array -- JLai
1572 int NumGoChains;                        //  Number of Go chain types -- JLai
1573
1574 // Declares and initializes Go variables
1575 void goInit();
1576
1577 // Builds explicit Gromacs pairs
1578 void build_gro_pair();
1579
1580 // Builds the initial Go parameters 
1581 void build_go_params(StringList *);
1582
1583 //  Read Go parameter file
1584 void read_go_file(char *);
1585
1586 //  Get Go cutoff for a given chain type pair
1587 Real get_go_cutoff(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].cutoff; };
1588
1589 //  Get Go epsilonRep for a given chain type pair
1590 Real get_go_epsilonRep(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilonRep; };
1591
1592 //  Get Go sigmaRep for a given chain type pair
1593 Real get_go_sigmaRep(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].sigmaRep; };
1594
1595 //  Get Go epsilon for a given chain type pair
1596 Real get_go_epsilon(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilon; };
1597
1598 //  Get Go exp_a for a given chain type pair
1599 int get_go_exp_a(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_a; };
1600
1601 //  Get Go exp_b for a given chain type pair
1602 int get_go_exp_b(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_b; };
1603
1604 //  Get Go exp_rep for a given chain type pair
1605 int get_go_exp_rep(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_rep; };
1606
1607 //  Whether residue IDs with this difference are restricted
1608 Bool go_restricted(int, int, int);
1609
1610 // Prints Go Params
1611 void print_go_params();
1612
1613 void initialize();
1614
1615 void send_GoMolecule(MOStream *);
1616 //  send the molecular structure 
1617 //  from the master to the clients
1618
1619 void receive_GoMolecule(MIStream *);
1620 //  receive the molecular structure
1621 //  from the master on a client
1622 };
1623
1624 #endif
1625