Molecule.h

Go to the documentation of this file.
00001 
00007 /*
00008    This class is used to store all of the structural       
00009    information for a simulation.  It reads in this information 
00010    from a .psf file, cross checks and obtains some information
00011    from the Parameters object that is passed in, and then     
00012    stores all this information for later use. 
00013 */
00014 
00015 
00016 #ifndef MOLECULE_H
00017 
00018 #define MOLECULE_H
00019 
00020 #include "parm.h"
00021 
00022 #include "common.h"
00023 #include "NamdTypes.h"
00024 #include "structures.h"
00025 #include "ConfigList.h"
00026 #include "Vector.h"
00027 #include "UniqueSet.h"
00028 #include "Hydrogen.h"
00029 #include "GromacsTopFile.h"
00030 /* BEGIN gf */
00031 #include "GridForceGrid.h"
00032 #include "Tensor.h"
00033 /* END gf */
00034 
00035 // Go -- JLai
00036 #define MAX_GO_CHAINS          10
00037 #define MAX_RESTRICTIONS       10
00038 // End of Go
00039 
00040 #include "molfile_plugin.h"
00041 
00042 #include <vector>
00043 
00044 class SimParameters;
00045 class Parameters;
00046 class ConfigList;
00047 class PDB;
00048 class MIStream;
00049 class MOStream;
00050 
00051 class ExclElem;
00052 class BondElem;
00053 class AngleElem;
00054 class DihedralElem;
00055 class ImproperElem;
00056 class TholeElem;  // Drude model
00057 class AnisoElem;  // Drude model
00058 class CrosstermElem;
00059 // JLai
00060 class GromacsPairElem;
00061 // End of JLai
00062 class ResidueLookupElem;
00063 
00064 struct OutputAtomRecord;
00065 
00066 template<class Type> class ObjectArena;
00067 
00068 class ExclusionCheck {
00069 public:
00070   int32 min,max;
00071   char *flags;
00072 
00073   ExclusionCheck(){
00074       min=0;
00075       max=-1;
00076       flags = NULL;
00077   }
00078 private:
00079   ExclusionCheck(const ExclusionCheck& chk){
00080   }
00081   ExclusionCheck &operator=(const ExclusionCheck& chk){
00082     return *this;
00083   }
00084 };
00085 #define EXCHCK_FULL 1
00086 #define EXCHCK_MOD 2
00087 
00088 class ResidueLookupElem
00089 {
00090 public:
00091   char mySegid[11];
00092   ResidueLookupElem *next;      // stored as a linked list
00093   int firstResid;               // valid resid is >= firstResid
00094   int lastResid;                // valid resid is <= lastResid
00095   ResizeArray<int> atomIndex;   // 0-based index for first atom in residue
00096 
00097   ResidueLookupElem(void) { next = 0; firstResid = -1; lastResid = -1; }
00098   ~ResidueLookupElem(void) { delete next; }
00099   int lookup(const char *segid, int resid, int *begin, int *end) const;
00100   ResidueLookupElem* append(const char *segid, int resid, int aid);
00101 };
00102 
00103 // Ported by JLai -- JE
00104 typedef struct go_val
00105 {
00106   Real epsilon;      // Epsilon
00107   int exp_a;         // First exponent for attractive L-J term
00108   int exp_b;         // Second exponent for attractive L-J term
00109   int exp_rep;       // Exponent for repulsive L-J term
00110   Real sigmaRep;     // Sigma for repulsive term
00111   Real epsilonRep;   // Epsilon for replusive term
00112   Real cutoff;       // Cutoff distance for Go calculation
00113   int  restrictions[MAX_RESTRICTIONS];  //  List of residue ID differences to be excluded from Go calculation
00114 } GoValue;
00115 
00116 typedef struct go_pair
00117 {
00118   int goIndxA;       // indexA
00119   int goIndxB;       // indexB
00120   double A;          // double A for the LJ pair
00121   double B;          // double B for the LJ pair
00122 } GoPair;
00123 // End of port - JL
00124 
00125 #define QMLSSMODEDIST 1
00126 #define QMLSSMODECOM 2
00127 
00128 #define QMFormatORCA 1
00129 #define QMFormatMOPAC 2
00130 #define QMFormatUSR 3
00131 
00132 #define QMCHRGNONE 0
00133 #define QMCHRGMULLIKEN 1
00134 #define QMCHRGCHELPG 2
00135 
00136 #define QMSCHEMECS 1
00137 #define QMSCHEMERCD 2
00138 #define QMSCHEMEZ1 3
00139 #define QMSCHEMEZ2 4
00140 #define QMSCHEMEZ3 5
00141 
00142 //only used for compressing the molecule information
00143 typedef struct seg_resid
00144 {
00145     char segname[11];
00146     int resid;
00147 }AtomSegResInfo;
00148 
00149 // List maintaining the global atom indicies sorted by helix groups.
00150 class Molecule
00151 {
00152  private:
00153 typedef struct constraint_params
00154 {
00155    Real k;    //  Force constant
00156    Vector refPos; //  Reference position for restraint
00157 } ConstraintParams;
00158 
00159 
00160 
00161 /* BEGIN gf */
00162 typedef struct gridfrc_params
00163 {
00164     Real k;     // force multiplier
00165     Charge q;   // charge
00166 } GridforceParams;
00167 /* END gf */
00168 
00169 
00170 typedef struct stir_params
00171 {
00172   Real startTheta;
00173   Vector refPos;  //  Reference position for stirring
00174 } StirParams;
00175 
00176 typedef struct movdrag_params
00177 {
00178   Vector v;            //  Linear velocity (A/step)
00179 } MovDragParams;
00180 
00181 
00182 typedef struct rotdrag_params
00183 {
00184    Real v;              //  Angular velocity (deg/step)
00185    Vector a;            //  Rotation axis
00186    Vector p;            //  Rotation pivot point
00187 } RotDragParams;
00188 
00189 typedef struct constorque_params
00190 {
00191    Real v;              //  "Torque" value (Kcal/(mol*A^2))
00192    Vector a;            //  Rotation axis
00193    Vector p;            //  Rotation pivot point
00194 } ConsTorqueParams;
00195 
00196 #ifdef MEM_OPT_VERSION
00197 //indicate a range of atoms from aid1 to aid2
00198 struct AtomSet{
00199     AtomID aid1;
00200     AtomID aid2;
00201 
00202     int operator < (const AtomSet &a) const{
00203                 return aid1 < a.aid1;
00204         }
00205 };
00206 typedef ResizeArray<AtomSet> AtomSetList;
00207 
00208   void load_atom_set(StringList *setfile, const char *setname,
00209         int *numAtomsInSet, Molecule::AtomSetList **atomsSet) const;
00210 
00211 #endif
00212 
00213 friend class ExclElem;
00214 friend class BondElem;
00215 friend class AngleElem;
00216 friend class DihedralElem;
00217 friend class ImproperElem;
00218 friend class TholeElem;  // Drude model
00219 friend class AnisoElem;  // Drude model
00220 friend class CrosstermElem;
00221 // JLai
00222 friend class GromacsPairElem;
00223 // End of JLai
00224 friend class WorkDistrib;
00225 
00226 private:
00227 
00228 #ifndef MEM_OPT_VERSION
00229         Atom *atoms;    //  Array of atom structures
00230         ObjectArena<char> *nameArena;
00231         AtomNameInfo *atomNames;//  Array of atom name info.  Only maintained on node 0 for VMD interface
00232         //replaced by atom signatures
00233         Bond *bonds;    //  Array of bond structures
00234         Angle *angles;    //  Array of angle structures
00235         Dihedral *dihedrals;  //  Array of dihedral structures
00236         Improper *impropers;  //  Array of improper structures                          
00237         Crossterm *crossterms;  //  Array of cross-term structures
00238         GromacsPair *gromacsPair; //  Array of gromacs-pair structures
00239 
00240         //These will be replaced by exclusion signatures
00241         Exclusion *exclusions;  //  Array of exclusion structures
00242         UniqueSet<Exclusion> exclusionSet;  //  Used for building
00243 
00244         int32 *cluster;   //  first atom of connected cluster
00245 
00246         ObjectArena<int32> *tmpArena;
00247         int32 **bondsWithAtom;  //  List of bonds involving each atom
00248         ObjectArena<int32> *arena;
00249 
00250         //function is replaced by atom signatures
00251         int32 **bondsByAtom;  //  List of bonds owned by each atom
00252         int32 **anglesByAtom;     //  List of angles owned by each atom
00253         int32 **dihedralsByAtom;  //  List of dihedrals owned by each atom
00254         int32 **impropersByAtom;  //  List of impropers owned by each atom
00255         int32 **crosstermsByAtom;  //  List of crossterms owned by each atom
00256 
00257         int32 **exclusionsByAtom; //  List of exclusions owned by each atom
00258         int32 **fullExclusionsByAtom; //  List of atoms excluded for each atom
00259         int32 **modExclusionsByAtom; //  List of atoms modified for each atom
00260 // JLai
00261         int32 **gromacsPairByAtom; // gromacsPair exclusion list which by definition should not have any exclusions (still not sure if it should be empty or zeroed out)
00262 // End of JLai
00263 
00264         ObjectArena<char> *exclArena;
00265         ExclusionCheck *all_exclusions;
00266 
00267         // DRUDE
00268         int32 **tholesByAtom;  // list of Thole correction terms owned by each atom
00269         int32 **anisosByAtom;  // list of anisotropic terms owned by each atom
00270         // DRUDE
00271 
00272 #else
00273         //Indexing to constant pools to save space
00274         AtomCstInfo *atoms;
00275         Index *eachAtomMass; //after initialization, this could be freed (possibly)
00276         Index *eachAtomCharge; //after initialization, this could be freed (possibly)
00277         AtomNameIdx *atomNames;
00278         ObjectArena<char> *nameArena; //the space for names to be allocated  
00279 
00280         //A generally true assumption: all atoms are arranged in the order of clusters. In other words,
00281         //all atoms for a cluster must appear before/after any atoms in other clusters
00282         //The first atom in the cluster (which has the lowest atom id) stores the cluster size
00283         //other atoms in the cluster stores -1
00284         int32 *clusterSigs;
00285         int32 numClusters;
00286         
00287         AtomSigID *eachAtomSig;
00288         ExclSigID *eachAtomExclSig;
00289 
00290     AtomSetList *fixedAtomsSet;    
00291     AtomSetList *constrainedAtomsSet;    
00292 #endif
00293  
00294   ResidueLookupElem *resLookup; // find residues by name
00295 
00296   AtomSegResInfo *atomSegResids; //only used for generating compressed molecule info
00297 
00298   Bond *donors;         //  Array of hydrogen bond donor structures
00299   Bond *acceptors;  //  Array of hydrogen bond acceptor
00300 
00301   // DRUDE
00302   DrudeConst *drudeConsts;  // supplement Atom data (length of Atom array)
00303   Thole *tholes;            // Thole (screened Coulomb) interactions
00304   Aniso *anisos;            // anisotropic terms
00305   Lphost *lphosts;          // lone pair hosts
00306   int32 *lphostIndexes;     // index for each LP into lphosts array
00307   // DRUDE
00308 
00309   int32 *consIndexes; //  Constraint indexes for each atom
00310   ConstraintParams *consParams;
00311 
00312 /* BEGIN gf */
00313   int32 **gridfrcIndexes;
00314   GridforceParams **gridfrcParams;
00315   GridforceGrid **gridfrcGrid;
00316 /* END gf */
00317 
00318         //  Parameters for each atom constrained
00319   int32 *stirIndexes; //Stirring indexes for each atoms
00320   StirParams *stirParams;
00321                           // Paramters for each atoms stirred
00322   int32 *movDragIndexes;  //  Moving drag indexes for each atom
00323   MovDragParams *movDragParams;
00324                                 //  Parameters for each atom moving-dragged
00325   int32 *rotDragIndexes;  //  Rotating drag indexes for each atom
00326   RotDragParams *rotDragParams;
00327                                 //  Parameters for each atom rotation-dragged
00328 
00329   Real *langevinParams;   //  b values for langevin dynamics
00330   int32 *fixedAtomFlags;  //  1 for fixed, -1 for fixed group, else 0
00331   int32 *exPressureAtomFlags; // 1 for excluded, -1 for excluded group.
00332 
00333   //In the memory optimized version: it will be NULL if the general
00334   //true assumption mentioned above holds. If not, its size is numClusters.
00335   //In the ordinary version: its size is numAtoms, and indicates the size
00336   //of connected cluster or 0.
00337   int32 *clusterSize;                            
00338 
00339         Real *rigidBondLengths;  //  if H, length to parent or 0. or
00340         //  if not H, length between children or 0.
00341 //fepb
00342         unsigned char *fepAtomFlags; 
00343 //fepe
00344 
00345   //occupancy and bfactor data from plugin-based IO implementation of loading structures
00346   float *occupancy;
00347   float *bfactor;
00348 
00349 
00350   // added during the trasition from 1x to 2
00351   SimParameters *simParams;
00352   Parameters *params;
00353 
00354 private:
00355         void initialize(SimParameters *, Parameters *param);
00356         // Sets most data to zero
00357 
00358   //LCPO
00359   int *lcpoParamType;
00360 
00361 #ifndef MEM_OPT_VERSION
00362         void read_psf_file(char *, Parameters *);
00363         //  Read in a .psf file given
00364         //  the filename and the parameter
00365         //  object to use          
00366   
00367   void read_atoms(FILE *, Parameters *);
00368         //  Read in atom info from .psf
00369   void read_bonds(FILE *, Parameters *);
00370         //  Read in bond info from .psf
00371   void read_angles(FILE *, Parameters *);
00372         //  Read in angle info from .psf
00373   void read_dihedrals(FILE *, Parameters *);
00374         //  Read in dihedral info from .psf
00375   void read_impropers(FILE *, Parameters *);
00376         //  Read in improper info from .psf
00377   void read_crossterms(FILE *, Parameters *);
00378         //  Read in cross-term info from .psf
00379   void read_donors(FILE *);
00380         //  Read in hydrogen bond donors from .psf
00381   void read_acceptors(FILE *);
00382         //  Read in hydrogen bond acceptors from .psf
00383   void read_exclusions(FILE *);
00384         //  Read in exclusion info from .psf
00385   // JLai August 16th, 2012 Modifications
00386   void read_exclusions(int*, int*, int);
00387   /* Read in exclusion array and sort entries */
00388   static bool goPairCompare (GoPair, GoPair);
00389   // End of JLai August 16th, 2012 Modifications
00390 
00391 
00392   // DRUDE: PSF reading
00393   void read_lphosts(FILE *);
00394         //  Read in lone pair hosts from Drude PSF
00395   void read_anisos(FILE *);
00396         //  Read in anisotropic terms from Drude PSF
00397   // DRUDE
00398 
00399   //LCPO
00400   //input type is Charmm/Amber/other
00401   //0 - Charmm/Xplor
00402   //1 - Amber TODO
00403   //2 - Plugin TODO
00404   //3 - Gromacs TODO
00405   void assignLCPOTypes(int inputType);
00406 
00407   //pluginIO-based loading atoms' structure
00408   void plgLoadAtomBasics(molfile_atom_t *atomarray);
00409   void plgLoadBonds(int *from, int *to); //atom index is 1-based in the parameters
00410   void plgLoadAngles(int *plgAngles);
00411   void plgLoadDihedrals(int *plgDihedrals);
00412   void plgLoadImpropers(int *plgImpropers);
00413   void plgLoadCrossterms(int *plgCterms);
00414 
00415         //  List of all exclusions, including
00416         //  explicit exclusions and those calculated
00417         //  from the bonded structure based on the
00418         //  exclusion policy.  Also includes 1-4's.
00419   void build_lists_by_atom();
00420         //  Build the list of structures by atom
00421 
00422   void build12excl(void);
00423   void build13excl(void);
00424   void build14excl(int);
00425 
00426   // DRUDE: extend exclusions for Drude and LP
00427   void build_inherited_excl(int);
00428   // DRUDE
00429   void stripFepExcl(void);
00430 
00431   void build_exclusions();
00432   // analyze the atoms, and determine which are oxygen, hb donors, etc.
00433   // this is called after a molecule is sent our (or received in)
00434   void build_atom_status(void);
00435 
00436 #else
00437   //the method to load the signatures of atoms etc. (i.e. reading the file in 
00438   //text fomrat of the compressed psf file)
00439   void read_mol_signatures(char *fname, Parameters *params, ConfigList *cfgList=0);     
00440   void load_one_inputatom(int aid, OutputAtomRecord *one, InputAtom *fAtom);
00441   void build_excl_check_signatures();  
00442 #endif
00443 
00444         void read_parm(const GromacsTopFile *);  
00445           
00446 public:
00447   // DRUDE
00448   int is_drude_psf;      // flag for reading Drude PSF
00449   int is_lonepairs_psf;  // flag for reading lone pairs from PSF
00450   // DRUDE
00451 
00452   // data for TIP4P
00453   Real r_om;
00454   Real r_ohc;
00455 
00456   // data for tail corrections
00457   BigReal tail_corr_ener;
00458   BigReal tail_corr_virial;
00459   BigReal tail_corr_dUdl_1;
00460   BigReal tail_corr_virial_1;
00461   void compute_LJcorrection();
00462   BigReal getEnergyTailCorr(const BigReal);
00463   BigReal getVirialTailCorr(const BigReal);
00464 
00465   int const * getLcpoParamType() {
00466     return lcpoParamType;
00467   }
00468 
00469   BigReal GetAtomAlpha(int i) const {
00470     return drudeConsts[i].alpha;
00471   }
00472 
00473 #ifdef MEM_OPT_VERSION
00474   AtomCstInfo *getAtoms() const { return atoms; }
00475   AtomNameIdx *getAtomNames() const { return atomNames; }
00476 #else
00477   Atom *getAtoms () const { return atoms; }
00478   AtomNameInfo *getAtomNames() const { return atomNames; }
00479 #endif
00480 
00481   AtomSegResInfo *getAtomSegResInfo() const { return atomSegResids; }
00482   
00483   // return number of fixed atoms, guarded by SimParameters
00484   int num_fixed_atoms() const {
00485     // local variables prefixed by s_
00486     int s_NumFixedAtoms = (simParams->fixedAtomsOn ? numFixedAtoms : 0);
00487     return s_NumFixedAtoms;  // value is "turned on" SimParameters
00488   }
00489 
00490   int num_fixed_groups() const {
00491     // local variables prefixed by s_
00492     int s_NumFixedAtoms = num_fixed_atoms();
00493     int s_NumFixedGroups = (s_NumFixedAtoms ? numFixedGroups : 0);
00494     return s_NumFixedGroups;  // value is "turned on" SimParameters
00495   }
00496 
00497   int64_t num_group_deg_freedom() const {
00498     // local variables prefixed by s_
00499     int64_t s_NumGroupDegFreedom = 3 * (int64_t) numHydrogenGroups;
00500     int64_t s_NumFixedAtoms = num_fixed_atoms();
00501     int s_NumFixedGroups = num_fixed_groups();
00502     if (s_NumFixedGroups) s_NumGroupDegFreedom -= 3 * s_NumFixedGroups;
00503     if ( ! (s_NumFixedAtoms || numConstraints
00504           || simParams->comMove || simParams->langevinOn) ) {
00505       s_NumGroupDegFreedom -= 3;
00506     }
00507     return s_NumGroupDegFreedom;
00508   }
00509 
00510   int64_t num_deg_freedom(int isInitialReport = 0) const {
00511     // local variables prefixed by s_
00512     int64_t s_NumDegFreedom = 3 * (int64_t) numAtoms;
00513     int64_t s_NumFixedAtoms = num_fixed_atoms();
00514     if (s_NumFixedAtoms) s_NumDegFreedom -= 3 * s_NumFixedAtoms;
00515     if (numLonepairs) s_NumDegFreedom -= 3 * numLonepairs;
00516     if ( ! (s_NumFixedAtoms || numConstraints
00517           || simParams->comMove || simParams->langevinOn) ) {
00518       s_NumDegFreedom -= 3;
00519     }
00520     if ( ! isInitialReport && simParams->pairInteractionOn) {
00521       //
00522       // DJH: a kludge?  We want to initially report system degrees of freedom
00523       //
00524       // this doesn't attempt to deal with fixed atoms or constraints
00525       s_NumDegFreedom = 3 * numFepInitial;
00526     }
00527     int s_NumFixedRigidBonds = 
00528       (simParams->fixedAtomsOn ? numFixedRigidBonds : 0);
00529     if (simParams->watmodel == WAT_TIP4) {
00530       // numLonepairs is subtracted here because all lonepairs have a rigid bond
00531       // to oxygen, but all of the LP degrees of freedom are dealt with above
00532       s_NumDegFreedom -= (numRigidBonds - s_NumFixedRigidBonds - numLonepairs);
00533     }
00534     else {
00535       // Non-polarized systems don't have LPs.
00536       // For Drude model, bonds that attach LPs are not counted as rigid;
00537       // LPs have already been subtracted from degrees of freedom.
00538       s_NumDegFreedom -= (numRigidBonds - s_NumFixedRigidBonds);
00539     }
00540     return s_NumDegFreedom;
00541   }
00542 
00543   int numAtoms;   //  Number of atoms                   
00544 
00545   int numRealBonds;   //  Number of bonds for exclusion determination
00546   int numBonds;   //  Number of bonds calculated, including extras
00547   int numAngles;    //  Number of angles
00548   int numDihedrals; //  Number of dihedrals
00549   int suspiciousAlchBonds;    //  angles dropped due to crossing FEP partitions
00550   int alchDroppedAngles;    //  angles dropped due to crossing FEP partitions
00551   int alchDroppedDihedrals; //  dihedrals dropped due to crossing FEP partitions
00552   int alchDroppedImpropers; //  impropers dropped due to crossing FEP partitions
00553   int numImpropers; //  Number of impropers
00554   int numCrossterms; //  Number of cross-terms
00555   int numDonors;          //  Number of hydrogen bond donors
00556   int numAcceptors; //  Number of hydrogen bond acceptors
00557   int numExclusions;  //  Number of exclusions
00558   int numTotalExclusions; //  Real Total Number of Exclusions // hack
00559 
00560   // DRUDE
00561   int numLonepairs; // Number of lone pairs
00562   int numDrudeAtoms;  // Number of Drude particles
00563   int numTholes;  // Number of Thole terms
00564   int numAnisos;  // Number of anisotropic terms
00565   int numLphosts;  // Number of lone pair hosts
00566   // DRUDE
00567   
00568   int numConstraints; //  Number of atoms constrained
00569 /* BEGIN gf */
00570   int numGridforceGrids;//  Number of gridforce grids
00571   int *numGridforces;   //  Number of atoms in gridforce file (array, one per grid)
00572 /* END gf */
00573   int numMovDrag;         //  Number of atoms moving-dragged
00574   int numRotDrag;         //  Number of atoms rotating-dragged
00575   int numConsTorque;  //  Number of atoms "constant"-torqued
00576   int numFixedAtoms;  //  Number of fixed atoms
00577   int numStirredAtoms;  //  Number of stirred atoms
00578   int numExPressureAtoms; //  Number of atoms excluded from pressure
00579   int numHydrogenGroups;  //  Number of hydrogen groups
00580   int maxHydrogenGroupSize;  //  Max atoms per hydrogen group
00581   int numMigrationGroups;  //  Number of migration groups
00582   int maxMigrationGroupSize;  //  Max atoms per migration group
00583   int numFixedGroups; //  Number of totally fixed hydrogen groups
00584   int numRigidBonds;  //  Number of rigid bonds
00585   int numFixedRigidBonds; //  Number of rigid bonds between fixed atoms
00586 //fepb
00587   int numFepInitial;  // no. of fep atoms with initial flag
00588   int numFepFinal;  // no. of fep atoms with final flag
00589 //fepe
00590 
00591   int numConsForce; //  Number of atoms that have constant force applied
00592   int32 *consForceIndexes;//  Constant force indexes for each atom
00593   Vector *consForce;  //  Constant force array
00594 
00595   int32 *consTorqueIndexes; //  "Constant" torque indexes for each atom
00596   ConsTorqueParams *consTorqueParams;
00597                                 //  Parameters for each atom "constant"-torqued
00598 
00599   // The following are needed for error checking because we
00600   // eliminate bonds, etc. which involve only fixed atoms
00601   int numCalcBonds; //  Number of bonds requiring calculation
00602   int numCalcAngles;  //  Number of angles requiring calculation
00603   int numCalcDihedrals; //  Number of dihedrals requiring calculation
00604   int numCalcImpropers; //  Number of impropers requiring calculation
00605   int numCalcCrossterms; //  Number of cross-terms requiring calculation
00606   int numCalcExclusions;  //  Number of exclusions requiring calculation
00607   int numCalcFullExclusions;  //  Number of full exclusions requiring calculation
00608 
00609   // DRUDE
00610   int numCalcTholes;  // Number of Thole correction terms requiring calculation
00611   int numCalcAnisos;  // Number of anisotropic terms requiring calculation
00612   // DRUDE
00613 
00614   //  Number of dihedrals with multiple periodicity
00615   int numMultipleDihedrals; 
00616   //  Number of impropers with multiple periodicity
00617   int numMultipleImpropers; 
00618   // indexes of "atoms" sorted by hydrogen groups
00619   HydrogenGroup hydrogenGroup;
00620 
00621   // Ported by JLai -- JE - Go
00622   int numGoAtoms;         //  Number of atoms subject to Go forces -- ported by JLai/ Original by JE
00623   int32 *atomChainTypes;  //  Go chain type for each atom; from 1 to MAX_GO_CHAINS
00624   int32 *goSigmaIndices;  //  Indices into goSigmas
00625   int32 *goResidIndices;  //  Indices into goSigmas
00626   Real  *goSigmas;        //  Sigma values for Go forces L-J type formula
00627   bool *goWithinCutoff;   //  Whether the reference atom-atom distance is within the Go cutoff
00628   Real  *goCoordinates;   //  Coordinates (x,y,z) for Go atoms in the native structure
00629   int *goResids;          //  Residue ID from PDB
00630   PDB *goPDB;             //  Pointer to PDB object to use
00631   // NAMD-Go2 calculation code
00632   int goNumLJPair;        //  Integer storing the total number of explicit pairs (LJ)
00633   int *goIndxLJA;         //  Pointer to the array of atom indices for LJ atom A
00634   int *goIndxLJB;         //  Pointer to the array of atom indices for LJ atom B
00635   double *goSigmaPairA;  //  Pointer to the array of A LJ parameters
00636   double *goSigmaPairB;  //  Pointer to the array of B LJ parameters
00637   int *pointerToGoBeg;    //  Array of pointers to array
00638   int *pointerToGoEnd;    //  Array of pointers to array
00639   // Gromacs LJ Pair list calculation code
00640   int numPair;            //  Integer storing the total number of explicit pairs (LJ + Gaussian)
00641   int numLJPair;          //  Integer storing the number of explicit LJ pairs
00642   int numCalcLJPair;         //  Number of explicit LJ pairs requiring calculation
00643   int *pointerToLJBeg;       //  Array of pointers to array 
00644   int *pointerToLJEnd;       //  Array of pointers to array B
00645   int *indxLJA;           //  Pointer to the array of atom indices for LJ atom A
00646   int *indxLJB;           //  Pointer to the array of atom indices for LJ atom B
00647   Real *pairC6;           //  Pointer to the array of C6 LJ parameters
00648   Real *pairC12;          //  Pointer to the array of C12 LJ parameters
00649   int *gromacsPair_type;   //  
00650   // Gromacs Gauss Pair list calculation code
00651   int *pointerToGaussBeg;    //  Array of pointers to array B
00652   int *pointerToGaussEnd;    //  Array of pointers to array B
00653   int numGaussPair;       //  Integer storing the number of explicit Gaussian pairs  
00654   int *indxGaussA;        //  Pointer to the array of atom indices for Gaussian atom A 
00655   int *indxGaussB;        //  Pointer to the array of atom indices for Gaussian atom B 
00656   Real *gA;               //  Pointer to the array of force constants to the Gaussian potential
00657   Real *gMu1;             //  Pointer to the array of mu (shifts Gaussian)
00658   Real *giSigma1;          //  Pointer to the array of sigma (controls spread of Gaussian)
00659   Real *gMu2;             //  Pointer to the array of mu (shifts Gaussian 2)
00660   Real *giSigma2;          //  Pointer to the array of sigma (controls spread of Gaussian 2)
00661   Real *gRepulsive;       //  Pointer to the a LJ-12 repulsive parameter that adds to the Gaussian
00662 
00663   // GO ENERGY CALCULATION CODE
00664   BigReal energyNative;    // Holds the energy value of the native structure
00665   BigReal energyNonnative; // Holds the energy value of the nonnative structure
00666   // GO ENERGY CALCULATION CODE
00667   // End of port - JL
00668 
00669   
00670 private:
00672   // Begin QM
00674     
00675   int qmDroppedBonds, qmDroppedAngles, qmDroppedDihedrals;
00676   int qmDroppedImpropers, qmDroppedCrossterms;
00677   ResizeArray<Bond> qmExtraBonds;
00678   
00679   Bool qmReplaceAll;              // Indicates if all forces should be replaced.
00680   int qmNumQMAtoms;           // Number of QM atoms, from all groups.
00681   
00682   // Array of abbreviated element names for all atoms.
00683   char **qmElementArray;
00684   // Array of elements for dummy atoms.
00685   char **qmDummyElement;
00686   
00687   // Number of QM atoms in each QM group
00688   int *qmGrpSizes;
00689   
00690   // QM group for each atom of the system. 0 means MM atom.
00691   Real *qmAtomGroup;
00692   // Number of different QM regions
00693   int qmNumGrps ;
00694   
00695   // QM Atom charges.
00696   Real *qmAtmChrg ;
00697   // global indices of all QM atoms.
00698   int *qmAtmIndx ;
00699   
00700   // Number of QM-MM bonds.
00701   int qmNumBonds ;
00702   // IDs of each group, that is, the value in the qmColumn, per group.
00703   Real *qmGrpID ;
00704   // The total charge of each QM group.
00705   Real *qmGrpChrg;
00706   // The multiplicity of QM groups
00707   Real *qmGrpMult;
00708   // Number of QM-MM bonds in each group.
00709   int *qmGrpNumBonds ;
00710   // Sequential list of bonds between a MM atom and a QM atom.
00711   // (with the bonded atoms ins this order: MM first, QM second).
00712   int **qmMMBond ;
00713   // List of bond indexes for each QM group.
00714   int **qmGrpBonds ;
00715   // List of atom indexes for MM atoms in QM-MM bonds, per group.
00716   // This is used to check if a point charge is actually an MM atom
00717   // that need to be ignored, and will be substituted by a dummy atom (for example).
00718   int **qmMMBondedIndx ;
00719   // List of indices for MM atoms which will receive fractions of charges from 
00720   // the MM atom of a QM-MM bond. This is used in the Charge Shift scheme.
00721   int **qmMMChargeTarget;
00722   // Number of target MM atoms per QM-MM bond. Targets will receive the partial
00723   // charge of the MM atom in a QM-MM bond. (in Charge shift scheme)
00724   int *qmMMNumTargs;
00725   // List of distances from thr QM atom where the dummy atom will be placed.
00726   BigReal *qmDummyBondVal;
00727   // Indicate if point charges will be used in QM calculations.
00728   Bool qmNoPC;
00729   
00731   // These variables are used in case mechanichal embedding is requested (NoPC = TRUE)
00732   // AND there are QM-MM bonds in the system.
00733   
00734   // Indicates the number of items in the arrays for Mechanical embedding with QM-MM bonds.
00735   int qmMeNumBonds;
00736   // Indicates the indices of MM atoms which participate in QM-MM bonds.
00737   int *qmMeMMindx;
00738   // Indicates the QM group to which the QM-MM bonds belongs.
00739   Real *qmMeQMGrp;
00741   
00742   // Indicates frequency of selection of point charges.
00743   int qmPCFreq;
00744   // Custom PC array;
00745   int *qmCustomPCIdxs;
00746   // Number of Indexes associated with each QM Group.
00747   int *qmCustPCSizes;
00748   // Total number of custom PC charges.
00749   int qmTotCustPCs;
00750   // Number of solvent molecules that will be selected and updated by NAMD, per group.
00751   int *qmLSSSize;
00752   // Number of atoms per solvent molecule. We need all molecules to have the same size.
00753   int qmLSSResidueSize;
00754   // IDs of all atoms belonging to QM solvent molecules, from all groups.
00755   int *qmLSSIdxs;
00756   // MAss of each atom in LSS solvent molecules.
00757   Mass *qmLSSMass;
00758   // Atom IDs of QM atoms which will be used to select solvent molecules by distance.
00759   int *qmLSSRefIDs;
00760   // Mass of each QM atom used as reference for solvent selection.
00761   Mass *qmLSSRefMass ;
00762   // Number of atom IDs/Mass per group.
00763   int *qmLSSRefSize;
00764   int qmLSSFreq;
00765   int qmLSSTotalNumAtms;
00766   // This will map each classical solvent atom with a global index of solvent residues,
00767   // so we don't depend on segment names and repeated residue indices.
00768   std::map<int,int> qmClassicSolv ;
00769   
00771   // Conditional SMD
00772     
00773     void read_qm_csdm_file(std::map<Real,int> &qmGrpIDMap) ;
00774     
00775     Bool qmcSMD;
00776         // Total number of cSMD instances.
00777         int cSMDnumInst;
00778     // Num instances per QM group
00779     int* cSMDindxLen;
00780     // Index of Conditional SMD guides per QM group
00781     int** cSMDindex;
00782     // Atom indices for Origin and Target of cSMD
00783     int** cSMDpairs;
00784     // Spring constants for cSMD
00785     Real* cSMDKs;
00786     // Speed of movement of guide particles for cSMD.
00787     Real* cSMDVels;
00788     // Distance cutoff for guide particles for cSMD.
00789     Real* cSMDcoffs;
00790   
00792   // end QM
00794   
00795 public:
00796   // QM
00797   void set_qm_replaceAll(Bool newReplaceAll) { qmReplaceAll = newReplaceAll; } ;
00798   
00799   const Real * get_qmAtomGroup() const {return qmAtomGroup; } ;
00800   Real get_qmAtomGroup(int indx) const {return qmAtomGroup[indx]; } ;
00801   
00802   Real *get_qmAtmChrg() {return qmAtmChrg; } ;
00803   const int *get_qmAtmIndx() {return qmAtmIndx; } ;
00804   
00805   int get_numQMAtoms() {return qmNumQMAtoms; } ;
00806   const char *const * get_qmElements() {return qmElementArray;} ;
00807   int get_qmNumGrps() const {return qmNumGrps; } ;
00808   const int * get_qmGrpSizes() {return qmGrpSizes; } ;
00809   Real * get_qmGrpID() { return qmGrpID; } ;
00810   Real *get_qmGrpChrg() {return qmGrpChrg; } ;
00811   Real *get_qmGrpMult() {return qmGrpMult; } ;
00812   
00813   int get_qmNumBonds() { return qmNumBonds; } ;
00814   const BigReal * get_qmDummyBondVal() { return qmDummyBondVal; } ;
00815   const int * get_qmGrpNumBonds() { return qmGrpNumBonds; } ;
00816   const int *const * get_qmMMBond() { return qmMMBond; } ;
00817   const int *const * get_qmGrpBonds() { return qmGrpBonds; } ;
00818   const int *const * get_qmMMBondedIndx() { return qmMMBondedIndx; } ;
00819   const char *const * get_qmDummyElement() { return qmDummyElement; } ;
00820   
00821   const int *const * get_qmMMChargeTarget() { return qmMMChargeTarget; } ;
00822   const int * get_qmMMNumTargs() { return qmMMNumTargs; } ;
00823   
00824   int* get_qmLSSSize() { return qmLSSSize;} ;
00825   int* get_qmLSSIdxs() { return qmLSSIdxs;} ;
00826   Mass *get_qmLSSMass() { return qmLSSMass; } ;
00827   int get_qmLSSFreq() { return qmLSSFreq; } ;
00828   int get_qmLSSResSize() { return qmLSSResidueSize; } ;
00829   int *get_qmLSSRefIDs() { return qmLSSRefIDs ; } ;
00830   int *get_qmLSSRefSize() { return qmLSSRefSize ; } ;
00831   Mass *get_qmLSSRefMass() {return qmLSSRefMass; } ;
00832   std::map<int,int> &get_qmMMSolv() { return qmClassicSolv;} ;
00833   
00834   Bool get_qmReplaceAll() {return qmReplaceAll; } ;
00835   
00836   Bool get_noPC() { return qmNoPC; } ;
00837   int get_qmMeNumBonds() { return qmMeNumBonds; } ;
00838   int *get_qmMeMMindx() { return qmMeMMindx; } ;
00839   Real *get_qmMeQMGrp() { return qmMeQMGrp; } ;
00840   
00841   int get_qmPCFreq() { return qmPCFreq; } ;
00842   
00843   int get_qmTotCustPCs() { return qmTotCustPCs; } ;
00844   int *get_qmCustPCSizes()  { return qmCustPCSizes; } ;
00845   int *get_qmCustomPCIdxs() { return qmCustomPCIdxs; } ;
00846   
00847   Bool get_qmcSMD() { return qmcSMD;} ;
00848   int get_cSMDnumInst() { return cSMDnumInst;} ;
00849   int const * get_cSMDindxLen() { return cSMDindxLen;} ;
00850   int const * const * get_cSMDindex() {return cSMDindex;} ;
00851   int const * const * get_cSMDpairs() {return cSMDpairs;} ;
00852   const Real* get_cSMDKs() {return cSMDKs;} ;
00853   const Real* get_cSMDVels() {return cSMDVels;} ;
00854   const Real* get_cSMDcoffs() {return cSMDcoffs;} ;
00855   
00856   void prepare_qm(const char *pdbFileName, Parameters *params, ConfigList *cfgList) ;
00857   void delete_qm_bonded() ;
00858   // end QM
00859   
00860   
00861   
00862   Molecule(SimParameters *, Parameters *param);
00863   Molecule(SimParameters *, Parameters *param, char *filename, ConfigList *cfgList=NULL);  
00864   Molecule(SimParameters *simParams, Parameters *param, molfile_plugin_t *pIOHdl, void *pIOFileHdl, int natoms);
00865   
00866   Molecule(SimParameters *, Parameters *, Ambertoppar *);
00867   void read_parm(Ambertoppar *);
00868 
00869   Molecule(SimParameters *, Parameters *, const GromacsTopFile *);
00870 
00871   ~Molecule();    //  Destructor
00872 
00873   void send_Molecule(MOStream *);
00874         //  send the molecular structure 
00875         //  from the master to the clients
00876 
00877   void receive_Molecule(MIStream *);
00878         //  receive the molecular structure
00879         //  from the master on a client
00880   
00881   void build_constraint_params(StringList *, StringList *, StringList *,
00882              PDB *, char *);
00883         //  Build the set of harmonic constraint 
00884         // parameters
00885 
00886 /* BEGIN gf */
00887   void build_gridforce_params(StringList *, StringList *, StringList *, StringList *, PDB *, char *);
00888         //  Build the set of gridForce-style force pars
00889 /* END gf */
00890 
00891   void build_movdrag_params(StringList *, StringList *, StringList *, 
00892           PDB *, char *);
00893         //  Build the set of moving drag pars
00894 
00895   void build_rotdrag_params(StringList *, StringList *, StringList *,
00896           StringList *, StringList *, StringList *,
00897           PDB *, char *);
00898         //  Build the set of rotating drag pars
00899 
00900   void build_constorque_params(StringList *, StringList *, StringList *,
00901              StringList *, StringList *, StringList *,
00902              PDB *, char *);
00903         //  Build the set of "constant" torque pars
00904 
00905 
00906   void build_constant_forces(char *);
00907         //  Build the set of constant forces
00908 
00909   void build_langevin_params(BigReal coupling, BigReal drudeCoupling,
00910       Bool doHydrogen);
00911   void build_langevin_params(StringList *, StringList *, PDB *, char *);
00912         //  Build the set of langevin dynamics parameters
00913 
00914 #ifdef MEM_OPT_VERSION
00915   void load_fixed_atoms(StringList *fixedFile);
00916   void load_constrained_atoms(StringList *constrainedFile);
00917 #endif
00918 
00919   void build_fixed_atoms(StringList *, StringList *, PDB *, char *);
00920         //  Determine which atoms are fixed (if any)
00921 
00922   void build_stirred_atoms(StringList *, StringList *, PDB *, char *);
00923         //  Determine which atoms are stirred (if any)
00924 
00925   void build_extra_bonds(Parameters *parameters, StringList *file);
00926 
00927 //fepb
00928         void build_fep_flags(StringList *, StringList *, PDB *, char *, const char *);
00929                                // selection of the mutant atoms
00930         void delete_alch_bonded(void);
00931 //fepe
00932 
00933   void build_exPressure_atoms(StringList *, StringList *, PDB *, char *);
00934         //  Determine which atoms are excluded from
00935                                 //  pressure (if any)
00936 
00937   // Ported by JLai -- Original JE - Go -- Change the unsigned int to ints
00938   void print_go_sigmas(); //  Print out Go sigma parameters
00939   void build_go_sigmas(StringList *, char *);
00940         //  Determine which atoms have Go forces applied
00941         //  calculate sigmas from distances between Go atom pairs
00942   void build_go_sigmas2(StringList *, char *);
00943         //  Determine which atoms have Go forces applied
00944         //  calculate sigmas from distances between Go atom pairs
00945   void build_go_arrays(StringList *, char *);
00946         //  Determine which atoms have Go forces applied
00947   BigReal get_gro_force(BigReal, BigReal, BigReal, int, int) const;
00948   BigReal get_gro_force2(BigReal, BigReal, BigReal, int, int, BigReal *, BigReal *) const;
00949   BigReal get_go_force(BigReal, int, int, BigReal *, BigReal *) const;
00950         //  Calculate the go force between a pair of atoms -- Modified to 
00951         //  output Go energies
00952   BigReal get_go_force_new(BigReal, int, int, BigReal *, BigReal *) const;
00953         //  Calculate the go force between a pair of atoms
00954   BigReal get_go_force2(BigReal, BigReal, BigReal, int, int, BigReal *, BigReal *) const;
00955         //  Calculate the go force between a pair of atoms
00956   Bool atoms_1to4(unsigned int, unsigned int);
00957 // End of port -- JL  
00958 
00959   void reloadCharges(float charge[], int n);
00960 
00961         Bool is_lp(int);     // return true if atom is a lone pair
00962         Bool is_drude(int);     // return true if atom is a Drude particle
00963         Bool is_hydrogen(int);     // return true if atom is hydrogen
00964         Bool is_oxygen(int);       // return true if atom is oxygen
00965   Bool is_hydrogenGroupParent(int); // return true if atom is group parent
00966   Bool is_water(int);        // return true if atom is part of water 
00967   int  get_groupSize(int);     // return # atoms in (hydrogen) group
00968         int get_mother_atom(int);  // return mother atom of a hydrogen
00969 
00970   #ifdef MEM_OPT_VERSION
00971   //the way to get the cluster size if the atom ids of the cluster are
00972   //contiguous. The input parameter is the atom id that leads the cluster.
00973   int get_cluster_size_con(int aid) const { return clusterSigs[aid]; }  
00974   //the way to get the cluster size if the atoms ids of the cluster are
00975   //not contiguous. The input parameter is the cluster index.
00976   int get_cluster_size_uncon(int cIdx) const { return clusterSize[cIdx]; }
00977   int get_cluster_idx(int aid) const { return clusterSigs[aid]; }
00978   int get_num_clusters() const { return numClusters; }
00979   #else
00980   int get_cluster(int anum) const { return cluster[anum]; }
00981   int get_clusterSize(int anum) const { return clusterSize[anum]; }
00982   #endif
00983 
00984 #ifndef MEM_OPT_VERSION
00985   const float *getOccupancyData() { return (const float *)occupancy; }
00986   void setOccupancyData(molfile_atom_t *atomarray);
00987   void freeOccupancyData() { delete [] occupancy; occupancy=NULL; }
00988 
00989   const float *getBFactorData() { return (const float *)bfactor; }
00990   void setBFactorData(molfile_atom_t *atomarray);
00991   void freeBFactorData() { delete [] bfactor; bfactor=NULL; }
00992 #endif
00993 
00994   //  Get the mass of an atom
00995   Real atommass(int anum) const
00996   {
00997     #ifdef MEM_OPT_VERSION
00998     return atomMassPool[eachAtomMass[anum]];
00999     #else
01000     return(atoms[anum].mass);
01001     #endif
01002   }
01003 
01004   //  Get the charge of an atom
01005   Real atomcharge(int anum) const
01006   {
01007     #ifdef MEM_OPT_VERSION
01008     return atomChargePool[eachAtomCharge[anum]];
01009     #else
01010     return(atoms[anum].charge);
01011     #endif
01012   }
01013   
01014   //  Get the vdw type of an atom
01015   Index atomvdwtype(int anum) const
01016   {      
01017       return(atoms[anum].vdw_type);
01018   }
01019 
01020   #ifndef MEM_OPT_VERSION
01021   //  Retrieve a bond structure
01022   Bond *get_bond(int bnum) const {return (&(bonds[bnum]));}
01023 
01024   //  Retrieve an angle structure
01025   Angle *get_angle(int anum) const {return (&(angles[anum]));}
01026 
01027   //  Retrieve an improper strutcure
01028   Improper *get_improper(int inum) const {return (&(impropers[inum]));}
01029 
01030   //  Retrieve a dihedral structure
01031   Dihedral *get_dihedral(int dnum) const {return (&(dihedrals[dnum]));}
01032 
01033   //  Retrieve a cross-term strutcure
01034   Crossterm *get_crossterm(int inum) const {return (&(crossterms[inum]));}
01035   #endif
01036 
01037   // DRUDE: retrieve lphost structure
01038   Lphost *get_lphost(int atomid) const {
01039     // don't call unless simParams->drudeOn == TRUE
01040     // otherwise lphostIndexes array doesn't exist!
01041     int index = lphostIndexes[atomid];
01042     return (index != -1 ? &(lphosts[index]) : NULL);
01043   }
01044   // DRUDE
01045 
01046   #ifndef MEM_OPT_VERSION
01047   Bond *getAllBonds() const {return bonds;}
01048   Angle *getAllAngles() const {return angles;}
01049   Improper *getAllImpropers() const {return impropers;}
01050   Dihedral *getAllDihedrals() const {return dihedrals;}
01051   Crossterm *getAllCrossterms() const {return crossterms;}
01052   #endif
01053 
01054   // DRUDE: retrieve entire lphosts array
01055   Lphost *getAllLphosts() const { return lphosts; }
01056   // DRUDE
01057 
01058   //  Retrieve a hydrogen bond donor structure
01059   Bond *get_donor(int dnum) const {return (&(donors[dnum]));}  
01060 
01061   //  Retrieve a hydrogen bond acceptor structure
01062   Bond *get_acceptor(int dnum) const {return (&(acceptors[dnum]));} 
01063 
01064   Bond *getAllDonors() const {return donors;}
01065   Bond *getAllAcceptors() const {return acceptors;}
01066 
01067   //  Retrieve an exclusion structure
01068   #ifndef MEM_OPT_VERSION
01069   Exclusion *get_exclusion(int ex) const {return (&(exclusions[ex]));}
01070   #endif
01071 
01072   //  Retrieve an atom type
01073   const char *get_atomtype(int anum) const
01074   {
01075     if (atomNames == NULL)
01076     {
01077       NAMD_die("Tried to find atom type on node other than node 0");
01078     }
01079 
01080     #ifdef MEM_OPT_VERSION    
01081     return atomTypePool[atomNames[anum].atomtypeIdx];
01082     #else
01083     return(atomNames[anum].atomtype);
01084     #endif
01085   }
01086 
01087   //  Lookup atom id from segment, residue, and name
01088   int get_atom_from_name(const char *segid, int resid, const char *aname) const;
01089 
01090   //  Lookup number of atoms in residue from segment and residue
01091   int get_residue_size(const char *segid, int resid) const;
01092 
01093   //  Lookup atom id from segment, residue, and index in residue
01094   int get_atom_from_index_in_residue(const char *segid, int resid, int index) const;
01095 
01096   
01097   //  The following routines are used to get the list of bonds
01098   //  for a given atom.  This is used when creating the bond lists
01099   //  for the force objects  
01100 
01101   #ifndef MEM_OPT_VERSION
01102   int32 *get_bonds_for_atom(int anum)
01103       { return bondsByAtom[anum]; } 
01104   int32 *get_angles_for_atom(int anum) 
01105       { return anglesByAtom[anum]; }
01106   int32 *get_dihedrals_for_atom(int anum) 
01107       { return dihedralsByAtom[anum]; }
01108   int32 *get_impropers_for_atom(int anum) 
01109       { return impropersByAtom[anum]; }  
01110   int32 *get_crossterms_for_atom(int anum) 
01111       { return crosstermsByAtom[anum]; }  
01112   int32 *get_exclusions_for_atom(int anum)
01113       { return exclusionsByAtom[anum]; }
01114   const int32 *get_full_exclusions_for_atom(int anum) const
01115       { return fullExclusionsByAtom[anum]; }
01116   const int32 *get_mod_exclusions_for_atom(int anum) const
01117       { return modExclusionsByAtom[anum]; }
01118   #endif
01119   
01120   //  Check for exclusions, either explicit or bonded.
01121         //  Returns 1 for full, 2 for 1-4 exclusions.
01122   #ifdef MEM_OPT_VERSION
01123   int checkExclByIdx(int idx1, int atom1, int atom2) const;
01124   const ExclusionCheck *get_excl_check_for_idx(int idx) const{      
01125       return &exclChkSigPool[idx];
01126   }
01127   #else
01128   int checkexcl(int atom1, int atom2) const;
01129 
01130   const ExclusionCheck *get_excl_check_for_atom(int anum) const{      
01131       return &all_exclusions[anum];             
01132   }
01133   #endif
01134 
01135 /* BEGIN gf */
01136   // Return true or false based on whether or not the atom
01137   // is subject to grid force
01138   Bool is_atom_gridforced(int atomnum, int gridnum) const
01139   {
01140       if (numGridforceGrids)
01141       {
01142           return(gridfrcIndexes[gridnum][atomnum] != -1);
01143       }
01144       else
01145       {
01146           return(FALSE);
01147       }
01148   }
01149 /* END gf */
01150 
01151 #ifndef MEM_OPT_VERSION
01152   //  Return true or false based on whether the specified atom
01153   //  is constrained or not.
01154   Bool is_atom_constrained(int atomnum) const
01155   {
01156     if (numConstraints)
01157     {
01158       //  Check the index to see if it is constrained
01159       return(consIndexes[atomnum] != -1);
01160     }
01161     else
01162     {
01163       //  No constraints at all, so just return FALSE
01164       return(FALSE);
01165     }
01166   }
01167 #endif
01168 
01169   //  Return true or false based on whether the specified atom
01170   //  is moving-dragged or not.
01171   Bool is_atom_movdragged(int atomnum) const
01172   {
01173     if (numMovDrag)
01174     {
01175       //  Check the index to see if it is constrained
01176       return(movDragIndexes[atomnum] != -1);
01177     }
01178     else
01179     {
01180       //  No constraints at all, so just return FALSE
01181       return(FALSE);
01182     }
01183   }
01184 
01185   //  Return true or false based on whether the specified atom
01186   //  is rotating-dragged or not.
01187   Bool is_atom_rotdragged(int atomnum) const
01188   {
01189     if (numRotDrag)
01190     {
01191       //  Check the index to see if it is constrained
01192       return(rotDragIndexes[atomnum] != -1);
01193     }
01194     else
01195     {
01196       //  No constraints at all, so just return FALSE
01197       return(FALSE);
01198     }
01199   }
01200 
01201   //  Return true or false based on whether the specified atom
01202   //  is "constant"-torqued or not.
01203   Bool is_atom_constorqued(int atomnum) const
01204   {
01205     if (numConsTorque)
01206     {
01207       //  Check the index to see if it is constrained
01208       return(consTorqueIndexes[atomnum] != -1);
01209     }
01210     else
01211     {
01212       //  No constraints at all, so just return FALSE
01213       return(FALSE);
01214     }
01215   }
01216 
01217 #ifndef MEM_OPT_VERSION
01218   //  Get the harmonic constraints for a specific atom
01219   void get_cons_params(Real &k, Vector &refPos, int atomnum) const
01220   {
01221     k = consParams[consIndexes[atomnum]].k;
01222     refPos = consParams[consIndexes[atomnum]].refPos;
01223   }
01224 #endif
01225 
01226 /* BEGIN gf */
01227   void get_gridfrc_params(Real &k, Charge &q, int atomnum, int gridnum) const
01228   {
01229       k = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].k;
01230       q = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].q;
01231   }
01232   
01233   GridforceGrid* get_gridfrc_grid(int gridnum) const
01234   {
01235       GridforceGrid *result = NULL;
01236       if (gridnum >= 0 && gridnum < numGridforceGrids) {
01237           result = gridfrcGrid[gridnum];
01238       }
01239       return result;
01240   }
01241   
01242   int set_gridfrc_grid(int gridnum, GridforceGrid *grid)
01243   {
01244       if (grid && gridnum >= 0 && gridnum < numGridforceGrids) {
01245           gridfrcGrid[gridnum] = grid;
01246           return 0;
01247       } else {
01248           return -1;
01249       }
01250   }
01251 /* END gf */
01252 
01253   Real langevin_param(int atomnum) const
01254   {
01255     return(langevinParams ? langevinParams[atomnum] : 0.);
01256   }
01257 
01258   //  Get the stirring constraints for a specific atom
01259   void get_stir_refPos(Vector &refPos, int atomnum) const
01260   {
01261     refPos = stirParams[stirIndexes[atomnum]].refPos;
01262   }
01263 
01264 
01265   void put_stir_startTheta(Real theta, int atomnum) const
01266   {
01267     stirParams[stirIndexes[atomnum]].startTheta = theta;
01268   }
01269 
01270 
01271   Real get_stir_startTheta(int atomnum) const
01272   {
01273     return stirParams[stirIndexes[atomnum]].startTheta;
01274   }
01275  
01276 
01277   //  Get the moving drag factor for a specific atom
01278   void get_movdrag_params(Vector &v, int atomnum) const
01279   {
01280     v = movDragParams[movDragIndexes[atomnum]].v;
01281   }
01282 
01283   //  Get the rotating drag pars for a specific atom
01284   void get_rotdrag_params(BigReal &v, Vector &a, Vector &p, 
01285         int atomnum) const
01286   {
01287     v = rotDragParams[rotDragIndexes[atomnum]].v;
01288     a = rotDragParams[rotDragIndexes[atomnum]].a;
01289     p = rotDragParams[rotDragIndexes[atomnum]].p;
01290   }
01291 
01292   //  Get the "constant" torque pars for a specific atom
01293   void get_constorque_params(BigReal &v, Vector &a, Vector &p, 
01294         int atomnum) const
01295   {
01296     v = consTorqueParams[consTorqueIndexes[atomnum]].v;
01297     a = consTorqueParams[consTorqueIndexes[atomnum]].a;
01298     p = consTorqueParams[consTorqueIndexes[atomnum]].p;
01299   }
01300 
01301 //fepb
01302   unsigned char get_fep_type(int anum) const
01303   {
01304     return(fepAtomFlags[anum]);
01305   }
01306 
01307   /* BKR - Get the FEP type (i.e. 0, 1, or 2) of a bonded _interaction_ based 
01308      on the atom indices of the atoms involved (internally converted to FEP 
01309      types) and the order of the bonded interaction (i.e. 2, 3, or 4). This 
01310      automatically accounts for whether or not purely alchemical interactions
01311      are being scaled (e.g. bonds between two atoms of fep_type 1).
01312 
01313      The logic here is admittedly a bit opaque. When the fep_type's are back
01314      mapped to -1,0,1, we can use the sum of all of the types to determine the 
01315      type of interaction for all bonded term types:
01316 
01317      0  - only non-alchemical atoms involved
01318      >0 - _atleast_ one appearing atom
01319      <0 - _atleast_ one disappearing atom
01320 
01321      If the magnitude of the sum is equal to the interaction order (i.e. 2, 3, 
01322      or 4), then it is a _purely_ alchemical interaction and it may be 
01323      desireable to retain it for sampling purposes (note that this adds a 
01324      constant to the free energy that will almost always differ at the 
01325      endpoints).  In order to avoid unexpected instability these interactions 
01326      are retained by default, but can be scaled along with everything else by 
01327      setting "alchBondDecouple on".
01328 
01329      NB: All pure alchemical interactions beyond order 2 are ALWAYS discarded
01330      by alchemify. This is good, because higher order interactions would break
01331      the above logic. For example, if we had an angle term between atoms of 
01332      types (1,1,-1) the sum would be 1, but this term should receive no scaling
01333      because it involves groups -1 and 1 but not 0.
01334   */
01335   int get_fep_bonded_type(const int *atomID, unsigned int order) const
01336   {
01337     int typeSum = 0;
01338     for ( int i=0; i < order; ++i ) {
01339       typeSum += (fepAtomFlags[atomID[i]] == 2 ? -1 : fepAtomFlags[atomID[i]]);
01340     }
01341     // Increase the cutoff if scaling purely alchemical bonds. 
01342     // This really only applies when order = 2.
01343     if ( simParams->alchBondDecouple ) order++;
01344 
01345     if ( typeSum == 0 ) return 0; // Most interactions get caught here.
01346     else if ( 0 < typeSum && typeSum < order ) return 1;
01347     else if ( 0 > typeSum && typeSum > -order ) return 2;
01348 
01349     if ( simParams->alchBondDecouple ) {
01350       // Alchemify should always keep this from bombing, but just in case...
01351       NAMD_die("Unexpected alchemical bonded interaction!");
01352     }
01353   }
01354 //fepe
01355 
01356 #ifndef MEM_OPT_VERSION
01357   Bool is_atom_fixed(int atomnum) const
01358   {
01359     return (numFixedAtoms && fixedAtomFlags[atomnum]);
01360   }
01361 #else
01362   //Since binary search is more expensive than direct array access,
01363   //and this function is usually called for consecutive atoms in this context,
01364   //the *listIdx returns the index to the range of atoms [aid1, aid2]
01365   //that are fixed. If the atom aid is fixed, then aid1=<aid<=aid2;
01366   //If the atom aid is not fixed, then aid1 indicates the smallest fixed atom
01367   //id that is larger than aid; so the listIdx could be equal the size of
01368   //fixedAtomsSet. --Chao Mei
01369   Bool is_atom_in_set(AtomSetList *localAtomsSet, int aid, int *listIdx) const;
01370   inline Bool is_atom_fixed(int aid, int *listIdx=NULL) const {
01371     return is_atom_in_set(fixedAtomsSet,aid,listIdx);
01372   }
01373   inline Bool is_atom_constrained(int aid, int *listIdx=NULL) const {
01374     return is_atom_in_set(constrainedAtomsSet,aid,listIdx);
01375   }
01376 #endif
01377         
01378   Bool is_atom_stirred(int atomnum) const
01379   {
01380     if (numStirredAtoms)
01381     {
01382       //  Check the index to see if it is constrained
01383       return(stirIndexes[atomnum] != -1);
01384     }
01385     else
01386     {
01387       //  No constraints at all, so just return FALSE
01388       return(FALSE);
01389     }
01390   }
01391   
01392 
01393   Bool is_group_fixed(int atomnum) const
01394   {
01395     return (numFixedAtoms && (fixedAtomFlags[atomnum] == -1));
01396   }
01397   Bool is_atom_exPressure(int atomnum) const
01398   {
01399     return (numExPressureAtoms && exPressureAtomFlags[atomnum]);
01400   }
01401   // 0 if not rigid or length to parent, for parent refers to H-H length
01402   // < 0 implies MOLLY but not SHAKE, > 1 implies both if MOLLY is on
01403   Real rigid_bond_length(int atomnum) const
01404   {
01405     return(rigidBondLengths[atomnum]);
01406   }
01407 
01408   void print_atoms(Parameters *); 
01409         //  Print out list of atoms
01410   void print_bonds(Parameters *); 
01411         //  Print out list of bonds
01412   void print_exclusions();//  Print out list of exclusions
01413 
01414 public:  
01415   int isOccupancyValid, isBFactorValid;
01416 
01417 #ifdef MEM_OPT_VERSION
01418   //read the per-atom file for the memory optimized version where the file 
01419   //name already exists the range of atoms to be read are 
01420   //[fromAtomID, toAtomID].
01421   void read_binary_atom_info(int fromAtomID, int toAtomID, InputAtomList &inAtoms);
01422 
01423   int getNumCalcExclusions(){return numCalcExclusions;}
01424   void setNumCalcExclusions(int x){numCalcExclusions= x;}
01425 
01426   Index getEachAtomMass(int i){return eachAtomMass[i];}
01427   Index getEachAtomCharge(int i){return eachAtomCharge[i];}
01428 
01429   ExclSigID getAtomExclSigId(int aid) const {
01430       return eachAtomExclSig[aid];
01431   }
01432 
01433   Real *getAtomMassPool(){return atomMassPool;}
01434   Real *getAtomChargePool(){return atomChargePool;}
01435   AtomCstInfo *getAtoms(){return atoms;}  
01436 
01437   int atomSigPoolSize;
01438   AtomSignature *atomSigPool;
01439 
01440   /* All the following are temporary variables for reading the compressed psf file */
01441   //declarations for atoms' constant information  
01442   int segNamePoolSize; //Its value is usually less than 5
01443   char **segNamePool; //This seems not to be important, but it only occupied very little space.
01444 
01445   int resNamePoolSize;
01446   char **resNamePool;
01447 
01448   int atomNamePoolSize;
01449   char **atomNamePool;
01450 
01451   int atomTypePoolSize;
01452   char **atomTypePool;
01453 
01454   int chargePoolSize;
01455   Real *atomChargePool;
01456 
01457   int massPoolSize;
01458   Real *atomMassPool;
01459 
01460   AtomSigID getAtomSigId(int aid) {
01461       return eachAtomSig[aid]; 
01462   }
01463 
01464   //Indicates the size of both exclSigPool and exclChkSigPool
01465   int exclSigPoolSize;
01466   //this will be deleted after build_lists_by_atom
01467   ExclusionSignature *exclSigPool;
01468   //This is the final data structure we want to store  
01469   ExclusionCheck *exclChkSigPool;
01470 
01471   void addNewExclSigPool(const std::vector<ExclusionSignature>&);  
01472 
01473   void delEachAtomSigs(){    
01474       //for NAMD-smp version, only one Molecule object is held
01475       //on each node, therefore, only one deletion operation should
01476       //be taken on a node, otherwise, there possibly would be some
01477       //wierd memory problems. The same reason applies to other deletion
01478       //operations inside the Molecule object.   
01479       if(CmiMyRank()) return;
01480 
01481       delete [] eachAtomSig;
01482       delete [] eachAtomExclSig;
01483       eachAtomSig = NULL;
01484       eachAtomExclSig = NULL;
01485   }
01486 
01487   void delChargeSpace(){
01488       if(CmiMyRank()) return;
01489 
01490       delete [] atomChargePool;
01491       delete [] eachAtomCharge;
01492       atomChargePool = NULL;
01493       eachAtomCharge = NULL;
01494   }
01495   
01496   void delMassSpace(){
01497       if(CmiMyRank()) return;
01498 
01499       delete [] atomMassPool;
01500       delete [] eachAtomMass;
01501       atomMassPool = NULL;
01502       eachAtomMass = NULL;
01503   }
01504   
01505   void delClusterSigs() {
01506       if(CmiMyRank()) return;      
01507 
01508       delete [] clusterSigs;
01509       clusterSigs = NULL;
01510   }
01511 
01512   void delAtomNames(){
01513       if(CmiMyRank()) return;
01514       delete [] atomNamePool;
01515       delete [] atomNames;
01516       atomNamePool = NULL;
01517       atomNames = NULL;
01518   }
01519 
01520   void delFixedAtoms(){
01521       if(CmiMyRank()) return;
01522       delete fixedAtomsSet;
01523       fixedAtomsSet = NULL;
01524   }
01525 
01526 private:
01527   Index insert_new_mass(Real newMass);
01528 
01529 #endif
01530 
01531 // Go stuff
01532 public:
01533 
01534 GoValue go_array[MAX_GO_CHAINS*MAX_GO_CHAINS];    //  Array of Go params -- JLai
01535 int go_indices[MAX_GO_CHAINS+1];        //  Indices from chainIDS to go array -- JLai
01536 int NumGoChains;                        //  Number of Go chain types -- JLai
01537 
01538 // Declares and initializes Go variables
01539 void goInit();
01540 
01541 // Builds explicit Gromacs pairs
01542 void build_gro_pair();
01543 
01544 // Builds the initial Go parameters 
01545 void build_go_params(StringList *);
01546 
01547 //  Read Go parameter file
01548 void read_go_file(char *);
01549 
01550 //  Get Go cutoff for a given chain type pair
01551 Real get_go_cutoff(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].cutoff; };
01552 
01553 //  Get Go epsilonRep for a given chain type pair
01554 Real get_go_epsilonRep(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilonRep; };
01555 
01556 //  Get Go sigmaRep for a given chain type pair
01557 Real get_go_sigmaRep(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].sigmaRep; };
01558 
01559 //  Get Go epsilon for a given chain type pair
01560 Real get_go_epsilon(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilon; };
01561 
01562 //  Get Go exp_a for a given chain type pair
01563 int get_go_exp_a(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_a; };
01564 
01565 //  Get Go exp_b for a given chain type pair
01566 int get_go_exp_b(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_b; };
01567 
01568 //  Get Go exp_rep for a given chain type pair
01569 int get_go_exp_rep(int chain1, int chain2) { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_rep; };
01570 
01571 //  Whether residue IDs with this difference are restricted
01572 Bool go_restricted(int, int, int);
01573 
01574 // Prints Go Params
01575 void print_go_params();
01576 
01577 void initialize();
01578 
01579 void send_GoMolecule(MOStream *);
01580 //  send the molecular structure 
01581 //  from the master to the clients
01582 
01583 void receive_GoMolecule(MIStream *);
01584 //  receive the molecular structure
01585 //  from the master on a client
01586 };
01587 
01588 #endif
01589 

Generated on Thu Sep 21 01:17:13 2017 for NAMD by  doxygen 1.4.7