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

Generated on Sat Jul 21 01:17:14 2018 for NAMD by  doxygen 1.4.7