Molecule Class Reference

#include <Molecule.h>

List of all members.

Classes

struct  constorque_params
struct  constraint_params
struct  gridfrc_params
struct  movdrag_params
struct  rotdrag_params
struct  stir_params

Public Member Functions

void compute_LJcorrection ()
BigReal getEnergyTailCorr (const BigReal, const int)
BigReal getVirialTailCorr (const BigReal)
int const * getLcpoParamType ()
BigReal GetAtomAlpha (int i) const
AtomgetAtoms () const
AtomNameInfogetAtomNames () const
AtomSegResInfogetAtomSegResInfo () const
int num_fixed_atoms () const
int num_fixed_groups () const
int64_t num_group_deg_freedom () const
int64_t num_deg_freedom (int isInitialReport=0) const
void set_qm_replaceAll (Bool newReplaceAll)
const Realget_qmAtomGroup () const
Real get_qmAtomGroup (int indx) const
Realget_qmAtmChrg ()
const int * get_qmAtmIndx ()
int get_numQMAtoms ()
const char *const * get_qmElements ()
int get_qmNumGrps () const
const int * get_qmGrpSizes ()
Realget_qmGrpID ()
Realget_qmGrpChrg ()
Realget_qmGrpMult ()
int get_qmNumBonds ()
const BigRealget_qmDummyBondVal ()
const int * get_qmGrpNumBonds ()
const int *const * get_qmMMBond ()
const int *const * get_qmGrpBonds ()
const int *const * get_qmMMBondedIndx ()
const char *const * get_qmDummyElement ()
const int *const * get_qmMMChargeTarget ()
const int * get_qmMMNumTargs ()
int * get_qmLSSSize ()
int * get_qmLSSIdxs ()
Massget_qmLSSMass ()
int get_qmLSSFreq ()
int get_qmLSSResSize ()
int * get_qmLSSRefIDs ()
int * get_qmLSSRefSize ()
Massget_qmLSSRefMass ()
std::map< int, int > & get_qmMMSolv ()
Bool get_qmReplaceAll ()
Bool get_noPC ()
int get_qmMeNumBonds ()
int * get_qmMeMMindx ()
Realget_qmMeQMGrp ()
int get_qmPCFreq ()
int get_qmTotCustPCs ()
int * get_qmCustPCSizes ()
int * get_qmCustomPCIdxs ()
Bool get_qmcSMD ()
int get_cSMDnumInst ()
int const * get_cSMDindxLen ()
int const *const * get_cSMDindex ()
int const *const * get_cSMDpairs ()
const Realget_cSMDKs ()
const Realget_cSMDVels ()
const Realget_cSMDcoffs ()
void prepare_qm (const char *pdbFileName, Parameters *params, ConfigList *cfgList)
void delete_qm_bonded ()
 Molecule (SimParameters *, Parameters *param)
 Molecule (SimParameters *, Parameters *param, char *filename, ConfigList *cfgList=NULL)
 Molecule (SimParameters *simParams, Parameters *param, molfile_plugin_t *pIOHdl, void *pIOFileHdl, int natoms)
 Molecule (SimParameters *, Parameters *, Ambertoppar *)
void read_parm (Ambertoppar *)
 Molecule (SimParameters *, Parameters *, const GromacsTopFile *)
 ~Molecule ()
void send_Molecule (MOStream *)
void receive_Molecule (MIStream *)
void build_constraint_params (StringList *, StringList *, StringList *, PDB *, char *)
void build_gridforce_params (StringList *, StringList *, StringList *, StringList *, PDB *, char *)
void build_movdrag_params (StringList *, StringList *, StringList *, PDB *, char *)
void build_rotdrag_params (StringList *, StringList *, StringList *, StringList *, StringList *, StringList *, PDB *, char *)
void build_constorque_params (StringList *, StringList *, StringList *, StringList *, StringList *, StringList *, PDB *, char *)
void build_constant_forces (char *)
void build_langevin_params (BigReal coupling, BigReal drudeCoupling, Bool doHydrogen)
void build_langevin_params (StringList *, StringList *, PDB *, char *)
void build_fixed_atoms (StringList *, StringList *, PDB *, char *)
void build_stirred_atoms (StringList *, StringList *, PDB *, char *)
void build_extra_bonds (Parameters *parameters, StringList *file)
void build_fep_flags (StringList *, StringList *, PDB *, char *, const char *)
void delete_alch_bonded (void)
void build_ss_flags (const StringList *ssfile, const StringList *sscol, PDB *initial_pdb, const char *cwd)
void build_exPressure_atoms (StringList *, StringList *, PDB *, char *)
void print_go_sigmas ()
void build_go_sigmas (StringList *, char *)
void build_go_sigmas2 (StringList *, char *)
void build_go_arrays (StringList *, char *)
BigReal get_gro_force (BigReal, BigReal, BigReal, int, int) const
BigReal get_gro_force2 (BigReal, BigReal, BigReal, int, int, BigReal *, BigReal *) const
BigReal get_go_force (BigReal, int, int, BigReal *, BigReal *) const
BigReal get_go_force_new (BigReal, int, int, BigReal *, BigReal *) const
BigReal get_go_force2 (BigReal, BigReal, BigReal, int, int, BigReal *, BigReal *) const
Bool atoms_1to4 (unsigned int, unsigned int)
void reloadCharges (float charge[], int n)
Bool is_lp (int)
Bool is_drude (int)
Bool is_hydrogen (int)
Bool is_oxygen (int)
Bool is_hydrogenGroupParent (int)
Bool is_water (int)
int get_groupSize (int)
int get_mother_atom (int)
int get_cluster (int anum) const
int get_clusterSize (int anum) const
const float * getOccupancyData ()
void setOccupancyData (molfile_atom_t *atomarray)
void freeOccupancyData ()
const float * getBFactorData ()
void setBFactorData (molfile_atom_t *atomarray)
void freeBFactorData ()
Real atommass (int anum) const
Real atomcharge (int anum) const
Index atomvdwtype (int anum) const
Bondget_bond (int bnum) const
Angleget_angle (int anum) const
Improperget_improper (int inum) const
Dihedralget_dihedral (int dnum) const
Crosstermget_crossterm (int inum) const
Lphostget_lphost (int atomid) const
BondgetAllBonds () const
AnglegetAllAngles () const
ImpropergetAllImpropers () const
DihedralgetAllDihedrals () const
CrosstermgetAllCrossterms () const
LphostgetAllLphosts () const
Bondget_donor (int dnum) const
Bondget_acceptor (int dnum) const
BondgetAllDonors () const
BondgetAllAcceptors () const
Exclusionget_exclusion (int ex) const
const char * get_atomtype (int anum) const
int get_atom_from_name (const char *segid, int resid, const char *aname) const
int get_residue_size (const char *segid, int resid) const
int get_atom_from_index_in_residue (const char *segid, int resid, int index) const
int32get_bonds_for_atom (int anum)
int32get_angles_for_atom (int anum)
int32get_dihedrals_for_atom (int anum)
int32get_impropers_for_atom (int anum)
int32get_crossterms_for_atom (int anum)
int32get_exclusions_for_atom (int anum)
const int32get_full_exclusions_for_atom (int anum) const
const int32get_mod_exclusions_for_atom (int anum) const
int checkexcl (int atom1, int atom2) const
const ExclusionCheckget_excl_check_for_atom (int anum) const
Bool is_atom_gridforced (int atomnum, int gridnum) const
Bool is_atom_constrained (int atomnum) const
Bool is_atom_movdragged (int atomnum) const
Bool is_atom_rotdragged (int atomnum) const
Bool is_atom_constorqued (int atomnum) const
void get_cons_params (Real &k, Vector &refPos, int atomnum) const
void get_gridfrc_params (Real &k, Charge &q, int atomnum, int gridnum) const
GridforceGridget_gridfrc_grid (int gridnum) const
int set_gridfrc_grid (int gridnum, GridforceGrid *grid)
Real langevin_param (int atomnum) const
void get_stir_refPos (Vector &refPos, int atomnum) const
void put_stir_startTheta (Real theta, int atomnum) const
Real get_stir_startTheta (int atomnum) const
void get_movdrag_params (Vector &v, int atomnum) const
void get_rotdrag_params (BigReal &v, Vector &a, Vector &p, int atomnum) const
void get_constorque_params (BigReal &v, Vector &a, Vector &p, int atomnum) const
unsigned char get_fep_type (int anum) const
unsigned char get_ss_type (int anum) const
int get_fep_bonded_type (const int *atomID, unsigned int order) const
Bool is_atom_fixed (int atomnum) const
Bool is_atom_stirred (int atomnum) const
Bool is_group_fixed (int atomnum) const
Bool is_atom_exPressure (int atomnum) const
Real rigid_bond_length (int atomnum) const
void print_atoms (Parameters *)
void print_bonds (Parameters *)
void print_exclusions ()
void goInit ()
void build_gro_pair ()
void build_go_params (StringList *)
void read_go_file (char *)
Real get_go_cutoff (int chain1, int chain2)
Real get_go_epsilonRep (int chain1, int chain2)
Real get_go_sigmaRep (int chain1, int chain2)
Real get_go_epsilon (int chain1, int chain2)
int get_go_exp_a (int chain1, int chain2)
int get_go_exp_b (int chain1, int chain2)
int get_go_exp_rep (int chain1, int chain2)
Bool go_restricted (int, int, int)
void print_go_params ()
void initialize ()
void send_GoMolecule (MOStream *)
void receive_GoMolecule (MIStream *)

Public Attributes

int ss_num_vdw_params
int * ss_vdw_type
int * ss_index
int is_drude_psf
int is_lonepairs_psf
Real r_om
Real r_ohc
BigReal tail_corr_ener
BigReal tail_corr_virial
BigReal tail_corr_dUdl_1
BigReal tail_corr_dUdl_2
BigReal tail_corr_virial_1
BigReal tail_corr_virial_2
int numAtoms
int numRealBonds
int numBonds
int numAngles
int numDihedrals
int suspiciousAlchBonds
int alchDroppedAngles
int alchDroppedDihedrals
int alchDroppedImpropers
int numImpropers
int numCrossterms
int numDonors
int numAcceptors
int numExclusions
int64 numTotalExclusions
int numLonepairs
 Number of lone pairs.
int numDrudeAtoms
 Number of Drude particles.
int numTholes
 Number of Thole terms.
int numAnisos
 Number of anisotropic terms.
int numLphosts
 Number of lone pair host records in PSF.
int numZeroMassAtoms
 Number of atoms with zero mass.
int numConstraints
int numGridforceGrids
int * numGridforces
int numMovDrag
int numRotDrag
int numConsTorque
int numFixedAtoms
int numStirredAtoms
int numExPressureAtoms
int numHydrogenGroups
int maxHydrogenGroupSize
int numMigrationGroups
int maxMigrationGroupSize
int numFixedGroups
int numRigidBonds
int numFixedRigidBonds
int numFepInitial
int numFepFinal
int numConsForce
int32consForceIndexes
VectorconsForce
int32consTorqueIndexes
ConsTorqueParams * consTorqueParams
int numCalcBonds
int numCalcAngles
int numCalcDihedrals
int numCalcImpropers
int numCalcCrossterms
int64 numCalcExclusions
int64 numCalcFullExclusions
int numCalcTholes
int numCalcAnisos
int numMultipleDihedrals
int numMultipleImpropers
HydrogenGroup hydrogenGroup
int numGoAtoms
int32atomChainTypes
int32goSigmaIndices
int32goResidIndices
RealgoSigmas
bool * goWithinCutoff
RealgoCoordinates
int * goResids
PDBgoPDB
int goNumLJPair
int * goIndxLJA
int * goIndxLJB
double * goSigmaPairA
double * goSigmaPairB
int * pointerToGoBeg
int * pointerToGoEnd
int numPair
int numLJPair
int numCalcLJPair
int * pointerToLJBeg
int * pointerToLJEnd
int * indxLJA
int * indxLJB
RealpairC6
RealpairC12
int * gromacsPair_type
int * pointerToGaussBeg
int * pointerToGaussEnd
int numGaussPair
int * indxGaussA
int * indxGaussB
RealgA
RealgMu1
RealgiSigma1
RealgMu2
RealgiSigma2
RealgRepulsive
BigReal energyNative
BigReal energyNonnative
int isOccupancyValid
int isBFactorValid
GoValue go_array [MAX_GO_CHAINS *MAX_GO_CHAINS]
int go_indices [MAX_GO_CHAINS+1]
int NumGoChains

Friends

class ExclElem
class BondElem
class AngleElem
class DihedralElem
class ImproperElem
class TholeElem
class AnisoElem
class CrosstermElem
class GromacsPairElem
class WorkDistrib

Detailed Description

Definition at line 150 of file Molecule.h.


Constructor & Destructor Documentation

Molecule::Molecule ( SimParameters simParams,
Parameters param 
)

Definition at line 423 of file Molecule.C.

References initialize().

00424 {
00425   initialize(simParams,param);
00426 }

Molecule::Molecule ( SimParameters simParams,
Parameters param,
char *  filename,
ConfigList cfgList = NULL 
)

Definition at line 436 of file Molecule.C.

References initialize(), SimParameters::LCPOOn, and SimParameters::useCompressedPsf.

00437 {
00438   initialize(simParams,param);
00439 
00440 #ifdef MEM_OPT_VERSION
00441   if(simParams->useCompressedPsf)
00442       read_mol_signatures(filename, param, cfgList);
00443 #else
00444         read_psf_file(filename, param); 
00445  //LCPO
00446   if (simParams->LCPOOn)
00447     assignLCPOTypes( 0 );
00448 #endif      
00449  }

Molecule::Molecule ( SimParameters simParams,
Parameters param,
molfile_plugin_t *  pIOHdl,
void *  pIOFileHdl,
int  natoms 
)

Definition at line 458 of file Molecule.C.

References initialize(), SimParameters::LCPOOn, NAMD_die(), numAngles, numAtoms, numBonds, numCrossterms, numDihedrals, numImpropers, numRealBonds, setBFactorData(), and setOccupancyData().

00459 {
00460 #ifdef MEM_OPT_VERSION
00461   NAMD_die("Sorry, plugin IO is not supported in the memory optimized version.");
00462 #else
00463     initialize(simParams, param);
00464     numAtoms = natoms;
00465     int optflags = MOLFILE_BADOPTIONS;
00466     molfile_atom_t *atomarray = (molfile_atom_t *) malloc(natoms*sizeof(molfile_atom_t));
00467     memset(atomarray, 0, natoms*sizeof(molfile_atom_t));
00468 
00469     //1a. read basic atoms information
00470     int rc = pIOHdl->read_structure(pIOFileHdl, &optflags, atomarray);
00471     if (rc != MOLFILE_SUCCESS && rc != MOLFILE_NOSTRUCTUREDATA) {
00472         free(atomarray);
00473         NAMD_die("ERROR: plugin failed reading structure data");
00474     }
00475     if(optflags == MOLFILE_BADOPTIONS) {
00476         free(atomarray);
00477         NAMD_die("ERROR: plugin didn't initialize optional data flags");
00478     }
00479     if(optflags & MOLFILE_OCCUPANCY) {
00480         setOccupancyData(atomarray);
00481     }
00482     if(optflags & MOLFILE_BFACTOR) {
00483         setBFactorData(atomarray);
00484     }
00485     //1b. load basic atoms information to the molecule object
00486     plgLoadAtomBasics(atomarray);    
00487     free(atomarray);
00488 
00489     //2a. read bonds
00490     //indices are one-based in read_bonds
00491     int *from, *to;
00492     float *bondorder;
00493     int *bondtype, nbondtypes;
00494     char **bondtypename;
00495     if(pIOHdl->read_bonds!=NULL) {
00496         if(pIOHdl->read_bonds(pIOFileHdl, &numBonds, &from, &to, &bondorder,
00497                                  &bondtype, &nbondtypes, &bondtypename)){
00498             NAMD_die("ERROR: failed reading bond information.");
00499         }
00500     }    
00501     //2b. load bonds information to the molecule object
00502     if(numBonds!=0) {
00503         plgLoadBonds(from,to);
00504     }
00505 
00506     //3a. read other bonded structures
00507     int *plgAngles, *plgDihedrals, *plgImpropers, *plgCterms;
00508     int ctermcols, ctermrows;
00509     int *angletypes, numangletypes, *dihedraltypes, numdihedraltypes;
00510     int *impropertypes, numimpropertypes; 
00511     char **angletypenames, **dihedraltypenames, **impropertypenames;
00512 
00513     plgAngles=plgDihedrals=plgImpropers=plgCterms=NULL;
00514     if(pIOHdl->read_angles!=NULL) {
00515         if(pIOHdl->read_angles(pIOFileHdl,
00516                   &numAngles, &plgAngles,
00517                   &angletypes, &numangletypes, &angletypenames,
00518                   &numDihedrals, &plgDihedrals,
00519                   &dihedraltypes, &numdihedraltypes, &dihedraltypenames,
00520                   &numImpropers, &plgImpropers,
00521                   &impropertypes, &numimpropertypes, &impropertypenames,
00522                   &numCrossterms, &plgCterms, &ctermcols, &ctermrows)) {
00523             NAMD_die("ERROR: failed reading angle information.");
00524         }
00525     }
00526     //3b. load other bonded structures to the molecule object
00527     if(numAngles!=0) plgLoadAngles(plgAngles);
00528     if(numDihedrals!=0) plgLoadDihedrals(plgDihedrals);
00529     if(numImpropers!=0) plgLoadImpropers(plgImpropers);
00530     if(numCrossterms!=0) plgLoadCrossterms(plgCterms);
00531 
00532   numRealBonds = numBonds;
00533   build_atom_status();
00534   //LCPO
00535   if (simParams->LCPOOn)
00536     assignLCPOTypes( 2 );
00537 #endif
00538 }

Molecule::Molecule ( SimParameters ,
Parameters ,
Ambertoppar  
)
Molecule::Molecule ( SimParameters ,
Parameters ,
const GromacsTopFile  
)
Molecule::~Molecule (  ) 

Definition at line 552 of file Molecule.C.

References atomSigPool, numAtoms, ss_index, and ss_vdw_type.

00553 {
00554   /*  Check to see if each array was ever allocated.  If it was   */
00555   /*  then free it            */
00556   if (atoms != NULL)
00557     delete [] atoms;
00558 
00559   if (atomNames != NULL)
00560   {
00561     // subarrarys allocated from arena - automatically deleted
00562     delete [] atomNames;
00563   }
00564   delete nameArena;
00565 
00566   if (resLookup != NULL)
00567     delete resLookup;
00568 
00569   // DRUDE: free arrays read from PSF
00570   if (drudeConsts != NULL) delete [] drudeConsts;
00571   if (lphosts != NULL) delete [] lphosts;
00572   if (anisos != NULL) delete [] anisos;
00573   if (tholes != NULL) delete [] tholes;
00574   if (lphostIndexes != NULL) delete [] lphostIndexes;
00575   // DRUDE
00576 
00577   //LCPO
00578   if (lcpoParamType != NULL) delete [] lcpoParamType;
00579 
00580   #ifdef MEM_OPT_VERSION
00581   if(eachAtomSig) delete [] eachAtomSig;
00582   if(atomSigPool) delete [] atomSigPool;
00583   #else
00584   if (bonds != NULL)
00585     delete [] bonds;
00586 
00587   if (angles != NULL)
00588     delete [] angles;
00589 
00590   if (dihedrals != NULL)
00591     delete [] dihedrals;
00592 
00593   if (impropers != NULL)
00594     delete [] impropers;
00595 
00596   if (crossterms != NULL)
00597     delete [] crossterms;
00598 
00599   if (exclusions != NULL)
00600     delete [] exclusions;
00601   #endif
00602 
00603   if (donors != NULL)
00604     delete [] donors;
00605 
00606   if (acceptors != NULL)
00607     delete [] acceptors;  
00608 
00609   #ifdef MEM_OPT_VERSION
00610   if(exclSigPool) delete [] exclSigPool;
00611   if(exclChkSigPool) delete [] exclChkSigPool;
00612   if(eachAtomExclSig) delete [] eachAtomExclSig;
00613   if(fixedAtomsSet) delete fixedAtomsSet;
00614   if(constrainedAtomsSet) delete constrainedAtomsSet;
00615   #else
00616   if (bondsByAtom != NULL)
00617        delete [] bondsByAtom;
00618   
00619   if (anglesByAtom != NULL)
00620        delete [] anglesByAtom;
00621   
00622   if (dihedralsByAtom != NULL)
00623        delete [] dihedralsByAtom;
00624   
00625   if (impropersByAtom != NULL)
00626        delete [] impropersByAtom;
00627   
00628   if (crosstermsByAtom != NULL)
00629        delete [] crosstermsByAtom;  
00630 
00631   if (exclusionsByAtom != NULL)
00632        delete [] exclusionsByAtom;
00633   
00634   if (fullExclusionsByAtom != NULL)
00635        delete [] fullExclusionsByAtom;
00636   
00637   if (modExclusionsByAtom != NULL)
00638        delete [] modExclusionsByAtom;
00639   
00640   if (all_exclusions != NULL)
00641        delete [] all_exclusions;
00642 
00643   // JLai
00644   if (gromacsPairByAtom != NULL)
00645       delete [] gromacsPairByAtom;
00646   // End of JLai
00647 
00648   // DRUDE
00649   if (tholesByAtom != NULL)
00650        delete [] tholesByAtom;
00651   if (anisosByAtom != NULL)
00652        delete [] anisosByAtom;
00653   // DRUDE
00654   #endif
00655 
00656   //LCPO
00657   if (lcpoParamType != NULL)
00658     delete [] lcpoParamType;
00659 
00660   if (fixedAtomFlags != NULL)
00661        delete [] fixedAtomFlags;
00662 
00663   if (stirIndexes != NULL)
00664     delete [] stirIndexes;
00665 
00666 
00667   #ifdef MEM_OPT_VERSION
00668   if(clusterSigs != NULL){      
00669       delete [] clusterSigs;
00670   }  
00671   #else
00672   if (cluster != NULL)
00673        delete [] cluster;  
00674   #endif
00675   if (clusterSize != NULL)
00676        delete [] clusterSize;
00677 
00678   if (exPressureAtomFlags != NULL)
00679        delete [] exPressureAtomFlags;
00680 
00681   if (rigidBondLengths != NULL)
00682        delete [] rigidBondLengths;
00683 
00684 //fepb
00685   if (fepAtomFlags != NULL)
00686        delete [] fepAtomFlags;
00687 //fepe
00688 
00689 //soluteScaling
00690   if (ssAtomFlags != NULL)
00691        delete [] ssAtomFlags;
00692   if (ss_vdw_type != NULL)
00693        delete [] ss_vdw_type;
00694   if (ss_index != NULL)
00695        delete [] ss_index;
00696 //soluteScaling
00697 
00698   if (qmAtomGroup != NULL)
00699        delete [] qmAtomGroup;
00700   
00701   if (qmAtmIndx != NULL)
00702        delete [] qmAtmIndx;
00703   
00704   if (qmAtmChrg != NULL)
00705        delete [] qmAtmChrg;
00706   
00707   
00708   if (qmGrpNumBonds != NULL)
00709        delete [] qmGrpNumBonds;
00710   
00711   if (qmGrpSizes != NULL)
00712        delete [] qmGrpSizes;
00713   
00714   if (qmDummyBondVal != NULL)
00715        delete [] qmDummyBondVal;
00716   
00717   if (qmMMNumTargs != NULL)
00718        delete [] qmMMNumTargs;
00719   
00720   if (qmGrpID != NULL)
00721        delete [] qmGrpID;
00722   
00723   if (qmGrpChrg != NULL)
00724        delete [] qmGrpChrg;
00725   
00726   if (qmGrpMult != NULL)
00727        delete [] qmGrpMult;
00728   
00729   if (qmMeMMindx != NULL)
00730        delete [] qmMeMMindx;
00731   
00732   if (qmMeQMGrp != NULL)
00733        delete [] qmMeQMGrp;
00734   
00735   if (qmLSSSize != NULL)
00736        delete [] qmLSSSize;
00737   
00738   if (qmLSSIdxs != NULL)
00739        delete [] qmLSSIdxs;
00740   
00741   if (qmLSSMass != NULL)
00742        delete [] qmLSSMass;
00743   
00744   if (qmLSSRefSize != NULL)
00745        delete [] qmLSSRefSize;
00746   
00747   if (qmLSSRefIDs != NULL)
00748        delete [] qmLSSRefIDs;
00749   
00750   if (qmLSSRefMass != NULL)
00751        delete [] qmLSSRefMass;
00752   
00753   if (qmMMBond != NULL) {
00754       for(int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
00755           if (qmMMBond[grpIndx] != NULL)
00756               delete [] qmMMBond[grpIndx];
00757       }
00758       delete [] qmMMBond;
00759   }
00760   
00761   if (qmGrpBonds != NULL) {
00762       for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
00763           if (qmGrpBonds[grpIndx] != NULL)
00764               delete [] qmGrpBonds[grpIndx];
00765       }
00766       delete [] qmGrpBonds;
00767   }
00768   
00769   if (qmMMBondedIndx != NULL) {
00770       for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
00771           if (qmMMBondedIndx[grpIndx] != NULL)
00772               delete [] qmMMBondedIndx[grpIndx];
00773       }
00774       delete [] qmMMBondedIndx;
00775   }
00776   
00777   if (qmMMChargeTarget != NULL) {
00778       for(int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
00779           if (qmMMChargeTarget[grpIndx] != NULL)
00780               delete [] qmMMChargeTarget[grpIndx];
00781       }
00782       delete [] qmMMChargeTarget;
00783   }
00784   
00785     if (qmElementArray != NULL){
00786         for(int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
00787             if (qmElementArray[atmInd] != NULL)
00788                 delete [] qmElementArray[atmInd];
00789         }
00790         delete [] qmElementArray;
00791     }
00792     
00793     if (qmDummyElement != NULL){
00794         for(int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
00795             if (qmDummyElement[atmInd] != NULL)
00796                 delete [] qmDummyElement[atmInd];
00797         }
00798         delete [] qmDummyElement;
00799     }
00800 
00801     if (qmCustPCSizes != NULL){
00802         delete [] qmCustPCSizes;
00803     }
00804     
00805     if (qmCustomPCIdxs != NULL){
00806         delete [] qmCustomPCIdxs;
00807     }
00808     
00809     if (cSMDindex != NULL) {
00810       for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
00811           if (cSMDindex[grpIndx] != NULL)
00812               delete [] cSMDindex[grpIndx];
00813       }
00814       delete [] cSMDindex;
00815     }
00816     
00817     if (cSMDpairs != NULL) {
00818       for(int instIndx = 0 ; instIndx < cSMDnumInst; instIndx++) {
00819           if (cSMDpairs[instIndx] != NULL)
00820               delete [] cSMDpairs[instIndx];
00821       }
00822       delete [] cSMDpairs;
00823     }
00824     
00825     if (cSMDindxLen != NULL)
00826         delete [] cSMDindxLen;
00827     if (cSMDKs != NULL)
00828         delete [] cSMDKs;
00829     if (cSMDVels != NULL)
00830         delete [] cSMDVels;
00831     if (cSMDcoffs != NULL)
00832         delete [] cSMDcoffs;
00833     
00834   #ifndef MEM_OPT_VERSION
00835   delete arena;
00836   delete exclArena;
00837   #endif
00838 }


Member Function Documentation

Real Molecule::atomcharge ( int  anum  )  const [inline]

Definition at line 1035 of file Molecule.h.

References atoms.

Referenced by WorkDistrib::createAtomLists(), NamdState::loadStructure(), Sequencer::reloadCharges(), Sequencer::rescaleSoluteCharges(), colvarproxy_namd::update_atom_properties(), and colvarproxy_namd::update_group_properties().

01036   {
01037     #ifdef MEM_OPT_VERSION
01038     return atomChargePool[eachAtomCharge[anum]];
01039     #else
01040     return(atoms[anum].charge);
01041     #endif
01042   }

Real Molecule::atommass ( int  anum  )  const [inline]

Definition at line 1025 of file Molecule.h.

References atoms.

Referenced by WorkDistrib::createAtomLists(), GlobalMasterFreeEnergy::getMass(), GlobalMasterEasy::getMass(), NamdState::loadStructure(), Tcl_centerOfMass(), Tcl_radiusOfGyration(), colvarproxy_namd::update_atom_properties(), and colvarproxy_namd::update_group_properties().

01026   {
01027     #ifdef MEM_OPT_VERSION
01028     return atomMassPool[eachAtomMass[anum]];
01029     #else
01030     return(atoms[anum].mass);
01031     #endif
01032   }

Bool Molecule::atoms_1to4 ( unsigned int  atom1,
unsigned int  atom2 
)

Definition at line 1537 of file GoMolecule.C.

References dihedral::atom1, angle::atom1, bond::atom1, dihedral::atom2, angle::atom2, bond::atom2, dihedral::atom3, angle::atom3, dihedral::atom4, DebugM, FALSE, get_angle(), get_angles_for_atom(), get_bond(), get_bonds_for_atom(), get_dihedral(), get_dihedrals_for_atom(), and TRUE.

Referenced by build_go_arrays(), build_go_sigmas(), and build_go_sigmas2().

01539 {
01540   int bondNum;   //  Bonds in bonded list
01541   int angleNum;  //  Angles in angle list
01542   int dihedralNum;   //  Dihedrals in dihedral list
01543   int *bonds;
01544   int *angles;
01545   int *dihedrals;
01546   Bond *bond;     //  Temporary bond structure
01547   Angle *angle;   //  Temporary angle structure
01548   Dihedral *dihedral; //  Temporary dihedral structure
01549 
01550   DebugM(2,"atoms_1to4(" << atom1 << "," << atom2 << ")" << std::endl);
01551 
01552   bonds = get_bonds_for_atom(atom1);
01553   bondNum = *bonds;
01554   while(bondNum != -1) {
01555     bond = get_bond(bondNum);
01556     DebugM(2,"bond  atom1:" << bond->atom1 << ", atom2:" << bond->atom2 << std::endl);
01557     if (atom2 == bond->atom1 || atom2 == bond->atom2) {
01558       return TRUE;
01559     }
01560     bondNum = *(++bonds);
01561   }
01562 
01563   bonds = get_bonds_for_atom(atom2);
01564   bondNum = *bonds;
01565   while(bondNum != -1) {
01566     bond = get_bond(bondNum);
01567     DebugM(2,"bond  atom1:" << bond->atom1 << ", atom2:" << bond->atom2 << std::endl);
01568     if (atom1 == bond->atom1 || atom1 == bond->atom2) {
01569       return TRUE;
01570     }
01571     bondNum = *(++bonds);
01572   }
01573 
01574   angles = get_angles_for_atom(atom1);
01575   angleNum = *angles;
01576   while(angleNum != -1) {
01577     angle = get_angle(angleNum);
01578     DebugM(2,"angle  atom1:" << angle->atom1 << ", atom2:" << angle->atom2 << ", atom3:" << angle->atom3 << std::endl);
01579     if (atom2 == angle->atom1 || atom2 == angle->atom2 || atom2 == angle->atom3) {
01580       return TRUE;
01581     }
01582     angleNum = *(++angles);
01583   }
01584 
01585   angles = get_angles_for_atom(atom2);
01586   angleNum = *angles;
01587   while(angleNum != -1) {
01588     angle = get_angle(angleNum);
01589     DebugM(2,"angle  atom1:" << angle->atom1 << ", atom2:" << angle->atom2 << ", atom3:" << angle->atom3 << std::endl);
01590     if (atom1 == angle->atom1 || atom1 == angle->atom2 || atom1 == angle->atom3) {
01591       return TRUE;
01592     }
01593     angleNum = *(++angles);
01594   }
01595 
01596   dihedrals = get_dihedrals_for_atom(atom1);
01597   dihedralNum = *dihedrals;
01598   while(dihedralNum != -1) {
01599     dihedral = get_dihedral(dihedralNum);
01600     DebugM(2,"dihedral  atom1:" << dihedral->atom1 << ", atom2:" << dihedral->atom2 << ", atom3:" << dihedral->atom3 << ", atom4:" << dihedral->atom4 << std::endl);
01601     if (atom2 == dihedral->atom1 || atom2 == dihedral->atom2 \
01602         || atom2 == dihedral->atom3 || atom2 == dihedral->atom4) {
01603       return TRUE;
01604     }
01605     dihedralNum = *(++dihedrals);
01606   }
01607 
01608   dihedrals = get_dihedrals_for_atom(atom2);
01609   dihedralNum = *dihedrals;
01610   while(dihedralNum != -1) {
01611     dihedral = get_dihedral(dihedralNum);
01612     DebugM(2,"dihedral  atom1:" << dihedral->atom1 << ", atom2:" << dihedral->atom2 << ", atom3:" << dihedral->atom3 << ", atom4:" << dihedral->atom4 << std::endl);
01613     if (atom1 == dihedral->atom1 || atom1 == dihedral->atom2 \
01614         || atom1 == dihedral->atom3 || atom1 == dihedral->atom4) {
01615       return TRUE;
01616     }
01617     dihedralNum = *(++dihedrals);
01618   }
01619   
01620   return FALSE;
01621 }

Index Molecule::atomvdwtype ( int  anum  )  const [inline]

Definition at line 1045 of file Molecule.h.

References atoms.

Referenced by WorkDistrib::createAtomLists(), and dumpbench().

01046   {      
01047       return(atoms[anum].vdw_type);
01048   }

void Molecule::build_constant_forces ( char *   ) 
void Molecule::build_constorque_params ( StringList ,
StringList ,
StringList ,
StringList ,
StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_constraint_params ( StringList ,
StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_exPressure_atoms ( StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_extra_bonds ( Parameters parameters,
StringList file 
)
void Molecule::build_fep_flags ( StringList ,
StringList ,
PDB ,
char *  ,
const char *   
)
void Molecule::build_fixed_atoms ( StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_go_arrays ( StringList goCoordFile,
char *  cwd 
)

goSigmas = new Real[numGoAtoms*numGoAtoms]; goWithinCutoff = new bool[numGoAtoms*numGoAtoms]; for (i=0; i<numGoAtoms; i++) { for (j=0; j<numGoAtoms; j++) { goSigmas[i*numGoAtoms + j] = 0.0; goWithinCutoff[i*numGoAtoms + j] = false; } }

Definition at line 950 of file GoMolecule.C.

References PDB::atom(), atomChainTypes, atoms_1to4(), StringList::data, DebugM, endi(), energyNative, energyNonnative, get_go_cutoff(), go_restricted(), goCoordinates, goPDB, goResids, goSigmaIndices, iINFO(), iout, j, NAMD_die(), StringList::next, PDB::num_atoms(), numAtoms, numGoAtoms, PDBAtom::residueseq(), PDBAtom::xcoor(), PDBAtom::ycoor(), and PDBAtom::zcoor().

Referenced by NamdState::loadStructure().

00952 {
00953   DebugM(3,"->build_go_arrays" << std::endl);
00954   //PDB *goPDB;      //  Pointer to PDB object to use
00955   int bcol = 4;      //  Column that data is in
00956   int32 chainType = 0;      //  b value from PDB file
00957   int i;         //  Loop counter
00958   int j;         //  Loop counter
00959   BigReal atomAtomDist;     //  Distance between two atoms -- JLai put back 
00960   int resid1;    //  Residue ID for first atom
00961   int resid2;    //  Residue ID for second atom
00962   int residDiff;     //  Difference between resid1 and resid2
00963   int goIndex;       //  Index into the goCoordinates array
00964   int goIndx;        //  Index into the goCoordinates array
00965   char filename[129];    //  Filename
00966   
00967   //JLai
00968   BigReal nativeEnergy, *native;
00969   BigReal nonnativeEnergy, *nonnative;
00970   nativeEnergy = 0;
00971   nonnativeEnergy = 0;
00972   native = &nativeEnergy;
00973   nonnative = &nonnativeEnergy;
00974 
00975   long nativeContacts = 0;
00976   long nonnativeContacts = 0;
00977 
00978   //  Get the PDB object that contains the Go coordinates.  If
00979   //  the user gave another file name, use it.  Otherwise, just use
00980   //  the PDB file that has the initial coordinates.  
00981   if (goCoordFile == NULL)
00982     {
00983       //goPDB = initial_pdb;
00984       NAMD_die("Error: goCoordFile is NULL - build_go_arrays");
00985     }
00986   else
00987   {
00988     if (goCoordFile->next != NULL)
00989       {
00990         NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
00991       }
00992     
00993     if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
00994       {
00995         strcpy(filename, goCoordFile->data);
00996       }
00997     else
00998       {
00999         strcpy(filename, cwd);
01000         strcat(filename, goCoordFile->data);
01001       }
01002     
01003     goPDB = new PDB(filename);
01004     if ( goPDB == NULL )
01005       {
01006         NAMD_die("goPDB memory allocation failed in Molecule::build_go_arrays");
01007       }
01008     
01009     if (goPDB->num_atoms() != numAtoms)
01010       {
01011         NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
01012       }
01013   }
01014   
01015   //  Allocate the array to hold Go atom indices into the sigma array
01016   goSigmaIndices = new int32[numAtoms];
01017   
01018   if (goSigmaIndices == NULL) {
01019     NAMD_die("goSigmaIndices memory allocation failed in Molecule::build_go_arrays");
01020   }
01021   
01022   numGoAtoms = 0;
01023   
01024   //  Loop through all the atoms and get the Go chain types
01025   for (i=0; i<numAtoms; i++) {
01026     chainType = (int32)(goPDB->atom(i))->occupancy();
01027     if ( chainType != 0 ) {
01028       DebugM(3,"build_go_arrays - atom:" << i << std::endl);
01029       goSigmaIndices[i] = numGoAtoms;
01030       numGoAtoms++;
01031     }
01032     else {
01033       goSigmaIndices[i] = -1;
01034     }
01035   }
01036 
01037   // Allocate the array to hold the sigma values for each Go atom pair
01049   //  Allocate the array to hold the chain types
01050   atomChainTypes = new int32[numGoAtoms];
01051 
01052   if (atomChainTypes == NULL) {
01053     NAMD_die("atomChainTypes memory allocation failed in Molecule::build_go_arrays");
01054   }
01055 
01056   // Allocate the array to hold (x,y,z) coordinates for all of the Go atoms
01057   goCoordinates = new Real[numGoAtoms*3];
01058 
01059   if (goCoordinates == NULL) {
01060     NAMD_die("goCoordinates memory allocation failed in Molecule::build_go_arrays");
01061   }
01062 
01063   goResids = new int[numGoAtoms];
01064 
01065   // Allocate the array to hold PDB residu IDs for all of the Go atoms
01066   if (goResids == NULL) {
01067     NAMD_die("goResids memory allocation failed in Molecule::build_go_arrays");
01068   }
01069   
01070   for (i=0; i<numAtoms; i++) {
01071     goIndex = goSigmaIndices[i];
01072     if (goIndex != -1) {
01073       //  Assign the chainType value!
01074       //  Get the chainType from the occupancy field
01075       atomChainTypes[goIndex] = (int32)(goPDB->atom(i))->occupancy();
01076       goCoordinates[goIndex*3] = goPDB->atom(i)->xcoor();
01077       goCoordinates[goIndex*3 + 1] = goPDB->atom(i)->ycoor();
01078       goCoordinates[goIndex*3 + 2] = goPDB->atom(i)->zcoor();
01079       goResids[goIndex] = goPDB->atom(i)->residueseq();
01080     }
01081   }
01082       // JLai
01083   energyNative = 0;
01084   energyNonnative = 0;
01085   //printf("INIT ENERGY: (N) %f (NN) %f\n", energyNative, energyNonnative);
01086   for (i=0; i<numAtoms-1; i++) {
01087     goIndex = goSigmaIndices[i];
01088     if (goIndex != -1) {
01089       for (j=i+1; j<numAtoms;j++) {
01090         goIndx = goSigmaIndices[j];
01091         if (goIndx != -1) {
01092           resid1 = (goPDB->atom(i))->residueseq();
01093           resid2 = (goPDB->atom(j))->residueseq();
01094           residDiff = resid2 - resid1;
01095           if (residDiff < 0) residDiff = -residDiff;
01096           if (atomChainTypes[goIndex] && atomChainTypes[goIndx] &&
01097               !(this->go_restricted(atomChainTypes[goIndex],atomChainTypes[goIndx],residDiff)) &&
01098               !atoms_1to4(i,j)) {
01099             atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
01100                                 pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
01101                                 pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
01102             if (atomAtomDist <= this->get_go_cutoff(atomChainTypes[goIndex],atomChainTypes[goIndx]) ) {
01103               nativeContacts++;
01104             } else {
01105               nonnativeContacts++;
01106             }    
01107           }
01108         }
01109       }
01110     }
01111   }
01112   iout << iINFO << "Number of UNIQUE    native contacts: " << nativeContacts     << "\n" << endi;
01113   iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts  << "\n" << endi;
01114   
01115   //  If we had to create a PDB object, delete it now
01116   if (goCoordFile != NULL) {
01117     delete goPDB;
01118   }
01119   
01120   return;
01121 }

void Molecule::build_go_params ( StringList g  ) 

Definition at line 79 of file GoMolecule.C.

References StringList::data, endi(), iINFO(), iout, NAMD_die(), StringList::next, and read_go_file().

Referenced by NamdState::loadStructure().

00079                                             {
00080   iout << iINFO << "Building Go Parameters" << "\n" << endi;
00081 #ifdef MEM_OPT_VERSION
00082   NAMD_die("Go forces are not supported in memory-optimized builds.");
00083 #else
00084   build_lists_by_atom();
00085 #endif
00086   int iterator = 0;
00087     do
00088     {
00089       iout << iINFO << "Reading Go File: " << iterator << "\n" << endi;
00090       read_go_file(g->data);
00091       g = g->next;
00092       iterator++;
00093     } while ( g != NULL && iterator < 100);    
00094 }

void Molecule::build_go_sigmas ( StringList goCoordFile,
char *  cwd 
)

Definition at line 577 of file GoMolecule.C.

References PDB::atom(), atomChainTypes, atoms_1to4(), StringList::data, DebugM, endi(), get_go_cutoff(), get_go_exp_a(), get_go_exp_b(), go_restricted(), goPDB, goSigmaIndices, goSigmas, goWithinCutoff, iINFO(), iout, j, NAMD_die(), StringList::next, PDB::num_atoms(), numAtoms, and numGoAtoms.

Referenced by NamdState::loadStructure().

00579 {
00580   DebugM(3,"->build_go_sigmas" << std::endl);
00581   PDB *goPDB;      //  Pointer to PDB object to use
00582   int bcol = 4;      //  Column that data is in
00583   int32 chainType = 0;      //  b value from PDB file
00584   int i;         //  Loop counter
00585   int j;         //  Loop counter
00586   int resid1;    //  Residue ID for first atom
00587   int resid2;    //  Residue ID for second atom
00588   int residDiff;     //  Difference between resid1 and resid2
00589   Real sigma;    //  Sigma calculated for a Go atom pair
00590   Real atomAtomDist;     //  Distance between two atoms
00591   Real exp_a;            //  First exponent in L-J like formula
00592   Real exp_b;            //  Second exponent in L-J like formula
00593   char filename[129];    //  Filename
00594   
00595   //JLai
00596   BigReal nativeEnergy, *native;
00597   BigReal nonnativeEnergy, *nonnative;
00598   nativeEnergy = 0;
00599   nonnativeEnergy = 0;
00600   native = &nativeEnergy;
00601   nonnative = &nonnativeEnergy;
00602 
00603   long nativeContacts = 0;
00604   long nonnativeContacts = 0;
00605 
00606   //  Get the PDB object that contains the Go coordinates.  If
00607   //  the user gave another file name, use it.  Otherwise, just use
00608   //  the PDB file that has the initial coordinates.  
00609   if (goCoordFile == NULL)
00610     {
00611       //goPDB = initial_pdb;
00612       NAMD_die("Error: goCoordFile is NULL - build_go_sigmas");
00613     }
00614   else
00615   {
00616     if (goCoordFile->next != NULL)
00617       {
00618         NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
00619       }
00620     
00621     if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
00622       {
00623         strcpy(filename, goCoordFile->data);
00624       }
00625     else
00626       {
00627         strcpy(filename, cwd);
00628         strcat(filename, goCoordFile->data);
00629       }
00630     
00631     goPDB = new PDB(filename);
00632     if ( goPDB == NULL )
00633       {
00634         NAMD_die("Memory allocation failed in Molecule::build_go_sigmas");
00635       }
00636     
00637     if (goPDB->num_atoms() != numAtoms)
00638       {
00639         NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
00640       }
00641   }
00642   //  Allocate the array to hold the chain types
00643   atomChainTypes = new int32[numAtoms];
00644   //  Allocate the array to hold Go atom indices into the sigma array
00645   goSigmaIndices = new int32[numAtoms];
00646   
00647   if (atomChainTypes == NULL) {
00648     NAMD_die("memory allocation failed in Molecule::build_go_sigmas");
00649   }
00650   
00651   numGoAtoms = 0;
00652   
00653   //  Loop through all the atoms and get the Go chain types
00654   for (i=0; i<numAtoms; i++) {
00655     //  Get the chainType from the occupancy field
00656     chainType = (int32)(goPDB->atom(i))->occupancy();
00657     //  Assign the chainType value
00658     if ( chainType != 0 ) {
00659       //DebugM(3,"build_go_sigmas - atom:" << i << ", chainType:" << chainType << std::endl);
00660       atomChainTypes[i] = chainType;
00661       goSigmaIndices[i] = numGoAtoms;
00662       numGoAtoms++;
00663     }
00664     else {
00665       atomChainTypes[i] = 0;
00666       goSigmaIndices[i] = -1;
00667     }
00668     //printf("CT: %d %d %d %d\n",i,numGoAtoms,atomChainTypes[i],goSigmaIndices[i]);
00669   }
00670 
00671   // Allocate the array to hold the sigma values for each Go atom pair
00672   goSigmas = new Real[numGoAtoms*numGoAtoms];
00673   goWithinCutoff = new bool[numGoAtoms*numGoAtoms];
00674   for (i=0; i<numGoAtoms; i++) {
00675     for (j=0; j<numGoAtoms; j++) {
00676       goSigmas[i*numGoAtoms + j] = -1.0;
00677       goWithinCutoff[i*numGoAtoms + j] = false;
00678     }
00679   }
00680   //  Loop through all atom pairs and calculate sigmas
00681   DebugM(3,"    numAtoms=" << numAtoms << std::endl);
00682   for (i=0; i<numAtoms; i++) {
00683     //DebugM(3,"    i=" << i << std::endl);
00684     resid1 = (goPDB->atom(i))->residueseq();
00685     //DebugM(3,"    resid1=" << resid1 << std::endl);
00686     //if ( goSigmaIndices[i] != -1) {
00687     //  goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[i]] = 0.0;
00688     //}
00689      for (j=i+1; j<numAtoms; j++) {
00690       //DebugM(3,"    j=" << j << std::endl);
00691       resid2 = (goPDB->atom(j))->residueseq();
00692       //printf("GSIi %d %d %d\n",i,numAtoms,goSigmaIndices[i]);
00693       //printf("SAN CHECK: %d\n",goSigmaIndices[37]);
00694       //printf("GSIj %d %d %d\n",j,numAtoms,goSigmaIndices[j]);
00695       //printf("ATOMS_1to4 %d\n",!atoms_1to4(i,j));
00696       //DebugM(3,"    resid2=" << resid2 << std::endl);
00697       //  if goSigmaIndices aren't defined, don't set anything in goSigmas
00698       if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 && !atoms_1to4(i,j) ) {
00699         //printf("TAKING DIFFERENCE\n");
00700         residDiff = resid2 - resid1;
00701         //printf("RESIDDIFF %d\n",residDiff);
00702         if (residDiff < 0) residDiff = -residDiff;
00703         //printf("RESIDDIFF2 %d\n",residDiff);
00704         //  if this is a Go atom pair that is not restricted
00705         //    calculate sigma
00706         //  sigmas are initially set to -1.0 if the atom pair fails go_restricted
00707         //printf("CHECKING LOOPING\n");
00708         if ( atomChainTypes[i] && atomChainTypes[j] &&
00709              !(this->go_restricted(atomChainTypes[i],atomChainTypes[j],residDiff)) ) {
00710           atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
00711                               pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
00712                               pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
00713           if ( atomAtomDist <= this->get_go_cutoff(atomChainTypes[i],atomChainTypes[j]) ) {
00714             exp_a = this->get_go_exp_a(atomChainTypes[i],atomChainTypes[j]);
00715             exp_b = this->get_go_exp_b(atomChainTypes[i],atomChainTypes[j]);
00716             sigma = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
00717             goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = sigma;
00718             goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = sigma;
00719             goWithinCutoff[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = true;
00720             goWithinCutoff[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = true;
00721             nativeContacts++;
00722           } else {
00723             goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = 0.0;
00724             goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = 0.0;
00725             nonnativeContacts++;
00726           }
00727         } else {
00728           goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = -1.0;
00729           goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = -1.0;
00730         }
00731       } 
00732     }
00733   }
00734 
00735   iout << iINFO << "Number of UNIQUE    native contacts: " << nativeContacts << "\n" << endi;
00736   iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts << "\n" << endi;
00737   
00738   //  If we had to create a PDB object, delete it now
00739   if (goCoordFile != NULL) {
00740     delete goPDB;
00741   }
00742   
00743   return;
00744 }

void Molecule::build_go_sigmas2 ( StringList goCoordFile,
char *  cwd 
)

Definition at line 747 of file GoMolecule.C.

References go_pair::A, ResizeArray< Elem >::add(), PDB::atom(), atomChainTypes, atoms_1to4(), go_pair::B, ResizeArray< Elem >::begin(), StringList::data, DebugM, ResizeArray< Elem >::end(), endi(), get_go_cutoff(), get_go_exp_a(), get_go_exp_b(), go_restricted(), go_pair::goIndxA, go_pair::goIndxB, goIndxLJA, goIndxLJB, goNumLJPair, goPDB, goResidIndices, goSigmaIndices, goSigmaPairA, goSigmaPairB, iINFO(), iout, j, NAMD_die(), StringList::next, PDB::num_atoms(), numAtoms, numGoAtoms, numLJPair, pointerToGoBeg, pointerToGoEnd, and sort.

Referenced by NamdState::loadStructure().

00749 {
00750   DebugM(3,"->build_go_sigmas2" << std::endl);
00751   PDB *goPDB;      //  Pointer to PDB object to use
00752   int bcol = 4;      //  Column that data is in
00753   int32 chainType = 0;      //  b value from PDB file
00754   int32 residType = 0;      //  resid value from PDB file
00755   int i;         //  Loop counter
00756   int j;         //  Loop counter
00757   int resid1;    //  Residue ID for first atom
00758   int resid2;    //  Residue ID for second atom
00759   int residDiff;     //  Difference between resid1 and resid2
00760   Real sigma;    //  Sigma calculated for a Go atom pair
00761   Real atomAtomDist;     //  Distance between two atoms
00762   Real exp_a;            //  First exponent in L-J like formula
00763   Real exp_b;            //  Second exponent in L-J like formula
00764   char filename[129];    //  Filename
00765   
00766   //JLai
00767   int numLJPair = 0;
00768   long nativeContacts = 0;
00769   long nonnativeContacts = 0;
00770 
00771   //  Get the PDB object that contains the Go coordinates.  If
00772   //  the user gave another file name, use it.  Otherwise, just use
00773   //  the PDB file that has the initial coordinates.  
00774   if (goCoordFile == NULL)
00775     {
00776       //goPDB = initial_pdb;
00777       NAMD_die("Error: goCoordFile is NULL - build_go_sigmas2");
00778     }
00779   else
00780   {
00781     if (goCoordFile->next != NULL)
00782       {
00783         NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
00784       }
00785     
00786     if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
00787       {
00788         strcpy(filename, goCoordFile->data);
00789       }
00790     else
00791       {
00792         strcpy(filename, cwd);
00793         strcat(filename, goCoordFile->data);
00794       }
00795     
00796     goPDB = new PDB(filename);
00797     if ( goPDB == NULL )
00798       {
00799         NAMD_die("Memory allocation failed in Molecule::build_go_sigmas2");
00800       }
00801     
00802     if (goPDB->num_atoms() != numAtoms)
00803       {
00804         NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
00805       }
00806   }
00807   //  Allocate the array to hold the chain types
00808   atomChainTypes = new int32[numAtoms];
00809   //  Allocate the array to hold Go atom indices into the sigma array
00810   goSigmaIndices = new int32[numAtoms];
00811   //  Allocate the array to hold resid 
00812   goResidIndices = new int32[numAtoms];
00813 
00814   if (atomChainTypes == NULL) {
00815     NAMD_die("memory allocation failed in Molecule::build_go_sigmas2");
00816   }
00817   
00818   numGoAtoms = 0;
00819   
00820   //  Loop through all the atoms and get the Go chain types
00821   for (i=0; i<numAtoms; i++) {
00822     //  Get the chainType from the occupancy field
00823     chainType = (int32)(goPDB->atom(i))->occupancy();
00824     //  Get the resid from the resid field
00825     residType = (int32)(goPDB->atom(i))->residueseq();
00826     //  Assign the chainType value
00827     if ( chainType != 0 ) {
00828       //DebugM(3,"build_go_sigmas2 - atom:" << i << ", chainType:" << chainType << std::endl);
00829       atomChainTypes[i] = chainType;
00830       goSigmaIndices[i] = numGoAtoms;
00831       goResidIndices[i] = residType;
00832       numGoAtoms++;
00833     }
00834     else {
00835       atomChainTypes[i] = 0;
00836       goSigmaIndices[i] = -1;
00837       goResidIndices[i] = -1;
00838     }
00839   }
00840 
00841   //Define:
00842   ResizeArray<GoPair> tmpGoPair;
00843   
00844   //  Loop through all atom pairs and calculate sigmas
00845   // Store sigmas into sorted array
00846   DebugM(3,"    numAtoms=" << numAtoms << std::endl);
00847   for (i=0; i<numAtoms; i++) {
00848     resid1 = (goPDB->atom(i))->residueseq();
00849      for (j=i+1; j<numAtoms; j++) {
00850       resid2 = (goPDB->atom(j))->residueseq();
00851       if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 && !atoms_1to4(i,j) ) {
00852         residDiff = resid2 - resid1;
00853         if (residDiff < 0) residDiff = -residDiff;
00854         if ( atomChainTypes[i] && atomChainTypes[j] &&
00855              !(this->go_restricted(atomChainTypes[i],atomChainTypes[j],residDiff)) ) {
00856           atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
00857                               pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
00858                               pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
00859           if ( atomAtomDist <= this->get_go_cutoff(atomChainTypes[i],atomChainTypes[j]) ) {
00860             exp_a = this->get_go_exp_a(atomChainTypes[i],atomChainTypes[j]);
00861             exp_b = this->get_go_exp_b(atomChainTypes[i],atomChainTypes[j]);
00862             sigma = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
00863             double tmpA = pow(sigma,exp_a);
00864             double tmpB = pow(sigma,exp_b);
00865             GoPair gp;
00866             GoPair gp2;
00867             gp.goIndxA = i;
00868             gp.goIndxB = j;
00869             gp.A = tmpA;
00870             gp.B = tmpB;
00871             tmpGoPair.add(gp);
00872             gp2.goIndxA = j;
00873             gp2.goIndxB = i;
00874             gp2.A = tmpA;
00875             gp2.B = tmpB;
00876             tmpGoPair.add(gp2);
00877             nativeContacts++;
00878           } else {
00879             nonnativeContacts++;
00880           }
00881         } 
00882       } 
00883     }
00884   }
00885 
00886   iout << iINFO << "Number of UNIQUE    native contacts: " << nativeContacts << "\n" << endi;
00887   iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts << "\n" << endi;
00888 
00889   // Copies the resizeArray into a block of continuous memory
00890   std::sort(tmpGoPair.begin(),tmpGoPair.end(),goPairCompare);
00891   goNumLJPair = 2*nativeContacts;
00892   goIndxLJA = new int[goNumLJPair];
00893   goIndxLJB = new int[goNumLJPair];
00894   goSigmaPairA = new double[goNumLJPair];
00895   goSigmaPairB = new double[goNumLJPair];
00896   for(i=0; i< goNumLJPair; i++) {
00897     goIndxLJA[i] = tmpGoPair[i].goIndxA;
00898     goIndxLJB[i] = tmpGoPair[i].goIndxB;
00899     goSigmaPairA[i] = tmpGoPair[i].A;
00900     goSigmaPairB[i] = tmpGoPair[i].B;
00901   }
00902 
00903   pointerToGoBeg = new int[numAtoms];
00904   pointerToGoEnd = new int[numAtoms];
00905   int oldIndex = -1;
00906   for(i=0; i<numAtoms; i++) {
00907     pointerToGoBeg[i] = -1;
00908     pointerToGoEnd[i] = -2;
00909   }
00910   for(i=0; i < goNumLJPair; i++) {
00911     if(pointerToGoBeg[goIndxLJA[i]] == -1) {
00912       pointerToGoBeg[goIndxLJA[i]] = i;
00913       oldIndex = goIndxLJA[i];
00914     }
00915     pointerToGoEnd[oldIndex] = i;
00916   }
00917 
00918   //  If we had to create a PDB object, delete it now
00919   if (goCoordFile != NULL) {
00920     delete goPDB;
00921   }
00922     
00923   return;
00924 }

void Molecule::build_gridforce_params ( StringList gridfrcfile,
StringList gridfrccol,
StringList gridfrcchrgcol,
StringList potfile,
PDB initial_pdb,
char *  cwd 
)

Definition at line 6171 of file Molecule.C.

References PDB::atom(), DebugM, endi(), MGridforceParamsList::get_first(), GridforceGrid::get_total_grids(), MGridforceParams::gridforceCol, MGridforceParams::gridforceFile, MGridforceParams::gridforceQcol, MGridforceParams::gridforceVfile, iout, iWARN(), SimParameters::mgridforcelist, NAMD_die(), GridforceGrid::new_grid(), MGridforceParams::next, PDB::num_atoms(), numGridforceGrids, and numGridforces.

Referenced by NamdState::loadStructure().

06177 {
06178     PDB *kPDB;
06179     register int i;             //  Loop counters
06180     register int j;
06181     register int k;
06182 
06183     DebugM(3,  "Entered build_gridforce_params multi...\n");
06184 //     DebugM(3, "\tgridfrcfile = " << gridfrcfile->data << endi);
06185 //     DebugM(3, "\tgridfrccol = " << gridfrccol->data << endi);
06186     
06187     MGridforceParams* mgridParams = simParams->mgridforcelist.get_first();
06188     numGridforceGrids = 0;
06189     while (mgridParams != NULL) {
06190         numGridforceGrids++;
06191         mgridParams = mgridParams->next;
06192     }
06193     
06194     DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n");
06195     gridfrcIndexes = new int32*[numGridforceGrids];
06196     gridfrcParams = new GridforceParams*[numGridforceGrids];
06197     gridfrcGrid = new GridforceGrid*[numGridforceGrids];
06198     numGridforces = new int[numGridforceGrids];
06199     
06200     int grandTotalGrids = 0;    // including all subgrids
06201     
06202     mgridParams = simParams->mgridforcelist.get_first();
06203     for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
06204         int current_index=0;    //  Index into values used
06205         int kcol = 5;           //  Column to look for force constant in
06206         int qcol = 0;           //  Column for charge (default 0: use electric charge)
06207         Real kval = 0;          //  Force constant value retreived
06208         char filename[129];     //  PDB filename
06209         char potfilename[129];  //  Potential file name
06210         
06211         if (mgridParams == NULL) {
06212             NAMD_die("Problem with mgridParams!");
06213         }
06214         
06215         // Now load values from mgridforcelist object
06216         if (mgridParams->gridforceFile == NULL)
06217         {
06218             if ( ! initial_pdb ) NAMD_die("Initial PDB file unavailable, gridforceFile required.");
06219             kPDB = initial_pdb;
06220         }
06221         else
06222         {
06223             DebugM(4, "mgridParams->gridforceFile = " << mgridParams->gridforceFile << "\n" << endi);
06224             
06225             if ( (cwd == NULL) || (mgridParams->gridforceFile[0] == '/') )
06226             {
06227                 strcpy(filename, mgridParams->gridforceFile);
06228             }
06229             else
06230             {
06231                 strcpy(filename, cwd);
06232                 strcat(filename, mgridParams->gridforceFile);
06233             }
06234         
06235             kPDB = new PDB(filename);
06236             if ( kPDB == NULL )
06237             {
06238                 NAMD_die("Memory allocation failed in Molecule::build_gridforce_params");
06239             }
06240            
06241             if (kPDB->num_atoms() != numAtoms)
06242             {
06243                 NAMD_die("Number of atoms in grid force PDB doesn't match coordinate PDB");
06244             }
06245         }
06246 
06247         //  Get the column that the force constant is going to be in.  It
06248         //  can be in any of the 5 floating point fields in the PDB, according
06249         //  to what the user wants.  The allowable fields are X, Y, Z, O, or
06250         //  B which correspond to the 1st, 2nd, ... 5th floating point fields.
06251         //  The default is the 5th field, which is beta (temperature factor)
06252         if (mgridParams->gridforceCol == NULL)
06253         {
06254             kcol = 5;
06255         }
06256         else
06257         {
06258             if (strcasecmp(mgridParams->gridforceCol, "X") == 0)
06259             {
06260                 kcol=1;
06261             }
06262             else if (strcasecmp(mgridParams->gridforceCol, "Y") == 0)
06263             {
06264                 kcol=2;
06265             }
06266             else if (strcasecmp(mgridParams->gridforceCol, "Z") == 0)
06267             {
06268                 kcol=3;
06269             }
06270             else if (strcasecmp(mgridParams->gridforceCol, "O") == 0)
06271             {
06272                 kcol=4;
06273             }
06274             else if (strcasecmp(mgridParams->gridforceCol, "B") == 0)
06275             {
06276                 kcol=5;
06277             }
06278             else
06279             {
06280                 NAMD_die("gridforcecol must have value of X, Y, Z, O, or B");
06281             }
06282         }
06283     
06284         //  Get the column that the charge is going to be in.
06285         if (mgridParams->gridforceQcol == NULL)
06286         {
06287             qcol = 0;   // Default: don't read charge from file, use electric charge
06288         }
06289         else
06290         {
06291             if (strcasecmp(mgridParams->gridforceQcol, "X") == 0)
06292             {
06293                 qcol=1;
06294             }
06295             else if (strcasecmp(mgridParams->gridforceQcol, "Y") == 0)
06296             {
06297                 qcol=2;
06298             }
06299             else if (strcasecmp(mgridParams->gridforceQcol, "Z") == 0)
06300             {
06301                 qcol=3;
06302             }
06303             else if (strcasecmp(mgridParams->gridforceQcol, "O") == 0)
06304             {
06305                 qcol=4;
06306             }
06307             else if (strcasecmp(mgridParams->gridforceQcol, "B") == 0)
06308             {
06309                 qcol=5;
06310             }
06311             else
06312             {
06313                 NAMD_die("gridforcechargecol must have value of X, Y, Z, O, or B");
06314             }
06315         }
06316     
06317         if (kcol == qcol) {
06318             NAMD_die("gridforcecol and gridforcechargecol cannot have same value");
06319         }
06320 
06321     
06322         //  Allocate an array that will store an index into the constraint
06323         //  parameters for each atom.  If the atom is not constrained, its
06324         //  value will be set to -1 in this array.
06325         gridfrcIndexes[gridnum] = new int32[numAtoms];
06326        
06327         if (gridfrcIndexes[gridnum] == NULL)
06328         {
06329             NAMD_die("memory allocation failed in Molecule::build_gridforce_params()");
06330         }
06331         
06332         //  Loop through all the atoms and find out which ones are constrained
06333         for (i=0; i<numAtoms; i++)
06334         {
06335             //  Get the k value based on where we were told to find it
06336             switch (kcol)
06337             {
06338             case 1:
06339                 kval = (kPDB->atom(i))->xcoor();
06340                 break;
06341             case 2:
06342                 kval = (kPDB->atom(i))->ycoor();
06343                 break;
06344             case 3:
06345                 kval = (kPDB->atom(i))->zcoor();
06346                 break;
06347             case 4:
06348                 kval = (kPDB->atom(i))->occupancy();
06349                 break;
06350             case 5:
06351                 kval = (kPDB->atom(i))->temperaturefactor();
06352                 break;
06353             }
06354            
06355             if (kval > 0.0)
06356             {
06357                 //  This atom is constrained
06358                 gridfrcIndexes[gridnum][i] = current_index;
06359                 current_index++;
06360             }
06361             else
06362             {
06363                 //  This atom is not constrained
06364                 gridfrcIndexes[gridnum][i] = -1;
06365             }
06366         }
06367     
06368         if (current_index == 0)
06369         {
06370             //  Constraints were turned on, but there weren't really any constrained
06371             iout << iWARN << "NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" << endi;
06372         }
06373         else
06374         {
06375             //  Allocate an array to hold the constraint parameters
06376             gridfrcParams[gridnum] = new GridforceParams[current_index];
06377             if (gridfrcParams[gridnum] == NULL)
06378             {
06379                 NAMD_die("memory allocation failed in Molecule::build_gridforce_params");
06380             }
06381         }
06382     
06383         numGridforces[gridnum] = current_index;
06384 
06385         //  Loop through all the atoms and assign the parameters for those
06386         //  that are constrained
06387         for (i=0; i<numAtoms; i++)
06388         {
06389             if (gridfrcIndexes[gridnum][i] != -1)
06390             {
06391                 //  This atom has grid force, so get the k value again
06392                 switch (kcol)
06393                 {
06394                 case 1:
06395                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->xcoor();
06396                     break;
06397                 case 2:
06398                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->ycoor();
06399                     break;
06400                 case 3:
06401                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->zcoor();
06402                     break;
06403                 case 4:
06404                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->occupancy();
06405                     break;
06406                 case 5:
06407                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->temperaturefactor();
06408                     break;
06409                 }
06410             
06411                 //  Also get charge column
06412                 switch (qcol)
06413                 {
06414                 case 0:
06415 #ifdef MEM_OPT_VERSION
06416                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]];
06417 #else
06418                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge;
06419 #endif
06420                     break;
06421                 case 1:
06422                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->xcoor();
06423                     break;
06424                 case 2:
06425                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->ycoor();
06426                     break;
06427                 case 3:
06428                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->zcoor();
06429                     break;
06430                 case 4:
06431                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->occupancy();
06432                     break;
06433                 case 5:
06434                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->temperaturefactor();
06435                     break;
06436                 }
06437             }
06438         }
06439        
06440         //  If we had to create new PDB objects, delete them now
06441         if (mgridParams->gridforceFile != NULL)
06442         {
06443             delete kPDB;
06444         }
06445     
06446         //  Now we fill in our grid information
06447     
06448         // Open potential file
06449         if ( (cwd == NULL) || (mgridParams->gridforceVfile[0] == '/') )
06450         {
06451             strcpy(potfilename, mgridParams->gridforceVfile);
06452         }
06453         else
06454         {
06455             strcpy(potfilename, cwd);
06456             strcat(potfilename, mgridParams->gridforceVfile);
06457         }
06458     
06459 //        iout << iINFO << "Allocating grid " << gridnum
06460 //             << "\n" << endi;
06461         
06462         DebugM(3, "allocating GridforceGrid(" << gridnum << ")\n");
06463         gridfrcGrid[gridnum] = GridforceGrid::new_grid(gridnum, potfilename, simParams, mgridParams);
06464         
06465         grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
06466         DebugM(4, "grandTotalGrids = " << grandTotalGrids << "\n" << endi);
06467         
06468         // Finally, get next mgridParams pointer
06469         mgridParams = mgridParams->next;
06470     }
06471 }

void Molecule::build_gro_pair (  ) 
void Molecule::build_langevin_params ( StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_langevin_params ( BigReal  coupling,
BigReal  drudeCoupling,
Bool  doHydrogen 
)
void Molecule::build_movdrag_params ( StringList ,
StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_rotdrag_params ( StringList ,
StringList ,
StringList ,
StringList ,
StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_ss_flags ( const StringList ssfile,
const StringList sscol,
PDB initial_pdb,
const char *  cwd 
)

Build the flags needed for solute scaling.

A PDB file is read, indicating which atoms are to be scaled, and an array is maintained marking which are to be included. Each marked atom then has its corresponding van der Waals type number reassigned to enable extending the LJTable with scaled interaction values.

Parameters:
ssfile config "soluteScalingFile = my.pdb" for PDB filename
sscol config "soluteScalingCol = O" for column of PDB ATOM records
initial_pdb the initial PDB file
cwd current working directory

Referenced by NamdState::loadStructure().

void Molecule::build_stirred_atoms ( StringList ,
StringList ,
PDB ,
char *   
)
int Molecule::checkexcl ( int  atom1,
int  atom2 
) const
void Molecule::compute_LJcorrection (  ) 
void Molecule::delete_alch_bonded ( void   ) 
void Molecule::delete_qm_bonded ( void   ) 

Definition at line 1547 of file MoleculeQM.C.

References crossterm::atom1, improper::atom1, dihedral::atom1, angle::atom1, bond::atom1, DebugM, endi(), SimParameters::extraBondsOn, iERROR(), iINFO(), iout, NAMD_die(), numAngles, numBonds, numCrossterms, numDihedrals, numImpropers, numRealBonds, and ResizeArray< Elem >::size().

01547                                      {
01548 
01549 #ifdef MEM_OPT_VERSION
01550     NAMD_die("QMMM interface is not supported in memory optimized builds.");
01551 #else
01552     
01553     DebugM(3,"Cleaning QM bonds, angles, dihedrals, impropers and crossterms.\n")
01554     
01555     // Bonds
01556     Bond *nonQMBonds;
01557     nonQMBonds = new Bond[numBonds] ;
01558     int nonQMBondCount = 0;
01559     qmDroppedBonds = 0;
01560     for (int i = 0; i < numBonds; i++) {
01561         
01562         int part1 = qmAtomGroup[bonds[i].atom1];
01563         int part2 = qmAtomGroup[bonds[i].atom2];
01564         if (part1 > 0 && part2 > 0 ) {
01565             
01566             // If the user defined extra bonds, we check if the QM-QM bond is an "extra"
01567             // bond, and do not delete it.
01568             if (simParams->extraBondsOn) {
01569 //                 std::cout << "Checking Bond: " << bonds[i].atom1 << "," << bonds[i].atom2 << std::endl ;
01570                 
01571                 for (int ebi=0; ebi < qmExtraBonds.size() ; ebi++) {
01572                     
01573                     if( (qmExtraBonds[ebi].atom1 == bonds[i].atom1 || 
01574                         qmExtraBonds[ebi].atom1 == bonds[i].atom2) && 
01575                     (qmExtraBonds[ebi].atom2 == bonds[i].atom1 || 
01576                         qmExtraBonds[ebi].atom2 == bonds[i].atom2) ) {
01577 //                         std::cout << "This is an extra bond! We will keep it." << std::endl;
01578                         nonQMBonds[nonQMBondCount++] = bonds[i];
01579                         break;
01580                     }
01581                 }
01582             }
01583             
01584             qmDroppedBonds++;
01585         } else {
01586             // Just a simple sanity check.
01587             // If there are no QM bonds defined by the user but we find a
01588             // bond between a QM and an MM atom, error out.
01589             if ((part1 > 0 || part2 > 0) && qmNumBonds == 0) {
01590                 iout << iERROR << " Atoms " << bonds[i].atom1
01591                 << " and " << bonds[i].atom2 << " form a QM-MM bond!\n" << endi ;
01592                 NAMD_die("No QM-MM bonds were defined by the user, but one was found.");
01593             }
01594             nonQMBonds[nonQMBondCount++] = bonds[i];
01595         }
01596     }
01597     numBonds = nonQMBondCount;
01598     delete [] bonds;
01599     bonds = new Bond[numBonds];
01600     for (int i = 0; i < nonQMBondCount; i++) {
01601         bonds[i]=nonQMBonds[i];
01602     }
01603     delete [] nonQMBonds;
01604   numRealBonds = numBonds;
01605     
01606   // Angles
01607   Angle *nonQMAngles;
01608   nonQMAngles = new Angle[numAngles];
01609   int nonQMAngleCount = 0;
01610   qmDroppedAngles = 0;
01611   for (int i = 0; i < numAngles; i++) {
01612     int part1 = qmAtomGroup[angles[i].atom1];
01613     int part2 = qmAtomGroup[angles[i].atom2];
01614     int part3 = qmAtomGroup[angles[i].atom3];
01615     if (part1 > 0 && part2 > 0 && part3 > 0) {
01616       //CkPrintf("-----ANGLE ATOMS %i %i %i partitions %i %i %i\n",angles[i].atom1, angles[i].atom2, angles[i].atom3, part1, part2, part3);
01617       qmDroppedAngles++;
01618     }
01619     else {
01620       nonQMAngles[nonQMAngleCount++] = angles[i];
01621     }
01622   }
01623   numAngles = nonQMAngleCount;
01624   delete [] angles;
01625   angles = new Angle[numAngles];
01626   for (int i = 0; i < nonQMAngleCount; i++) {
01627     angles[i]=nonQMAngles[i];
01628   }
01629   delete [] nonQMAngles;
01630 
01631 
01632   // Dihedrals
01633   Dihedral *nonQMDihedrals;
01634   nonQMDihedrals = new Dihedral[numDihedrals];
01635   int nonQMDihedralCount = 0;
01636   qmDroppedDihedrals = 0;
01637   for (int i = 0; i < numDihedrals; i++) {
01638     int part1 = qmAtomGroup[dihedrals[i].atom1];
01639     int part2 = qmAtomGroup[dihedrals[i].atom2];
01640     int part3 = qmAtomGroup[dihedrals[i].atom3];
01641     int part4 = qmAtomGroup[dihedrals[i].atom4];
01642     if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0) {
01643       //CkPrintf("-----i %i DIHEDRAL ATOMS %i %i %i %i partitions %i %i %i %i\n",i,dihedrals[i].atom1, dihedrals[i].atom2, dihedrals[i].atom3, dihedrals[i].atom4, part1, part2, part3,part4);
01644       qmDroppedDihedrals++;
01645     }
01646     else {
01647       nonQMDihedrals[nonQMDihedralCount++] = dihedrals[i];
01648     }
01649   }
01650   numDihedrals = nonQMDihedralCount;
01651   delete [] dihedrals;
01652   dihedrals = new Dihedral[numDihedrals];
01653   for (int i = 0; i < numDihedrals; i++) {
01654     dihedrals[i]=nonQMDihedrals[i];
01655   }
01656   delete [] nonQMDihedrals;
01657 
01658   // Impropers
01659   Improper *nonQMImpropers;
01660   nonQMImpropers = new Improper[numImpropers];
01661   int nonQMImproperCount = 0;
01662   qmDroppedImpropers = 0;
01663   for (int i = 0; i < numImpropers; i++) {
01664     int part1 = qmAtomGroup[impropers[i].atom1];
01665     int part2 = qmAtomGroup[impropers[i].atom2];
01666     int part3 = qmAtomGroup[impropers[i].atom3];
01667     int part4 = qmAtomGroup[impropers[i].atom4];
01668     if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0) {
01669       //CkPrintf("-----i %i IMPROPER ATOMS %i %i %i %i partitions %i %i %i %i\n",i,impropers[i].atom1, impropers[i].atom2, impropers[i].atom3, impropers[i].atom4, part1, part2, part3,part4);
01670       qmDroppedImpropers++;
01671     }
01672     else {
01673       nonQMImpropers[nonQMImproperCount++] = impropers[i];
01674     }
01675   }
01676   numImpropers = nonQMImproperCount;
01677   delete [] impropers;
01678   impropers = new Improper[numImpropers];
01679   for (int i = 0; i < numImpropers; i++) {
01680     impropers[i]=nonQMImpropers[i];
01681   }
01682   delete [] nonQMImpropers;
01683   
01684   // Crossterms 
01685   Crossterm *nonQMCrossterms;
01686   nonQMCrossterms = new Crossterm[numCrossterms];
01687   int nonQMCrosstermCount = 0;
01688   qmDroppedCrossterms = 0 ;
01689   for (int i = 0; i < numCrossterms; i++) {
01690     int part1 = qmAtomGroup[crossterms[i].atom1];
01691     int part2 = qmAtomGroup[crossterms[i].atom2];
01692     int part3 = qmAtomGroup[crossterms[i].atom3];
01693     int part4 = qmAtomGroup[crossterms[i].atom4];
01694     int part5 = qmAtomGroup[crossterms[i].atom5];
01695     int part6 = qmAtomGroup[crossterms[i].atom6];
01696     int part7 = qmAtomGroup[crossterms[i].atom7];
01697     int part8 = qmAtomGroup[crossterms[i].atom8];
01698     if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0 &&
01699         part5 > 0 && part6 > 0 && part7 > 0 && part8 > 0) {
01700       //CkPrintf("-----i %i IMPROPER ATOMS %i %i %i %i partitions %i %i %i %i\n",i,impropers[i].atom1, impropers[i].atom2, impropers[i].atom3, impropers[i].atom4, part1, part2, part3,part4);
01701       qmDroppedCrossterms++;
01702     }
01703     else {
01704       nonQMCrossterms[nonQMCrosstermCount++] = crossterms[i];
01705     }
01706   }
01707   numCrossterms = nonQMCrosstermCount;
01708   delete [] crossterms;
01709   crossterms = new Crossterm[numCrossterms] ;
01710   for (int i=0; i < numCrossterms; i++){
01711       crossterms[i] = nonQMCrossterms[i] ;
01712   }
01713   delete [] nonQMCrossterms ;
01714   
01715   if (!CkMyPe()) {
01716       iout << iINFO << "The QM region will remove " << qmDroppedBonds << " bonds, " << 
01717           qmDroppedAngles << " angles, " << qmDroppedDihedrals << " dihedrals, "
01718           << qmDroppedImpropers << " impropers and " << qmDroppedCrossterms 
01719           << " crossterms.\n" << endi ;
01720   }
01721   
01722 #endif
01723 }

void Molecule::freeBFactorData (  )  [inline]

Definition at line 1021 of file Molecule.h.

Referenced by NamdState::loadStructure().

01021 { delete [] bfactor; bfactor=NULL; }

void Molecule::freeOccupancyData (  )  [inline]

Definition at line 1017 of file Molecule.h.

Referenced by NamdState::loadStructure().

01017 { delete [] occupancy; occupancy=NULL; }

Bond* Molecule::get_acceptor ( int  dnum  )  const [inline]

Definition at line 1092 of file Molecule.h.

01092 {return (&(acceptors[dnum]));} 

Angle* Molecule::get_angle ( int  anum  )  const [inline]

Definition at line 1055 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01055 {return (&(angles[anum]));}

int32* Molecule::get_angles_for_atom ( int  anum  )  [inline]

Definition at line 1134 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01135       { return anglesByAtom[anum]; }

int Molecule::get_atom_from_index_in_residue ( const char *  segid,
int  resid,
int  index 
) const

Definition at line 158 of file Molecule.C.

References ResidueLookupElem::lookup(), and NAMD_die().

Referenced by GlobalMasterFreeEnergy::getAtomID(), and GlobalMasterEasy::getAtomID().

00159                                                        {
00160 
00161   if (atomNames == NULL || resLookup == NULL)
00162   {
00163     NAMD_die("Tried to find atom from name on node other than node 0");
00164   }
00165   int i = 0;
00166   int end = 0;
00167   if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
00168   if ( index >= 0 && index < ( end - i ) ) return ( index + i );
00169   return -1;
00170 }

int Molecule::get_atom_from_name ( const char *  segid,
int  resid,
const char *  aname 
) const

Definition at line 121 of file Molecule.C.

References atomNamePool, ResidueLookupElem::lookup(), and NAMD_die().

Referenced by colvarproxy_namd::check_atom_id(), GlobalMasterFreeEnergy::getAtomID(), and GlobalMasterEasy::getAtomID().

00122                                                                {
00123 
00124   if (atomNames == NULL || resLookup == NULL)
00125   {
00126     NAMD_die("Tried to find atom from name on node other than node 0");
00127   }
00128 
00129   int i = 0;
00130   int end = 0;
00131   if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
00132   for ( ; i < end; ++i ) {
00133     #ifdef MEM_OPT_VERSION    
00134     Index idx = atomNames[i].atomnameIdx;
00135     if(!strcasecmp(aname, atomNamePool[idx])) return i;
00136     #else
00137     if ( ! strcasecmp(aname,atomNames[i].atomname) ) return i;
00138     #endif
00139   }
00140   return -1;
00141 }

const char* Molecule::get_atomtype ( int  anum  )  const [inline]

Definition at line 1103 of file Molecule.h.

References atomTypePool, and NAMD_die().

01104   {
01105     if (atomNames == NULL)
01106     {
01107       NAMD_die("Tried to find atom type on node other than node 0");
01108     }
01109 
01110     #ifdef MEM_OPT_VERSION    
01111     return atomTypePool[atomNames[anum].atomtypeIdx];
01112     #else
01113     return(atomNames[anum].atomtype);
01114     #endif
01115   }

Bond* Molecule::get_bond ( int  bnum  )  const [inline]

Definition at line 1052 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01052 {return (&(bonds[bnum]));}

int32* Molecule::get_bonds_for_atom ( int  anum  )  [inline]

Definition at line 1132 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01133       { return bondsByAtom[anum]; } 

int Molecule::get_cluster ( int  anum  )  const [inline]

Definition at line 1010 of file Molecule.h.

Referenced by wrap_coor_int().

01010 { return cluster[anum]; }

int Molecule::get_clusterSize ( int  anum  )  const [inline]

Definition at line 1011 of file Molecule.h.

Referenced by wrap_coor_int().

01011 { return clusterSize[anum]; }

void Molecule::get_cons_params ( Real k,
Vector refPos,
int  atomnum 
) const [inline]

Definition at line 1249 of file Molecule.h.

Referenced by ComputeRestraints::doForce().

01250   {
01251     k = consParams[consIndexes[atomnum]].k;
01252     refPos = consParams[consIndexes[atomnum]].refPos;
01253   }

void Molecule::get_constorque_params ( BigReal v,
Vector a,
Vector p,
int  atomnum 
) const [inline]

Definition at line 1323 of file Molecule.h.

References consTorqueIndexes, and consTorqueParams.

Referenced by ComputeConsTorque::doForce().

01325   {
01326     v = consTorqueParams[consTorqueIndexes[atomnum]].v;
01327     a = consTorqueParams[consTorqueIndexes[atomnum]].a;
01328     p = consTorqueParams[consTorqueIndexes[atomnum]].p;
01329   }

Crossterm* Molecule::get_crossterm ( int  inum  )  const [inline]

Definition at line 1064 of file Molecule.h.

01064 {return (&(crossterms[inum]));}

int32* Molecule::get_crossterms_for_atom ( int  anum  )  [inline]

Definition at line 1140 of file Molecule.h.

01141       { return crosstermsByAtom[anum]; }  

const Real* Molecule::get_cSMDcoffs (  )  [inline]

Definition at line 866 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00866 {return cSMDcoffs;} ;

int const* const* Molecule::get_cSMDindex (  )  [inline]

Definition at line 862 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00862 {return cSMDindex;} ;

int const* Molecule::get_cSMDindxLen (  )  [inline]

Definition at line 861 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00861 { return cSMDindxLen;} ;

const Real* Molecule::get_cSMDKs (  )  [inline]

Definition at line 864 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00864 {return cSMDKs;} ;

int Molecule::get_cSMDnumInst (  )  [inline]

Definition at line 860 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00860 { return cSMDnumInst;} ;

int const* const* Molecule::get_cSMDpairs (  )  [inline]

Definition at line 863 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00863 {return cSMDpairs;} ;

const Real* Molecule::get_cSMDVels (  )  [inline]

Definition at line 865 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00865 {return cSMDVels;} ;

Dihedral* Molecule::get_dihedral ( int  dnum  )  const [inline]

Definition at line 1061 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01061 {return (&(dihedrals[dnum]));}

int32* Molecule::get_dihedrals_for_atom ( int  anum  )  [inline]

Definition at line 1136 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01137       { return dihedralsByAtom[anum]; }

Bond* Molecule::get_donor ( int  dnum  )  const [inline]

Definition at line 1089 of file Molecule.h.

01089 {return (&(donors[dnum]));}  

const ExclusionCheck* Molecule::get_excl_check_for_atom ( int  anum  )  const [inline]

Definition at line 1160 of file Molecule.h.

Referenced by dumpbench().

01160                                                                {      
01161       return &all_exclusions[anum];             
01162   }

Exclusion* Molecule::get_exclusion ( int  ex  )  const [inline]

Definition at line 1099 of file Molecule.h.

References exclusions.

01099 {return (&(exclusions[ex]));}

int32* Molecule::get_exclusions_for_atom ( int  anum  )  [inline]

Definition at line 1142 of file Molecule.h.

References exclusionsByAtom.

01143       { return exclusionsByAtom[anum]; }

int Molecule::get_fep_bonded_type ( const int *  atomID,
unsigned int  order 
) const [inline]

Definition at line 1372 of file Molecule.h.

References NAMD_die(), and simParams.

Referenced by ImproperElem::computeForce(), DihedralElem::computeForce(), BondElem::computeForce(), AnisoElem::computeForce(), and AngleElem::computeForce().

01373   {
01374     int typeSum = 0;
01375     for ( int i=0; i < order; ++i ) {
01376       typeSum += (fepAtomFlags[atomID[i]] == 2 ? -1 : fepAtomFlags[atomID[i]]);
01377     }
01378     // Increase the cutoff if scaling purely alchemical bonds. 
01379     // This really only applies when order = 2.
01380     if ( simParams->alchBondDecouple ) order++;
01381 
01382     // The majority of interactions are type 0, so catch those first.
01383     if ( typeSum == 0 || abs(typeSum) == order ) return 0;
01384     else if ( 0 < typeSum && typeSum < order ) return 1;
01385     else if ( -order < typeSum && typeSum < 0 ) return 2;
01386 
01387     // Alchemify should always keep this from bombing, but just in case...
01388     NAMD_die("Unexpected alchemical bonded interaction!");
01389     return 0;
01390   }

unsigned char Molecule::get_fep_type ( int  anum  )  const [inline]
const int32* Molecule::get_full_exclusions_for_atom ( int  anum  )  const [inline]

Definition at line 1144 of file Molecule.h.

Referenced by ComputeNonbondedCUDA::build_exclusions().

01145       { return fullExclusionsByAtom[anum]; }

Real Molecule::get_go_cutoff ( int  chain1,
int  chain2 
) [inline]

Definition at line 1588 of file Molecule.h.

References go_val::cutoff, go_array, and MAX_GO_CHAINS.

Referenced by build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), get_go_force(), get_go_force2(), and get_go_force_new().

01588 { return go_array[MAX_GO_CHAINS*chain1 + chain2].cutoff; };

Real Molecule::get_go_epsilon ( int  chain1,
int  chain2 
) [inline]

Definition at line 1597 of file Molecule.h.

References go_val::epsilon, go_array, and MAX_GO_CHAINS.

Referenced by get_go_force(), get_go_force2(), and get_go_force_new().

01597 { return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilon; };

Real Molecule::get_go_epsilonRep ( int  chain1,
int  chain2 
) [inline]

Definition at line 1591 of file Molecule.h.

References go_val::epsilonRep, go_array, and MAX_GO_CHAINS.

Referenced by get_go_force(), get_go_force2(), and get_go_force_new().

01591 { return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilonRep; };

int Molecule::get_go_exp_a ( int  chain1,
int  chain2 
) [inline]

Definition at line 1600 of file Molecule.h.

References go_val::exp_a, go_array, and MAX_GO_CHAINS.

Referenced by build_go_sigmas(), build_go_sigmas2(), get_go_force(), get_go_force2(), and get_go_force_new().

01600 { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_a; };

int Molecule::get_go_exp_b ( int  chain1,
int  chain2 
) [inline]

Definition at line 1603 of file Molecule.h.

References go_val::exp_b, go_array, and MAX_GO_CHAINS.

Referenced by build_go_sigmas(), build_go_sigmas2(), get_go_force(), get_go_force2(), and get_go_force_new().

01603 { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_b; };

int Molecule::get_go_exp_rep ( int  chain1,
int  chain2 
) [inline]

Definition at line 1606 of file Molecule.h.

References go_val::exp_rep, go_array, and MAX_GO_CHAINS.

Referenced by get_go_force(), get_go_force2(), and get_go_force_new().

01606 { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_rep; };

BigReal Molecule::get_go_force ( BigReal  r,
int  atom1,
int  atom2,
BigReal goNative,
BigReal goNonnative 
) const

Definition at line 1260 of file GoMolecule.C.

References atomChainTypes, get_go_cutoff(), get_go_epsilon(), get_go_epsilonRep(), get_go_exp_a(), get_go_exp_b(), get_go_exp_rep(), get_go_sigmaRep(), goSigmaIndices, goSigmas, goWithinCutoff, and numGoAtoms.

01265 {
01266   BigReal goForce = 0.0;
01267   Real pow1;
01268   Real pow2;
01269   //  determine which Go chain pair we are working with
01270   //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ")" << std::endl);
01271   int32 chain1 = atomChainTypes[atom1];
01272   int32 chain2 = atomChainTypes[atom2];
01273 
01274   //DebugM(3,"  chain1:" << chain1 << ", chain2:" << chain2 << std::endl);
01275   if (chain1 == 0 || chain2 == 0)  return 0.0;
01276 
01277   //  retrieve Go cutoff for this chain pair
01278   //TMP// JLai -- I'm going to replace this with a constant accessor.  This is just a temporary thing
01279   Real goCutoff = const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2);
01280   //DebugM(3,"  goCutoff:" << goCutoff << std::endl);
01281   if (goCutoff == 0)  return 0.0;
01282   //  if repulsive then calculate repulsive
01283   //  sigmas are initially set to -1.0 if the atom pair fails go_restricted
01284   if (goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]] != -1.0) {
01285     if (!goWithinCutoff[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]]) {
01286       Real epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1,chain2);
01287       Real sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1,chain2);
01288       int exp_rep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
01289       pow1 = pow(sigmaRep/r,exp_rep);
01290       goForce = 4*((exp_rep/(r*r)) * epsilonRep * pow1);
01291       *goNative = 0.0;
01292       *goNonnative = (4 * epsilonRep * pow1 );
01293       //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", epsilonRep:" << epsilonRep << ", sigmaRep:" << sigmaRep << ", r:" << r << ", goForce:" << goForce << std::endl);
01294     }
01295     //  if attractive then calculate attractive
01296     else {
01297       int goSigmaIndex1 = goSigmaIndices[atom1];
01298       int goSigmaIndex2 = goSigmaIndices[atom2];
01299       if (goSigmaIndex1 != -1 && goSigmaIndex2 != -1) {
01300         Real epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
01301         int exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
01302         int exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
01303         Real sigma_ij = goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]];
01304         // Positive gradient of potential, not negative gradient of potential
01305         pow1 = pow(sigma_ij/r,exp_a);
01306         pow2 = pow(sigma_ij/r,exp_b);
01307         goForce = ((4/(r*r)) * epsilon * (exp_a * pow1 - exp_b * pow2));
01308         //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", sigma_ij:" << sigma_ij << ", r:" << r << ", goForce:" << goForce << std::endl);
01309         *goNative = (4 * epsilon * ( pow1 -  pow2 ) );
01310         *goNonnative = 0.0;
01311       }
01312     }
01313   }
01314   //DebugM(3,"goForce:" << goForce << std::endl);
01315   return goForce;
01316 }

BigReal Molecule::get_go_force2 ( BigReal  x,
BigReal  y,
BigReal  z,
int  atom1,
int  atom2,
BigReal goNative,
BigReal goNonnative 
) const

Definition at line 1456 of file GoMolecule.C.

References atomChainTypes, get_go_cutoff(), get_go_epsilon(), get_go_epsilonRep(), get_go_exp_a(), get_go_exp_b(), get_go_exp_rep(), get_go_sigmaRep(), go_restricted(), goIndxLJB, goResidIndices, goSigmaPairA, goSigmaPairB, pointerToGoBeg, and pointerToGoEnd.

01463 {
01464  
01465   // Check to see if restricted.  If so, escape function early
01466   int32 chain1 = atomChainTypes[atom1];
01467   int32 chain2 = atomChainTypes[atom2];
01468 
01469   if(chain1 == 0 || chain2 == 0) return 0.0;
01470   Molecule *mol = const_cast<Molecule*>(this);
01471   Real goCutoff = mol->get_go_cutoff(chain1,chain2);
01472   if(goCutoff == 0) return 0.0;
01473 
01474   int resid1 = goResidIndices[atom1];
01475   int resid2 = goResidIndices[atom2];
01476   int residDiff = abs(resid1 - resid2);
01477   if((mol->go_restricted(chain1,chain2,residDiff))) {
01478     return 0.0;
01479   }
01480 
01481   int LJIndex = -1;
01482   int LJbegin = pointerToGoBeg[atom1];
01483   int LJend = pointerToGoEnd[atom1];
01484   for(int i = LJbegin; i <= LJend; i++) {
01485     if(goIndxLJB[i] == atom2) {
01486       LJIndex = i;
01487     }
01488   }
01489   
01490   BigReal r2 = x*x + y*y + z*z;
01491   BigReal r = sqrt(r2);
01492 
01493   if (LJIndex == -1) {
01494     int exp_rep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
01495     BigReal epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1, chain2);
01496     BigReal sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1, chain2);
01497     double sigmaRepPow = pow(sigmaRep,exp_rep);
01498     BigReal LJ = (4*epsilonRep*exp_rep*sigmaRepPow*pow(r,-(exp_rep+1)));
01499     *goNative = 0;
01500     *goNonnative = (4*epsilonRep*sigmaRepPow*pow(r,-exp_rep));
01501     //*goNonnative = (4*epsilonRep * pow(sigmaRep/r,exp_rep));
01502     return (LJ/r);
01503   } else {
01504     // Code to calculate distances because the pair was found in one of the lists
01505     int exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
01506     int exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
01507     // We want the force, so we have to take the n+1 derivative
01508     BigReal powA = pow(r,-(exp_a + 1));
01509     BigReal powB = pow(r,-(exp_b + 1));
01510     BigReal powaa = pow(r,-exp_a);
01511     BigReal powbb = pow(r,-exp_b);
01512     BigReal epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
01513     BigReal LJ = 4 * epsilon * (exp_a*goSigmaPairA[LJIndex]*powA - exp_b*goSigmaPairB[LJIndex]*powB);
01514     *goNative =  4 * epsilon * (goSigmaPairA[LJIndex]*powaa - goSigmaPairB[LJIndex]*powbb);
01515     *goNonnative = 0;
01516     return (LJ/r);
01517   }
01518   return 0;
01519 }

BigReal Molecule::get_go_force_new ( BigReal  r,
int  atom1,
int  atom2,
BigReal goNative,
BigReal goNonnative 
) const

Definition at line 1334 of file GoMolecule.C.

References atomChainTypes, DebugM, get_go_cutoff(), get_go_epsilon(), get_go_epsilonRep(), get_go_exp_a(), get_go_exp_b(), get_go_exp_rep(), get_go_sigmaRep(), go_restricted(), goCoordinates, goResids, and goSigmaIndices.

01339 {
01340   int resid1;
01341   int resid2;
01342   int residDiff;
01343   Real xcoorDiff;
01344   Real ycoorDiff;
01345   Real zcoorDiff;
01346   Real atomAtomDist;
01347   Real exp_a;
01348   Real exp_b;
01349   Real sigma_ij;
01350   Real epsilon;
01351   Real epsilonRep;
01352   Real sigmaRep;
01353   Real expRep;
01354   Real pow1;
01355   Real pow2;
01356   
01357   BigReal goForce = 0.0;
01358   *goNative = 0;
01359   *goNonnative = 0;
01360 
01361   //  determine which Go chain pair we are working with
01362   DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ")" << std::endl);
01363   int goIndex1 = goSigmaIndices[atom1];
01364   int goIndex2 = goSigmaIndices[atom2];
01365 
01366   int32 chain1 = atomChainTypes[goIndex1];
01367   int32 chain2 = atomChainTypes[goIndex2];
01368 
01369   DebugM(3,"  chain1:" << chain1 << ", chain2:" << chain2 << std::endl);
01370   if (chain1 == 0 || chain2 == 0)  return 0.0;
01371 
01372   //  retrieve Go cutoff for this chain pair
01373   Real goCutoff = const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2);
01374   DebugM(3,"  goCutoff:" << goCutoff << std::endl);
01375   if (goCutoff == 0)  return 0.0;
01376 
01377   //  sigmas are initially set to -1.0 if the atom pair fails go_restricted
01378   //  no goSigmas array anymore
01379   //Real sigma_ij = goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]];
01380 
01381   // XXX - used to be a condition for the following if
01382   //if the atoms are within 4 of each other
01383   //if ( !atoms_1to4(atom1,atom2) ) {
01384 
01385   //  if goSigmaIndices aren't defined, don't calculate forces
01386   if ( goIndex1 != -1 && goIndex2 != -1 ) {
01387     resid1 = goResids[goIndex1];
01388     resid2 = goResids[goIndex2];
01389     residDiff = resid2 - resid1;
01390     if (residDiff < 0) residDiff = -residDiff;
01391     //  if this is a Go atom pair that is not restricted
01392     if ( !(const_cast<Molecule*>(this)->go_restricted(chain1,chain2,residDiff)) ) {
01393       xcoorDiff = goCoordinates[goIndex1*3] - goCoordinates[goIndex2*3];
01394       ycoorDiff = goCoordinates[goIndex1*3 + 1] - goCoordinates[goIndex2*3 + 1];
01395       zcoorDiff = goCoordinates[goIndex1*3 + 2] - goCoordinates[goIndex2*3 + 2];
01396       atomAtomDist = sqrt(xcoorDiff*xcoorDiff + ycoorDiff*ycoorDiff + zcoorDiff*zcoorDiff);
01397       
01398       //  if attractive then calculate attractive
01399       if ( atomAtomDist <= const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2) ) {
01400         exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
01401         exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
01402         sigma_ij = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
01403         
01404         // [JLai] print out atoms involved in native contacts
01405         // printf("ATOM1: %d C1: %d ATOM2: %d C2: %d\n", atom1,chain1,atom2,chain2);
01406 
01407         epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
01408         pow1 = pow(sigma_ij/r,static_cast<double>(exp_a));
01409         pow2 = pow(sigma_ij/r,static_cast<double>(exp_b));
01410         //goForce = ((4/r) * epsilon * (exp_a * pow(sigma_ij/r,exp_a) - exp_b * pow(sigma_ij/r,exp_b)));
01411         goForce = ((4/(r*r)) * epsilon * (exp_a * pow1 - exp_b * pow2));
01412         DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", exp_a:" << exp_a << ", exp_b:" << exp_b << ", sigma_ij:" << sigma_ij << ", r:" << r << ", goForce:" << goForce << std::endl);
01413         //goEnergy = (4 * epsilon * ( pow(sigma_ij/r,exp_a) -  pow(sigma_ij/r,exp_b) ) ); // JLai I changed some of the expressions
01414         *goNative = (4 * epsilon * ( pow1 -  pow2 ) ); 
01415         *goNonnative = 0;
01416       }
01417       
01418       //  if repulsive then calculate repulsive
01419       else {
01420         epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1,chain2);
01421         sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1,chain2);
01422         expRep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
01423         pow1 = pow(sigmaRep/r,(BigReal)expRep);
01424         //goForce = ((12.0/r) * epsilonRep * pow(sigmaRep/r,12.0));
01425         goForce = ((4/(r*r)) * expRep * epsilonRep * pow1);
01426         DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", epsilonRep:" << epsilonRep << ", sigmaRep:" << sigmaRep << ", r:" << r << ", goForce:" << goForce << std::endl);
01427         //goEnergy = (4 * epsilonRep * pow(sigmaRep/r,12.0)); // JLai I changed some of the expressions
01428         *goNonnative = (4 * epsilonRep * pow1); 
01429         *goNative = 0;
01430       }
01431     }
01432   }
01433   
01434   //DebugM(3,"goForce:" << goForce << std::endl);
01435   return goForce;
01436 }

Real Molecule::get_go_sigmaRep ( int  chain1,
int  chain2 
) [inline]

Definition at line 1594 of file Molecule.h.

References go_array, MAX_GO_CHAINS, and go_val::sigmaRep.

Referenced by get_go_force(), get_go_force2(), and get_go_force_new().

01594 { return go_array[MAX_GO_CHAINS*chain1 + chain2].sigmaRep; };

GridforceGrid* Molecule::get_gridfrc_grid ( int  gridnum  )  const [inline]

Definition at line 1263 of file Molecule.h.

References numGridforceGrids.

Referenced by ComputeGridForce::doForce(), Node::reloadGridforceGrid(), and Node::updateGridScale().

01264   {
01265       GridforceGrid *result = NULL;
01266       if (gridnum >= 0 && gridnum < numGridforceGrids) {
01267           result = gridfrcGrid[gridnum];
01268       }
01269       return result;
01270   }

void Molecule::get_gridfrc_params ( Real k,
Charge q,
int  atomnum,
int  gridnum 
) const [inline]

Definition at line 1257 of file Molecule.h.

Referenced by ComputeGridForce::do_calc().

01258   {
01259       k = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].k;
01260       q = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].q;
01261   }

BigReal Molecule::get_gro_force ( BigReal  ,
BigReal  ,
BigReal  ,
int  ,
int   
) const
BigReal Molecule::get_gro_force2 ( BigReal  x,
BigReal  y,
BigReal  z,
int  atom1,
int  atom2,
BigReal pairLJEnergy,
BigReal pairGaussEnergy 
) const

Definition at line 1173 of file GoMolecule.C.

References A, gA, giSigma1, giSigma2, gMu1, gMu2, gRepulsive, indxGaussB, indxLJB, pairC12, pairC6, pointerToGaussBeg, pointerToGaussEnd, pointerToLJBeg, and pointerToLJEnd.

01180 {
01181   //Initialize return energies to zero
01182   *pairLJEnergy = 0.0;
01183   *pairGaussEnergy = 0.0;
01184 
01185   // Linear search for LJ data
01186   int LJIndex = -1;
01187   int LJbegin = pointerToLJBeg[atom1];
01188   int LJend = pointerToLJEnd[atom1];
01189   for(int i = LJbegin; i <= LJend; i++) {
01190     if(indxLJB[i] == atom2) {
01191       LJIndex = i;
01192       break;
01193     }
01194   }
01195 
01196   // Linear search for Gaussian data
01197   int GaussIndex = -1;
01198   int Gaussbegin = pointerToGaussBeg[atom1];
01199   int Gaussend = pointerToGaussEnd[atom1];
01200   for(int i = Gaussbegin; i <= Gaussend; i++) {
01201     if(indxGaussB[i] == atom2) {
01202       GaussIndex = i;
01203       break;
01204     }
01205   }
01206 
01207   if( LJIndex == -1 && GaussIndex == -1) {
01208     return 0;
01209   } else {
01210     // Code to calculate distances because the pair was found in one of the lists
01211     BigReal r2 = x*x + y*y + z*z;
01212     BigReal r = sqrt(r2);
01213     BigReal ri = 1/r;
01214     BigReal ri6 = (ri*ri*ri*ri*ri*ri);
01215     BigReal ri12 = ri6*ri6;
01216     BigReal ri13 = ri12*ri;
01217     BigReal LJ = 0;
01218     BigReal Gauss = 0;
01219     // Code to calculate LJ 
01220     if (LJIndex != -1) {
01221       BigReal ri7 = ri6*ri;
01222       LJ = (12*(pairC12[LJIndex]*ri13) - 6*(pairC6[LJIndex]*ri7));
01223       *pairLJEnergy = (pairC12[LJIndex]*ri12 - pairC6[LJIndex]*ri6);
01224       //std::cout << pairC12[LJIndex] << " " << pairC6[LJIndex] << " " << ri13 << " " << ri7 << " " << LJ << " " << r << "\n";
01225     }
01226     // Code to calculate Gaussian
01227     if (GaussIndex != -1) {
01228       BigReal gr = 12*gRepulsive[GaussIndex]*ri13;
01229       BigReal r1prime = r - gMu1[GaussIndex];
01230       BigReal tmp1 = r1prime * r1prime;
01231       BigReal r2prime = r - gMu2[GaussIndex];
01232       BigReal tmp2 = r2prime * r2prime;
01233       BigReal tmp_gauss1 = 0;
01234       BigReal one_gauss1 = 1;
01235       BigReal tmp_gauss2 = 0;
01236       BigReal one_gauss2 = 1;
01237       if (giSigma1[GaussIndex] != 0) {
01238         tmp_gauss1 = exp(-tmp1*giSigma1[GaussIndex]);
01239         one_gauss1 = 1 - tmp_gauss1;
01240       }
01241       if (giSigma2[GaussIndex] != 0) {
01242         tmp_gauss2 = exp(-tmp2*giSigma2[GaussIndex]);
01243         one_gauss2 = 1 - tmp_gauss2;
01244       } 
01245       BigReal A = gA[GaussIndex];
01246       Gauss = gr*one_gauss1*one_gauss2 - 2*A*tmp_gauss1*one_gauss2*r1prime*giSigma1[GaussIndex] \
01247         - 2*tmp_gauss1*one_gauss2*r1prime*giSigma1[GaussIndex]*gRepulsive[GaussIndex]*ri12 - 2*A*tmp_gauss2*one_gauss1*r2prime*giSigma2[GaussIndex] \
01248         - 2*tmp_gauss2*one_gauss1*r2prime*giSigma2[GaussIndex]*gRepulsive[GaussIndex]*ri12;
01249       *pairGaussEnergy = A*(-1+(one_gauss1)*(one_gauss2)*(1+gRepulsive[GaussIndex]*ri12/A));
01250     }
01251     //std::cout << "Net force: " << (LJ + Gauss) << " with ri " << (LJ + Gauss)*ri << "\n";
01252     return (LJ + Gauss)*ri;
01253   }
01254   return 0;
01255 }

int Molecule::get_groupSize ( int   ) 
Improper* Molecule::get_improper ( int  inum  )  const [inline]

Definition at line 1058 of file Molecule.h.

Referenced by dumpbench().

01058 {return (&(impropers[inum]));}

int32* Molecule::get_impropers_for_atom ( int  anum  )  [inline]

Definition at line 1138 of file Molecule.h.

Referenced by dumpbench().

01139       { return impropersByAtom[anum]; }  

Lphost* Molecule::get_lphost ( int  atomid  )  const [inline]

Definition at line 1068 of file Molecule.h.

01068                                        {
01069     // don't call unless simParams->drudeOn == TRUE
01070     // otherwise lphostIndexes array doesn't exist!
01071     int index = lphostIndexes[atomid];
01072     return (index != -1 ? &(lphosts[index]) : NULL);
01073   }

const int32* Molecule::get_mod_exclusions_for_atom ( int  anum  )  const [inline]

Definition at line 1146 of file Molecule.h.

01147       { return modExclusionsByAtom[anum]; }

int Molecule::get_mother_atom ( int   ) 
void Molecule::get_movdrag_params ( Vector v,
int  atomnum 
) const [inline]

Definition at line 1308 of file Molecule.h.

Referenced by Sequencer::addMovDragToPosition().

01309   {
01310     v = movDragParams[movDragIndexes[atomnum]].v;
01311   }

Bool Molecule::get_noPC (  )  [inline]

Definition at line 848 of file Molecule.h.

Referenced by ComputeQM::initialize().

00848 { return qmNoPC; } ;

int Molecule::get_numQMAtoms (  )  [inline]
Real* Molecule::get_qmAtmChrg (  )  [inline]
const int* Molecule::get_qmAtmIndx (  )  [inline]
Real Molecule::get_qmAtomGroup ( int  indx  )  const [inline]

Definition at line 812 of file Molecule.h.

00812 {return qmAtomGroup[indx]; } ;

const Real* Molecule::get_qmAtomGroup (  )  const [inline]
Bool Molecule::get_qmcSMD (  )  [inline]

Definition at line 859 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00859 { return qmcSMD;} ;

int* Molecule::get_qmCustomPCIdxs (  )  [inline]

Definition at line 857 of file Molecule.h.

Referenced by ComputeQM::initialize().

00857 { return qmCustomPCIdxs; } ;

int* Molecule::get_qmCustPCSizes (  )  [inline]

Definition at line 856 of file Molecule.h.

Referenced by ComputeQM::initialize().

00856 { return qmCustPCSizes; } ;

const BigReal* Molecule::get_qmDummyBondVal (  )  [inline]

Definition at line 826 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00826 { return qmDummyBondVal; } ;

const char* const* Molecule::get_qmDummyElement (  )  [inline]

Definition at line 831 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00831 { return qmDummyElement; } ;

const char* const* Molecule::get_qmElements (  )  [inline]

Definition at line 818 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00818 {return qmElementArray;} ;

const int* const* Molecule::get_qmGrpBonds (  )  [inline]

Definition at line 829 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00829 { return qmGrpBonds; } ;

Real* Molecule::get_qmGrpChrg (  )  [inline]

Definition at line 822 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00822 {return qmGrpChrg; } ;

Real* Molecule::get_qmGrpID (  )  [inline]

Definition at line 821 of file Molecule.h.

Referenced by ComputeQM::initialize(), and ComputeQMMgr::recvPartQM().

00821 { return qmGrpID; } ;

Real* Molecule::get_qmGrpMult (  )  [inline]

Definition at line 823 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00823 {return qmGrpMult; } ;

const int* Molecule::get_qmGrpNumBonds (  )  [inline]

Definition at line 827 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00827 { return qmGrpNumBonds; } ;

const int* Molecule::get_qmGrpSizes (  )  [inline]

Definition at line 820 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00820 {return qmGrpSizes; } ;

int Molecule::get_qmLSSFreq (  )  [inline]

Definition at line 839 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00839 { return qmLSSFreq; } ;

int* Molecule::get_qmLSSIdxs (  )  [inline]

Definition at line 837 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00837 { return qmLSSIdxs;} ;

Mass* Molecule::get_qmLSSMass (  )  [inline]

Definition at line 838 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00838 { return qmLSSMass; } ;

int* Molecule::get_qmLSSRefIDs (  )  [inline]

Definition at line 841 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00841 { return qmLSSRefIDs ; } ;

Mass* Molecule::get_qmLSSRefMass (  )  [inline]

Definition at line 843 of file Molecule.h.

00843 {return qmLSSRefMass; } ;

int* Molecule::get_qmLSSRefSize (  )  [inline]

Definition at line 842 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00842 { return qmLSSRefSize ; } ;

int Molecule::get_qmLSSResSize (  )  [inline]

Definition at line 840 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00840 { return qmLSSResidueSize; } ;

int* Molecule::get_qmLSSSize (  )  [inline]

Definition at line 836 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00836 { return qmLSSSize;} ;

int* Molecule::get_qmMeMMindx (  )  [inline]

Definition at line 850 of file Molecule.h.

Referenced by ComputeQM::initialize().

00850 { return qmMeMMindx; } ;

int Molecule::get_qmMeNumBonds (  )  [inline]

Definition at line 849 of file Molecule.h.

Referenced by ComputeQM::initialize(), and ComputeQMMgr::recvPartQM().

00849 { return qmMeNumBonds; } ;

Real* Molecule::get_qmMeQMGrp (  )  [inline]

Definition at line 851 of file Molecule.h.

Referenced by ComputeQM::initialize().

00851 { return qmMeQMGrp; } ;

const int* const* Molecule::get_qmMMBond (  )  [inline]

Definition at line 828 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00828 { return qmMMBond; } ;

const int* const* Molecule::get_qmMMBondedIndx (  )  [inline]

Definition at line 830 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00830 { return qmMMBondedIndx; } ;

const int* const* Molecule::get_qmMMChargeTarget (  )  [inline]

Definition at line 833 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00833 { return qmMMChargeTarget; } ;

const int* Molecule::get_qmMMNumTargs (  )  [inline]

Definition at line 834 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00834 { return qmMMNumTargs; } ;

std::map<int,int>& Molecule::get_qmMMSolv (  )  [inline]

Definition at line 844 of file Molecule.h.

00844 { return qmClassicSolv;} ;

int Molecule::get_qmNumBonds (  )  [inline]

Definition at line 825 of file Molecule.h.

00825 { return qmNumBonds; } ;

int Molecule::get_qmNumGrps (  )  const [inline]

Definition at line 819 of file Molecule.h.

Referenced by ComputeQM::initialize(), ComputeQMMgr::recvPartQM(), and ComputeQM::saveResults().

00819 {return qmNumGrps; } ;

int Molecule::get_qmPCFreq (  )  [inline]

Definition at line 853 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00853 { return qmPCFreq; } ;

Bool Molecule::get_qmReplaceAll (  )  [inline]

Definition at line 846 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00846 {return qmReplaceAll; } ;

int Molecule::get_qmTotCustPCs (  )  [inline]

Definition at line 855 of file Molecule.h.

Referenced by ComputeQM::initialize().

00855 { return qmTotCustPCs; } ;

int Molecule::get_residue_size ( const char *  segid,
int  resid 
) const

Definition at line 144 of file Molecule.C.

References ResidueLookupElem::lookup(), and NAMD_die().

Referenced by GlobalMasterFreeEnergy::getNumAtoms(), GlobalMasterEasy::getNumAtoms(), and prepare_qm().

00145                                             {
00146 
00147   if (atomNames == NULL || resLookup == NULL)
00148   {
00149     NAMD_die("Tried to find atom from name on node other than node 0");
00150   }
00151   int i = 0;
00152   int end = 0;
00153   if ( resLookup->lookup(segid,resid,&i,&end) ) return 0;
00154   return ( end - i );
00155 }

void Molecule::get_rotdrag_params ( BigReal v,
Vector a,
Vector p,
int  atomnum 
) const [inline]

Definition at line 1314 of file Molecule.h.

Referenced by Sequencer::addRotDragToPosition().

01316   {
01317     v = rotDragParams[rotDragIndexes[atomnum]].v;
01318     a = rotDragParams[rotDragIndexes[atomnum]].a;
01319     p = rotDragParams[rotDragIndexes[atomnum]].p;
01320   }

unsigned char Molecule::get_ss_type ( int  anum  )  const [inline]

Definition at line 1338 of file Molecule.h.

Referenced by ComputeHomeTuples< TholeElem, Thole, TholeValue >::loadTuples(), and Sequencer::rescaleSoluteCharges().

01339         {
01340                 return(ssAtomFlags[anum]);
01341         }

void Molecule::get_stir_refPos ( Vector refPos,
int  atomnum 
) const [inline]

Definition at line 1289 of file Molecule.h.

01290   {
01291     refPos = stirParams[stirIndexes[atomnum]].refPos;
01292   }

Real Molecule::get_stir_startTheta ( int  atomnum  )  const [inline]

Definition at line 1301 of file Molecule.h.

Referenced by ComputeStir::doForce().

01302   {
01303     return stirParams[stirIndexes[atomnum]].startTheta;
01304   }

Bond* Molecule::getAllAcceptors (  )  const [inline]

Definition at line 1095 of file Molecule.h.

01095 {return acceptors;}

Angle* Molecule::getAllAngles (  )  const [inline]

Definition at line 1078 of file Molecule.h.

Referenced by buildAngleData().

01078 {return angles;}

Bond* Molecule::getAllBonds (  )  const [inline]

Definition at line 1077 of file Molecule.h.

Referenced by buildBondData().

01077 {return bonds;}

Crossterm* Molecule::getAllCrossterms (  )  const [inline]

Definition at line 1081 of file Molecule.h.

Referenced by buildCrosstermData().

01081 {return crossterms;}

Dihedral* Molecule::getAllDihedrals (  )  const [inline]

Definition at line 1080 of file Molecule.h.

Referenced by buildDihedralData().

01080 {return dihedrals;}

Bond* Molecule::getAllDonors (  )  const [inline]

Definition at line 1094 of file Molecule.h.

01094 {return donors;}

Improper* Molecule::getAllImpropers (  )  const [inline]

Definition at line 1079 of file Molecule.h.

Referenced by buildImproperData().

01079 {return impropers;}

Lphost* Molecule::getAllLphosts (  )  const [inline]

Definition at line 1085 of file Molecule.h.

01085 { return lphosts; }

BigReal Molecule::GetAtomAlpha ( int  i  )  const [inline]

Definition at line 480 of file Molecule.h.

References drude_constants::alpha.

00480                                     {
00481     return drudeConsts[i].alpha;
00482   }

AtomNameInfo* Molecule::getAtomNames (  )  const [inline]

Definition at line 489 of file Molecule.h.

Referenced by buildAtomData().

00489 { return atomNames; }

Atom* Molecule::getAtoms (  )  const [inline]

Definition at line 488 of file Molecule.h.

References atoms.

Referenced by buildAtomData(), WorkDistrib::createAtomLists(), and outputCompressedFile().

00488 { return atoms; }

AtomSegResInfo* Molecule::getAtomSegResInfo (  )  const [inline]

Definition at line 492 of file Molecule.h.

Referenced by buildAtomData().

00492 { return atomSegResids; }

const float* Molecule::getBFactorData (  )  [inline]

Definition at line 1019 of file Molecule.h.

Referenced by NamdState::loadStructure(), and outputCompressedFile().

01019 { return (const float *)bfactor; }

BigReal Molecule::getEnergyTailCorr ( const   BigReal,
const   int 
)
int const* Molecule::getLcpoParamType (  )  [inline]

Definition at line 476 of file Molecule.h.

Referenced by HomePatch::setLcpoType().

00476                                  {
00477     return lcpoParamType;
00478   }

const float* Molecule::getOccupancyData (  )  [inline]

Definition at line 1015 of file Molecule.h.

Referenced by NamdState::loadStructure(), and outputCompressedFile().

01015 { return (const float *)occupancy; }

BigReal Molecule::getVirialTailCorr ( const   BigReal  ) 
Bool Molecule::go_restricted ( int  chain1,
int  chain2,
int  rDiff 
)

Definition at line 525 of file GoMolecule.C.

References FALSE, go_array, MAX_GO_CHAINS, MAX_RESTRICTIONS, and TRUE.

Referenced by build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), get_go_force2(), and get_go_force_new().

00526 {
00527   int i;      //  Loop counter
00528   for(i=0; i<MAX_RESTRICTIONS; i++) {
00529     if (go_array[(MAX_GO_CHAINS*chain1) + chain2].restrictions[i]  == rDiff) {
00530       return TRUE;
00531     } else if (go_array[(MAX_GO_CHAINS*chain1) + chain2].restrictions[i] == -1) {
00532       return FALSE;
00533     }
00534   }
00535   return FALSE;
00536 }

void Molecule::goInit (  ) 

Definition at line 54 of file GoMolecule.C.

References atomChainTypes, energyNative, energyNonnative, goCoordinates, goPDB, goResids, goSigmaIndices, goSigmas, goWithinCutoff, and numGoAtoms.

00054                       {
00055   numGoAtoms=0;
00056   energyNative=0;
00057   energyNonnative=0;
00058   atomChainTypes=NULL;
00059   goSigmaIndices=NULL;
00060   goSigmas=NULL;
00061   goWithinCutoff=NULL;
00062   goCoordinates=NULL;
00063   goResids=NULL;
00064   goPDB=NULL;
00065 }

void Molecule::initialize (  ) 

Referenced by Molecule().

Bool Molecule::is_atom_constorqued ( int  atomnum  )  const [inline]

Definition at line 1233 of file Molecule.h.

References consTorqueIndexes, FALSE, and numConsTorque.

01234   {
01235     if (numConsTorque)
01236     {
01237       //  Check the index to see if it is constrained
01238       return(consTorqueIndexes[atomnum] != -1);
01239     }
01240     else
01241     {
01242       //  No constraints at all, so just return FALSE
01243       return(FALSE);
01244     }
01245   }

Bool Molecule::is_atom_constrained ( int  atomnum  )  const [inline]

Definition at line 1184 of file Molecule.h.

References FALSE, and numConstraints.

Referenced by ComputeRestraints::doForce().

01185   {
01186     if (numConstraints)
01187     {
01188       //  Check the index to see if it is constrained
01189       return(consIndexes[atomnum] != -1);
01190     }
01191     else
01192     {
01193       //  No constraints at all, so just return FALSE
01194       return(FALSE);
01195     }
01196   }

Bool Molecule::is_atom_exPressure ( int  atomnum  )  const [inline]

Definition at line 1434 of file Molecule.h.

References numExPressureAtoms.

Referenced by Sequencer::langevinPiston().

01435   {
01436     return (numExPressureAtoms && exPressureAtomFlags[atomnum]);
01437   }

Bool Molecule::is_atom_fixed ( int  atomnum  )  const [inline]

Definition at line 1394 of file Molecule.h.

References numFixedAtoms.

Referenced by WorkDistrib::createAtomLists().

01395   {
01396     return (numFixedAtoms && fixedAtomFlags[atomnum]);
01397   }

Bool Molecule::is_atom_gridforced ( int  atomnum,
int  gridnum 
) const [inline]

Definition at line 1168 of file Molecule.h.

References FALSE, and numGridforceGrids.

Referenced by ComputeGridForce::do_calc().

01169   {
01170       if (numGridforceGrids)
01171       {
01172           return(gridfrcIndexes[gridnum][atomnum] != -1);
01173       }
01174       else
01175       {
01176           return(FALSE);
01177       }
01178   }

Bool Molecule::is_atom_movdragged ( int  atomnum  )  const [inline]

Definition at line 1201 of file Molecule.h.

References FALSE, and numMovDrag.

Referenced by Sequencer::addMovDragToPosition().

01202   {
01203     if (numMovDrag)
01204     {
01205       //  Check the index to see if it is constrained
01206       return(movDragIndexes[atomnum] != -1);
01207     }
01208     else
01209     {
01210       //  No constraints at all, so just return FALSE
01211       return(FALSE);
01212     }
01213   }

Bool Molecule::is_atom_rotdragged ( int  atomnum  )  const [inline]

Definition at line 1217 of file Molecule.h.

References FALSE, and numRotDrag.

Referenced by Sequencer::addRotDragToPosition().

01218   {
01219     if (numRotDrag)
01220     {
01221       //  Check the index to see if it is constrained
01222       return(rotDragIndexes[atomnum] != -1);
01223     }
01224     else
01225     {
01226       //  No constraints at all, so just return FALSE
01227       return(FALSE);
01228     }
01229   }

Bool Molecule::is_atom_stirred ( int  atomnum  )  const [inline]

Definition at line 1415 of file Molecule.h.

References FALSE, and numStirredAtoms.

Referenced by ComputeStir::doForce().

01416   {
01417     if (numStirredAtoms)
01418     {
01419       //  Check the index to see if it is constrained
01420       return(stirIndexes[atomnum] != -1);
01421     }
01422     else
01423     {
01424       //  No constraints at all, so just return FALSE
01425       return(FALSE);
01426     }
01427   }

Bool Molecule::is_drude ( int   ) 
Bool Molecule::is_group_fixed ( int  atomnum  )  const [inline]

Definition at line 1430 of file Molecule.h.

References numFixedAtoms.

01431   {
01432     return (numFixedAtoms && (fixedAtomFlags[atomnum] == -1));
01433   }

Bool Molecule::is_hydrogen ( int   ) 
Bool Molecule::is_hydrogenGroupParent ( int   ) 
Bool Molecule::is_lp ( int   ) 
Bool Molecule::is_oxygen ( int   ) 
Bool Molecule::is_water ( int   ) 
Real Molecule::langevin_param ( int  atomnum  )  const [inline]

Definition at line 1283 of file Molecule.h.

Referenced by WorkDistrib::createAtomLists().

01284   {
01285     return(langevinParams ? langevinParams[atomnum] : 0.);
01286   }

int64_t Molecule::num_deg_freedom ( int  isInitialReport = 0  )  const [inline]

Definition at line 521 of file Molecule.h.

References num_fixed_atoms(), numAtoms, numConstraints, numFepInitial, numFixedRigidBonds, numLonepairs, numRigidBonds, simParams, and WAT_TIP4.

Referenced by ParallelIOMgr::bcastHydroBasedCounter(), Controller::Controller(), NamdState::loadStructure(), and Controller::receivePressure().

00521                                                          {
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   }

int Molecule::num_fixed_atoms (  )  const [inline]

Definition at line 495 of file Molecule.h.

References numFixedAtoms, and simParams.

Referenced by num_deg_freedom(), num_fixed_groups(), num_group_deg_freedom(), and Controller::receivePressure().

00495                               {
00496     // local variables prefixed by s_
00497     int s_NumFixedAtoms = (simParams->fixedAtomsOn ? numFixedAtoms : 0);
00498     return s_NumFixedAtoms;  // value is "turned on" SimParameters
00499   }

int Molecule::num_fixed_groups (  )  const [inline]

Definition at line 501 of file Molecule.h.

References num_fixed_atoms(), and numFixedGroups.

Referenced by num_group_deg_freedom(), and Controller::receivePressure().

00501                                {
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   }

int64_t Molecule::num_group_deg_freedom (  )  const [inline]

Definition at line 508 of file Molecule.h.

References num_fixed_atoms(), num_fixed_groups(), numConstraints, numHydrogenGroups, and simParams.

Referenced by Controller::receivePressure().

00508                                         {
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   }

void Molecule::prepare_qm ( const char *  pdbFileName,
Parameters params,
ConfigList cfgList 
)

Prepares Live Solvent Selection

Data gathering from PDB to find QM atom and bond info

Multiplicity of each QM group

Charge of each QM group

Populate arrays that are used throughout the the calculations.

Overides Link Atom element with user selection.

Bond Schemes. Prepares for treatment of QM-MM bonds in ComputeQM.C

Live Solvent Selection

Custom Point Charge selection

Topology preparation

Definition at line 109 of file MoleculeQM.C.

References ResizeArray< Elem >::add(), Parameters::assign_vdw_index(), qmSolvData::atmIDs, PDB::atom(), bond::atom1, bond::atom2, ResizeArray< Elem >::begin(), bond::bond_type, atom_constants::charge, StringList::data, DebugM, PDBAtom::element(), ResizeArray< Elem >::end(), endi(), SimParameters::extraBondsOn, SortedArray< Elem >::find(), ConfigList::find(), SimParameters::fixedAtomsOn, get_residue_size(), ObjectArena< Type >::getNewArray(), iERROR(), iINFO(), SortedArray< Elem >::insert(), iout, iWARN(), ResidueLookupElem::lookup(), atom_constants::mass, NAMD_blank_string(), NAMD_die(), NAMD_read_line(), StringList::next, PDB::num_atoms(), numAtoms, numBonds, numRealBonds, PDBAtom::occupancy(), SimParameters::PMEOn, SimParameters::qmBondDist, SimParameters::qmBondOn, SimParameters::qmBondScheme, SimParameters::qmChrgFromPSF, SimParameters::qmColumn, SimParameters::qmCSMD, SimParameters::qmCustomPCSel, SimParameters::qmElecEmbed, SimParameters::qmFormat, QMFormatMOPAC, QMFormatORCA, SimParameters::qmLSSFreq, SimParameters::qmLSSMode, QMLSSMODECOM, QMLSSMODEDIST, SimParameters::qmLSSOn, SimParameters::qmLSSResname, SimParameters::qmMOPACAddConfigChrg, SimParameters::qmNoPC, SimParameters::qmPCSelFreq, QMSCHEMECS, QMSCHEMERCD, QMSCHEMEZ1, QMSCHEMEZ2, QMSCHEMEZ3, SimParameters::qmVDW, PDBAtom::residuename(), PDBAtom::residueseq(), PDBAtom::segmentname(), ResizeArray< Elem >::size(), split(), SimParameters::stepsPerCycle, and PDBAtom::temperaturefactor().

Referenced by NamdState::loadStructure().

00110                                                                         {
00111 #ifdef MEM_OPT_VERSION
00112     NAMD_die("QMMM interface is not supported in memory optimized builds.");
00113 #else
00114     
00115     PDB *pdbP;      //  Pointer to PDB object to use
00116     
00117     qmNumQMAtoms = 0;
00118     
00119     qmNumGrps = 0 ;
00120     
00121     iout << iINFO << "Using the following PDB file for QM parameters: " << 
00122     pdbFileName << "\n" << endi;
00123     if (pdbFileName)
00124         pdbP = new PDB(pdbFileName);
00125     else
00126         iout << "pdbFile name not defined!\n\n" << endi ;
00127     
00128     if (pdbP->num_atoms() != numAtoms)
00129     {
00130        NAMD_die("Number of atoms in QM parameter PDB file doesn't match coordinate PDB");
00131     }
00132     
00133     qmElementArray = new char*[numAtoms] ;
00134     
00135     qmAtomGroup = new Real[numAtoms] ;
00136     
00137     BigReal bondVal;
00138     Real atmGrp;
00139     
00140     std::set<Real> atmGrpSet;
00141     
00142     std::vector<Real> grpChrgVec;
00143     
00144     // Value in the bond column fro QM-MM bonds.
00145     std::vector<BigReal> qmBondValVec;
00146     // Identifiers of different QM regions
00147     std::vector<Real> qmGroupIDsVec;
00148     // This maps the qm group ID with the serail index of this group in 
00149     // various arrays.
00150     std::map<Real,int> qmGrpIDMap ;
00151     
00152     std::vector<int> qmAtmIndxVec, qmGrpSizeVec;
00153     
00154     // Used to parse user selection in different places.
00155     std::vector<std::string> strVec;
00156     
00157     qmNoPC = simParams->qmNoPC;
00158     
00159     // We set to zero by default, So there is no need for extra processing.
00160     qmPCFreq = 0;
00161     if (simParams->qmPCSelFreq > 1)
00162         qmPCFreq = simParams->qmPCSelFreq;
00163     
00164     
00167     
00168     
00169     qmLSSTotalNumAtms = 0;
00170     SortedArray< qmSolvData> lssGrpRes;
00171     std::map<Real,std::vector<int> > lssGrpRefIDs;
00172     refSelStrMap lssRefUsrSel;
00173     int totalNumRefAtms = 0;
00174     int lssClassicResIndx = 0 ;
00175     int lssCurrClassResID = -1 ;
00176     char lssCurrClassResSegID[5];
00177     if (simParams->qmLSSOn) {
00178         DebugM(4, "LSS is ON! Processing QM solvent.\n") ;
00179         
00180         if (resLookup == NULL)
00181             NAMD_die("We need residue data to conduct LSS.");
00182          
00183         if (simParams->qmLSSMode == QMLSSMODECOM) {
00184             
00185             StringList *current = cfgList->find("QMLSSRef");
00186             for ( ; current; current = current->next ) {
00187                 
00188                 strVec = split( std::string(current->data) , " ");
00189                 
00190                 if (strVec.size() != 3 ) {
00191                     iout << iERROR << "Format error in QM LSS size: " 
00192                     << current->data
00193                     << "\n" << endi;
00194                     NAMD_die("Error processing QM information.");
00195                 }
00196                 
00197                 std::stringstream storConv ;
00198                 
00199                 storConv << strVec[0] ;
00200                 Real grpID ;
00201                 storConv >> grpID;
00202                 if (storConv.fail()) {
00203                     iout << iERROR << "Error parsing QMLSSRef selection: " 
00204                     << current->data
00205                     << "\n" << endi;
00206                     NAMD_die("Error processing QM information.");
00207                 }
00208                 
00209                 std::string segName = strVec[1].substr(0,4);
00210                 
00211                 storConv.clear() ;
00212                 storConv << strVec[2];
00213                 int resID ;
00214                 storConv >> resID;
00215                 if (storConv.fail()) {
00216                     iout << iERROR << "Error parsing QMLSSRef selection: " 
00217                     << current->data
00218                     << "\n" << endi;
00219                     NAMD_die("Error processing QM information.");
00220                 }
00221                 
00222                 auto it = lssRefUsrSel.find(grpID) ;
00223                 if (it == lssRefUsrSel.end())
00224                     lssRefUsrSel.insert(refSelStrPair(grpID,refSelStrVec()));
00225                 
00226                 lssRefUsrSel[grpID].push_back(refSelStr(segName, resID));
00227             }
00228             
00229             for (auto it=lssRefUsrSel.begin(); it!=lssRefUsrSel.end(); it++) {
00230                 iout << iINFO << "LSS user selection for COM composition of group "
00231                 << it->first << ":\n" << endi ;
00232                 for (int i=0; i<it->second.size();i++) {
00233                     iout << iINFO << "Segment " << it->second[i].segid 
00234                     << " ; residue " << it->second[i].resid << "\n" << endi ;
00235                 }
00236             }
00237         }
00238     }
00239     
00240     
00241     
00244     
00245     
00246     for (int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
00247         
00248         // If we are looking for QM-MM bonds, then we need to store extra info.
00249         if (simParams->qmBondOn) {
00250             
00251             // We store both the qm group and the bond value
00252             if ( strcmp(simParams->qmColumn,"beta") == 0 ){
00253                 atmGrp = pdbP->atom(atmInd)->temperaturefactor() ;
00254                 
00255                 bondVal = pdbP->atom(atmInd)->occupancy() ;
00256             }
00257             else {
00258                 atmGrp = pdbP->atom(atmInd)->occupancy() ;
00259                 
00260                 bondVal = pdbP->atom(atmInd)->temperaturefactor() ;
00261             }
00262             
00263             qmBondValVec.push_back(bondVal);
00264         }
00265         else {
00266             
00267             if ( strcmp(simParams->qmColumn,"beta") == 0 ){
00268                 atmGrp = pdbP->atom(atmInd)->temperaturefactor() ;
00269             }
00270             else {
00271                 atmGrp = pdbP->atom(atmInd)->occupancy() ;
00272             }
00273         }
00274         
00275         // We store all atom QM Group IDs in an array for 
00276         // future transport to all PEs.
00277         qmAtomGroup[atmInd] = atmGrp;
00278         
00279         // We store all atom's elements for quick access in the QM code.
00280         // Elements are used to tell the QM code the atom type.
00281         qmElementArray[atmInd] = new char[3];
00282         strncpy(qmElementArray[atmInd],pdbP->atom(atmInd)->element(),3);
00283         
00284         // For QM atoms, we keep global atom indices and parse different QM Group
00285         // IDs, keeping track of each groups size
00286         if (atmGrp > 0){
00287             
00288             if (simParams->fixedAtomsOn) {
00289                 if ( fixedAtomFlags[atmInd] == 1 ) {
00290                     iout << iERROR << "QM atom cannot be fixed in space!\n" 
00291                     << endi;
00292                     NAMD_die("Error processing QM information.");
00293                 }
00294             }
00295             
00296             // Changes the VdW type of QM atoms.
00297             // This is may be used to counter balance the overpolarization 
00298             // that QM atoms sufer.
00299             if (simParams->qmVDW) {
00300                 // Creating a new type string
00301                 std::string qmType(qmElementArray[atmInd]) ;
00302                 // Erases empty spaces
00303                 qmType.erase(
00304                         std::remove_if(qmType.begin(), qmType.end(), isspace ),
00305                     qmType.end());
00306                 // pre-pends a "q" to the element name.
00307                 qmType = std::string("q") + qmType;
00308             
00309 //                 iout << "QM VdW type: " << atoms[atmInd].vdw_type 
00310 //                 << " atom type: " << atomNames[atmInd].atomtype 
00311 //                 << " new type "<< qmType << "\n" << endi;
00312                 
00313                 /*  Look up the vdw constants for this atom    */
00314                 // I am casting a non-const here because the function doesn't actually
00315                 // change the string value, but it doesn't ask for a const char* as 
00316                 // an argument.
00317                 params->assign_vdw_index(const_cast<char*>( qmType.c_str()), 
00318                                          &(atoms[atmInd]));
00319                 
00320 //                 iout << "----> new VdW type: " << atoms[atmInd].vdw_type << "\n" << endi;
00321             }
00322             
00323             // Indexes all global indices of QM atoms.
00324             qmAtmIndxVec.push_back(atmInd);
00325             
00326             auto retVal = atmGrpSet.insert(atmGrp);
00327             
00328             // If a new QM group is found, store the reverse reference in a map
00329             // and increase the total count.
00330             if (retVal.second) {
00331                 // This mak makes the reverse identification from the group ID
00332                 // to the sequential order in which each group was first found.
00333                 qmGrpIDMap.insert(std::pair<BigReal,int>(atmGrp, qmNumGrps));
00334                 
00335                 // This vector keeps track of the group ID for each group
00336                 qmGroupIDsVec.push_back(atmGrp);
00337                 
00338                 // This counter is used to keep track of the sequential order in
00339                 // which QM groups were first seen in the reference PDB file.
00340                 qmNumGrps++ ;
00341                 
00342                 // If this is a new QM group, initialize a new variable in the 
00343                 // vector to keep track of the number of atoms in this group.
00344                 qmGrpSizeVec.push_back(1);
00345                 
00346                 // For each new QM group, a new entry in the total charge
00347                 // vector is created
00348                 grpChrgVec.push_back( atoms[atmInd].charge );
00349             }
00350             else {
00351                 // If the QM group has already been seen, increase the count
00352                 // of atoms in that group.
00353                 qmGrpSizeVec[ qmGrpIDMap[atmGrp] ]++ ;
00354                 
00355                 // If the QM group exists, adds current atom charge to total charge.
00356                 grpChrgVec[ qmGrpIDMap[atmGrp] ] += atoms[atmInd].charge;
00357             }
00358             
00359             // In case we are using live solvent selection
00360             if(simParams->qmLSSOn) {
00361                 
00362                 int resID = pdbP->atom(atmInd)->residueseq();
00363                 char segName[5];
00364                 strncpy(segName, pdbP->atom(atmInd)->segmentname(),5);
00365                 
00366                 int resSize = get_residue_size(segName,resID);
00367                 
00368                 int i =-1, end =-1;
00369                 
00370                 resLookup->lookup(segName,resID,&i, &end);
00371                 
00372                 if ( strncasecmp(pdbP->atom(atmInd)->residuename(),
00373                            simParams->qmLSSResname, 4) == 0) {
00374                     
00375                     // We try to insert every residue from every atom
00376                     qmSolvData *resP = lssGrpRes.find(qmSolvData(resID, i, resSize));
00377                     
00378                     if (resP != NULL) {
00379                         resP->atmIDs.push_back(atmInd) ;
00380 //                         DebugM(3, "Existing residue " << resP->resID 
00381 //                         << " from segID " << resP->segName
00382 //                         << " received atom "
00383 //                         << atmInd << "\n" );
00384                     }
00385                     else {
00386                         lssGrpRes.insert(qmSolvData(resID, i, resSize, segName, atmGrp));
00387 //                         DebugM(3, lssGrpRes.size() << ") Inserted new residue: " 
00388 //                         << resID << " from segID " << segName
00389 //                         << " with atom "
00390 //                         << i << "\n" ) ;
00391                     }
00392                 }
00393                 else {
00394                     // We store the QM atoms, per group, which are NOT part of
00395                     // solvent molecules.
00396                     
00397                     // Checks if we have a vector for this QM group.
00398                     auto itNewGrp = lssGrpRefIDs.find(atmGrp) ;
00399                     if (itNewGrp == lssGrpRefIDs.end()) {
00400                         lssGrpRefIDs.insert(std::pair<Real, std::vector<int> >(
00401                             atmGrp, std::vector<int>() ) );
00402                     }
00403                     
00404                     switch (simParams->qmLSSMode)
00405                     {
00406                     
00407                     case QMLSSMODECOM:
00408                     {
00409                         // If we are using COM selection, checks if the atom
00410                         // is part of the user selected residues
00411                         for(int i=0; i<lssRefUsrSel[atmGrp].size(); i++) {
00412                             if (lssRefUsrSel[atmGrp][i].resid == resID &&
00413                                 strncmp(lssRefUsrSel[atmGrp][i].segid.c_str(),
00414                                         segName,5) == 0 ) {
00415                                 
00416                                 lssGrpRefIDs[atmGrp].push_back(atmInd);
00417                                 totalNumRefAtms++;
00418                             }
00419                         }
00420                         
00421                     } break;
00422                     
00423                     case QMLSSMODEDIST:
00424                     {
00425                     
00426                         lssGrpRefIDs[atmGrp].push_back(atmInd);
00427                         totalNumRefAtms++;
00428                     
00429                     } break;
00430                     
00431                     }
00432                 }
00433                 
00434             }
00435             
00436         }
00437         else if (atmGrp == 0) {
00438             
00439             if(simParams->qmLSSOn) {
00440                 
00441                 if ( strncasecmp(pdbP->atom(atmInd)->residuename(),
00442                            simParams->qmLSSResname, 4) == 0) {
00443                     
00444                     int resID = pdbP->atom(atmInd)->residueseq();
00445                     char segName[5];
00446                     strncpy(segName, pdbP->atom(atmInd)->segmentname(),5);
00447                     
00448                     if (lssCurrClassResID < 0) {
00449                         lssCurrClassResID = resID ;
00450                         strncpy(lssCurrClassResSegID, segName,5);
00451                         lssClassicResIndx = 0;
00452                     }
00453                     else if (lssCurrClassResID != resID ||
00454                             strcmp(lssCurrClassResSegID, segName) != 0 ) {
00455                         lssCurrClassResID = resID ;
00456                         strncpy(lssCurrClassResSegID, segName,5);
00457                         lssClassicResIndx++;
00458                     }
00459                     
00460                     qmClassicSolv.insert(std::pair<int,int>(atmInd,lssClassicResIndx));
00461                     
00462                 }
00463             }
00464         }
00465         else if(atmGrp < 0) {
00466             iout << iERROR << "QM group ID cannot be less than zero!\n" << endi;
00467             NAMD_die("Error processing QM information.");
00468         }
00469     }
00470     
00471     // Sanity check
00472     if (simParams->qmLSSOn) {
00473         if (lssGrpRes.size() == 0)
00474             NAMD_die("LSS was activated but there are no QM solvent molecules!\n");
00475         
00476         for (auto it=lssRefUsrSel.begin(); it!=lssRefUsrSel.end(); it++) {
00477             auto itTarget = qmGrpIDMap.find(it->first);
00478             if (itTarget == qmGrpIDMap.end()) {
00479                 iout << iERROR << "QM group ID for LSS could not be found in input!"
00480                 << " QM group ID: " << it->first << "\n" << endi;
00481                 NAMD_die("Error processing QM information.");
00482             }
00483         }
00484         
00485         DebugM(3,"We found " << lssClassicResIndx << " classical solvent molecules totalling "
00486         << qmClassicSolv.size() << " atoms.\n" );
00487         
00488     }
00489     
00490     qmNumQMAtoms = qmAtmIndxVec.size();
00491     
00492     if (qmNumQMAtoms == 0)
00493         NAMD_die("No QM atoms were found in this QM simulation!") ;
00494     
00495     iout << iINFO << "Number of QM atoms (excluding Dummy atoms): " << 
00496     qmNumQMAtoms << "\n" << endi;
00497     
00498     qmAtmChrg = new Real[qmNumQMAtoms] ;
00499     qmAtmIndx = new int[qmNumQMAtoms] ;
00500     for (int i=0; i<qmNumQMAtoms; i++) {
00501         // qmAtmChrg gets initialized with the PSF charges at the end of this
00502         // function, but values may change as QM steps are taken.
00503         qmAtmChrg[i] = 0;  
00504         qmAtmIndx[i] = qmAtmIndxVec[i] ;
00505     }
00506     
00507     // This map relates the QM group index with a vector of pairs
00508     // of bonded MM-QM atoms (with the bonded atoms ins this order: 
00509     // MM first, QM second).
00510     std::map<int, std::vector<std::pair<int,int> > > qmGrpIDBonds ;
00511     int bondCounter = 0;
00512     if (simParams->qmBondOn) {
00513         
00514         // Checks all bonds for QM-MM pairs.
00515         // Makes sanity checks against QM-QM pairs and MM-MM pairs that
00516         // are flagged by the user to be bonds between QM and MM regions.
00517         // QM-QM bonds will be removed in another function.
00518         for (int bndIt = 0; bndIt < numBonds; bndIt++) {
00519             
00520             bond curr = bonds[bndIt] ;
00521             
00522             // In case either atom has a non-zero
00523             if ( qmBondValVec[curr.atom1] != 0 &&
00524                  qmBondValVec[curr.atom2] != 0 ) {
00525                 
00526                 // Sanity checks (1 of 2)
00527                 if (qmAtomGroup[curr.atom1] != 0 &&
00528                     qmAtomGroup[curr.atom2] != 0) {
00529                     iout << iERROR << "Atoms " << curr.atom1 << " and " << 
00530                     curr.atom2 << " are assigned as QM atoms.\n" << endi;
00531                     NAMD_die("Error in QM-MM bond assignment.") ;
00532                 }
00533                 
00534                 // Sanity checks (2 of 2)
00535                 if (qmAtomGroup[curr.atom1] == 0 &&
00536                     qmAtomGroup[curr.atom2] == 0) {
00537                     iout << iERROR << "Atoms " << curr.atom1 << " and " << 
00538                     curr.atom2 << " are assigned as MM atoms.\n" << endi;
00539                     NAMD_die("Error in QM-MM bond assignment.") ;
00540                 }
00541                 
00542                 int currGrpID = 0;
00543                 std::pair<int,int> newPair(0,0);
00544                 
00545                 // We create a QM-MM pair with the MM atom first
00546                 if (qmAtomGroup[curr.atom1] != 0) {
00547                     newPair = std::pair<int,int>(curr.atom2, curr.atom1) ;
00548                     currGrpID = qmAtomGroup[curr.atom1];
00549                 } else {
00550                     newPair = std::pair<int,int>(curr.atom1, curr.atom2) ;
00551                     currGrpID = qmAtomGroup[curr.atom2];
00552                 } 
00553                 
00554                 int grpIndx = qmGrpIDMap[currGrpID] ;
00555                 
00556                 // Stores bonds in vectors key-ed by the QM group ID of the QM atom.
00557                 auto retIter = qmGrpIDBonds.find(grpIndx) ;
00558                 // In case thi QM-MM bonds belong to a QM group we have not seen 
00559                 // yet, stores a new vector in the map.
00560                 if (retIter == qmGrpIDBonds.end()) {
00561                     qmGrpIDBonds.insert(std::pair<BigReal,std::vector<std::pair<int,int> > >(
00562                     grpIndx, std::vector<std::pair<int,int> >() ) );
00563                 }
00564                 
00565                 qmGrpIDBonds[grpIndx].push_back( newPair );
00566                 
00567                 bondCounter++ ;
00568             }
00569             
00570             
00571         }
00572         
00573 //         iout << "Finished processing QM-MM bonds.\n" << endi ;
00574         
00575         if (bondCounter == 0)
00576             iout << iWARN << "We found " << bondCounter << " QM-MM bonds.\n" << endi ;
00577         else
00578             iout << iINFO << "We found " << bondCounter << " QM-MM bonds.\n" << endi ;
00579         
00580     }
00581     
00582     // Initializes several arrays used to setup the QM simulation.
00583     qmNumBonds = bondCounter ;
00584     
00585     qmMMBond = new int*[qmNumBonds];
00586     
00587     qmDummyBondVal = new BigReal[qmNumBonds];
00588     
00589     qmMMChargeTarget = new int*[qmNumBonds] ;
00590     qmMMNumTargs = new int[qmNumBonds] ;
00591     
00592     qmDummyElement = new char*[qmNumBonds] ;
00593     
00594     
00595     qmNumGrps = atmGrpSet.size();
00596     
00597     qmGrpSizes = new int[qmNumGrps];
00598     
00599     qmGrpID = new Real[qmNumGrps];
00600     
00601     qmGrpChrg = new Real[qmNumGrps];
00602     
00603     qmGrpMult = new Real[qmNumGrps] ;
00604     
00605     qmGrpNumBonds = new int[qmNumGrps];
00606     
00607     qmGrpBonds = new int*[qmNumGrps];
00608     qmMMBondedIndx = new int*[qmNumGrps];
00609     
00610     
00613     
00614     
00615     // We first initialize the multiplicity vector with 
00616     // default values, then populate it with user defined values.
00617     for (int grpIter = 0; grpIter < qmNumGrps; grpIter++) {
00618         qmGrpMult[grpIter] = 0;
00619     }
00620     
00621     int multCount = 0;
00622     StringList *current = cfgList->find("QMMult");
00623     for ( ; current; current = current->next ) {
00624         
00625         auto strVec = split( std::string(current->data) , " ");
00626         
00627         if (strVec.size() != 2 ) {
00628             iout << iERROR << "Format error in QM Multiplicity string: " 
00629             << current->data
00630             << "\n" << endi;
00631             NAMD_die("Error processing QM information.");
00632         }
00633         
00634         std::stringstream storConv ;
00635         
00636         storConv << strVec[0] ;
00637         Real grpID ;
00638         storConv >> grpID;
00639         if (storConv.fail()) {
00640             iout << iERROR << "Error parsing QMMult selection: " 
00641             << current->data
00642             << "\n" << endi;
00643             NAMD_die("Error processing QM information.");
00644         }
00645         
00646         storConv.clear() ;
00647         storConv << strVec[1];
00648         Real multiplicity ;
00649         storConv >> multiplicity;
00650         if (storConv.fail()) {
00651             iout << iERROR << "Error parsing QMMult selection: " 
00652             << current->data
00653             << "\n" << endi;
00654             NAMD_die("Error processing QM information.");
00655         }
00656         
00657         auto it = qmGrpIDMap.find(grpID);
00658         
00659         if (it == qmGrpIDMap.end()) {
00660             iout << iERROR << "Error parsing QMMult selection. Could not find QM group ID: " 
00661             << grpID
00662             << "\n" << endi;
00663             NAMD_die("Error processing QM information.");
00664         }
00665         else {
00666             iout << iINFO << "Applying user defined multiplicity "
00667             << multiplicity << " to QM group ID " << grpID << "\n" << endi;
00668             qmGrpMult[it->second] = multiplicity;
00669         }
00670         
00671         multCount++;
00672     }
00673     
00674     if (multCount != qmNumGrps && simParams->qmFormat == QMFormatORCA) {
00675         NAMD_die("ORCA neds multiplicity values for all QM regions.");
00676     }
00677     
00678     
00681     
00682     
00683     for (size_t grpIndx = 0; grpIndx < qmGrpSizeVec.size(); grpIndx++) {
00684         
00685         bool nonInteger = true;
00686         if ((fabsf(roundf(grpChrgVec[grpIndx]) - grpChrgVec[grpIndx]) <= 0.001f) ) {
00687             grpChrgVec[grpIndx] = roundf(grpChrgVec[grpIndx]) ;
00688             nonInteger = false;
00689         }
00690         
00691         iout << iINFO << grpIndx + 1 << ") Group ID: " << qmGroupIDsVec[grpIndx]
00692         << " ; Group size: " << qmGrpSizeVec[grpIndx] << " atoms"
00693         << " ; Total PSF charge: " << grpChrgVec[grpIndx] << "\n" << endi ;
00694         
00695         if (nonInteger && simParams->PMEOn)
00696             NAMD_die("QM atoms do not add up to a whole charge, which is needed for PME.") ;
00697     }
00698     
00699     int chrgCount = 0;
00700     current = cfgList->find("QMCharge");
00701     for ( ; current; current = current->next ) {
00702         
00703         auto strVec = split( std::string(current->data) , " ");
00704         
00705         if (strVec.size() != 2 ) {
00706             iout << iERROR << "Format error in QM Charge string: " 
00707             << current->data
00708             << "\n" << endi;
00709             NAMD_die("Error processing QM information.");
00710         }
00711         
00712         std::stringstream storConv ;
00713         
00714         storConv << strVec[0] ;
00715         Real grpID ;
00716         storConv >> grpID;
00717         if (storConv.fail()) {
00718             iout << iERROR << "Error parsing QMCharge selection: " 
00719             << current->data
00720             << "\n" << endi;
00721             NAMD_die("Error processing QM information.");
00722         }
00723         
00724         storConv.clear() ;
00725         storConv << strVec[1];
00726         Real charge ;
00727         storConv >> charge;
00728         if (storConv.fail()) {
00729             iout << iERROR << "Error parsing QMCharge selection: " 
00730             << current->data
00731             << "\n" << endi;
00732             NAMD_die("Error processing QM information.");
00733         }
00734         
00735         auto it = qmGrpIDMap.find(grpID);
00736         
00737         if (it == qmGrpIDMap.end()) {
00738             iout << iERROR << "Error parsing QMCharge selection. Could not find QM group ID: " 
00739             << grpID
00740             << "\n" << endi;
00741             NAMD_die("Error processing QM information.");
00742         }
00743         else {
00744             iout << iINFO << "Found user defined charge "
00745             << charge << " for QM group ID " << grpID << ". Will ignore PSF charge.\n" << endi;
00746             grpChrgVec[it->second] = charge;
00747         }
00748         
00749         chrgCount++;
00750     }
00751     
00752     simParams->qmMOPACAddConfigChrg = false;
00753     // Checks if QM group charges can be defined for MOPAC.
00754     // Since charge can be defined in QMConfigLine, we need extra logic.
00755     if (simParams->qmFormat == QMFormatMOPAC) {
00756         
00757         // Checks if group charge was supplied in config line for MOPAC.
00758         std::string::size_type chrgLoc = std::string::npos ;
00759         
00760         // This will hold a sting with the first user supplied configuration line for MOPAC.
00761         std::string configLine(cfgList->find("QMConfigLine")->data) ;
00762         chrgLoc = configLine.find("CHARGE") ;
00763         
00764         if ( chrgLoc != std::string::npos ) {
00765             iout << iINFO << "Found user defined charge in command line. This \
00766 will be used for all QM regions and will take precedence over all other charge \
00767 definitions.\n" << endi;
00768         }
00769         else if ( chrgLoc == std::string::npos && (simParams->qmChrgFromPSF || chrgCount == qmNumGrps) ) {
00770         // If no charge was defined in the configuration line, gets from PSF and/or 
00771         // from user defined charges (through the QMCharge keyword).
00772             simParams->qmMOPACAddConfigChrg = true;
00773         }
00774         else
00775         {
00776         // If we could nont find a charge definition in the config line AND
00777         // no specific charge was selected for each QM region through QMCharge AND
00778         // the QMChargeFromPSF was not turned ON, the we stop NAMD and scream at the user.
00779 //         if ( chrgLoc == std::string::npos && (chrgCount != qmNumGrps ) && !simParams->qmChrgFromPSF) 
00780             iout << iERROR << "Could not find charge for all QM groups. For MOPAC, \
00781 charges can be defined through QMConfigLine, QMCharge or QMChargeFromPSF keywords.\n" << endi;
00782             NAMD_die("Error processing QM information.");
00783         }
00784         
00785     }
00786     
00787     if (simParams->qmFormat == QMFormatORCA ) {
00788         if ((chrgCount != qmNumGrps ) && !simParams->qmChrgFromPSF) {
00789             // If we are not supposed to get charges from the PSF, and not 
00790             // enough charges were set to cover all QM regions, 
00791             // we stop NAMD and scream at the user.
00792             iout << iERROR << "Could not find charge for all QM groups. For ORCA, \
00793 charges can be defined through QMCharge or QMChargeFromPSF keywords.\n" << endi;
00794             NAMD_die("Error processing QM information.");
00795         }
00796     }
00797     
00798     // If mechanichal embedding was requested but we have QM-MM bonds, we need 
00799     // to send extra info to ComputeQM to preserve calculation speed.
00800     if (qmNumBonds > 0 && qmNoPC) {
00801         qmMeNumBonds = qmNumBonds;
00802         qmMeMMindx = new int[qmMeNumBonds] ;
00803         qmMeQMGrp = new Real[qmMeNumBonds] ;
00804     }
00805     else {
00806         qmMeNumBonds = 0 ;
00807         qmMeMMindx = NULL;
00808         qmMeQMGrp = NULL;
00809     }
00810     
00811     
00814     
00815     
00816     bondCounter = 0;
00817     for (int grpIter = 0; grpIter < qmNumGrps; grpIter++) {
00818         
00819         qmGrpID[grpIter] = qmGroupIDsVec[grpIter] ;
00820         
00821         qmGrpSizes[grpIter] = qmGrpSizeVec[grpIter] ;
00822         
00823         qmGrpChrg[grpIter] = grpChrgVec[grpIter];
00824         
00825 //         iout << "Loaded " << qmGrpSizes[grpIter] << " atoms to this group.\n" << endi;
00826         
00827         int currNumbBonds = qmGrpIDBonds[grpIter].size() ;
00828         
00829         // Assigns the number of bonds that the current QM group has.
00830         qmGrpNumBonds[grpIter] = currNumbBonds;
00831         
00832         if (currNumbBonds > 0) {
00833             
00834             qmGrpBonds[grpIter] = new int[currNumbBonds];
00835             qmMMBondedIndx[grpIter] = new int[currNumbBonds];
00836             
00837             for (int bndIter = 0; bndIter<currNumbBonds; bndIter++) {
00838                 
00839                 // Adds the bonds to the overall sequential list.
00840                 qmMMBond[bondCounter] = new int[2] ;
00841                 qmMMBond[bondCounter][0] = qmGrpIDBonds[grpIter][bndIter].first ;
00842                 qmMMBond[bondCounter][1] = qmGrpIDBonds[grpIter][bndIter].second ;
00843                 
00844                 // For the current QM group, and the current bond, gets the bond index.
00845                 qmGrpBonds[grpIter][bndIter] =  bondCounter;
00846                 
00847                 // For the current QM group, and the current bond, gets the MM atom.
00848                 qmMMBondedIndx[grpIter][bndIter] = qmGrpIDBonds[grpIter][bndIter].first ;
00849                 
00850                 // Assign the default value of dummy element
00851                 qmDummyElement[bondCounter] = new char[3];
00852                 strcpy(qmDummyElement[bondCounter],"H\0");
00853                 
00854                 // Determines the distance that will separate the new Dummy atom
00855                 // and the Qm atom to which it will be bound.
00856                 bondVal = 0;
00857                 if (simParams->qmBondDist) {
00858                     if (qmBondValVec[qmMMBond[bondCounter][0]] != 
00859                         qmBondValVec[qmMMBond[bondCounter][1]]
00860                     ) {
00861                         iout << iERROR << "qmBondDist is ON but the values in the bond column are different for atom "
00862                         << qmMMBond[bondCounter][0] << " and " << qmMMBond[bondCounter][1]
00863                         << ".\n" << endi ;
00864                         NAMD_die("Error in QM-MM bond processing.");
00865                     }
00866                     
00867                     bondVal = qmBondValVec[qmMMBond[bondCounter][0]] ;
00868                 } 
00869                 else {
00870                     
00871                     if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"C") == 0 ) {
00872                         bondVal = 1.09 ;
00873                     }
00874                     else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"O") == 0 ) {
00875                         bondVal = 0.98 ;
00876                     }
00877                     else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"N") == 0 ) {
00878                         bondVal = 0.99 ;
00879                     }
00880                     else {
00881                         iout << iERROR << "qmBondDist is OFF but the QM atom has no default distance value. Atom "
00882                         << qmMMBond[bondCounter][1] << ", with element: " 
00883                         << qmElementArray[qmMMBond[bondCounter][1]] 
00884                         <<  ".\n" << endi ;
00885                         NAMD_die("Error in QM-MM bond processing.");
00886                     }
00887                     
00888                 }
00889                 
00890                 iout << iINFO << "MM-QM pair: " << qmMMBond[bondCounter][0] << ":"
00891                 << qmMMBond[bondCounter][1] 
00892                 << " -> Value (distance or ratio): " << bondVal
00893                 << " (QM Group " << grpIter << " ID " << qmGrpID[grpIter] << ")"
00894                 << "\n" << endi ;
00895                 
00896                 qmDummyBondVal[bondCounter] = bondVal;
00897                 
00898                 // In case we are preparing for a mechanical embedding simulation
00899                 // with no point charges, populate the following vectors
00900                 if (qmMeNumBonds > 0) {
00901                     qmMeMMindx[bondCounter] = qmMMBond[bondCounter][0];
00902                     qmMeQMGrp[bondCounter] = qmGrpID[grpIter];
00903                 }
00904                 
00905                 bondCounter++ ;
00906             }
00907         }
00908         
00909     }
00910     
00911     
00914     
00915     current = NULL;
00916     if (qmNumBonds > 0)
00917         current = cfgList->find("QMLinkElement");
00918     
00919     int numParsedLinkElem = 0;
00920     for ( ; current != NULL; current = current->next ) {
00921         
00922         DebugM(3,"Parsing link atom element: " << current->data << "\n" );
00923         
00924         strVec = split( std::string(current->data) , " ");
00925         
00926         // We need the two atoms that compose the QM-MM bonds and 
00927         // then the element.
00928         if (strVec.size() != 3 && qmNumBonds > 1) {
00929             iout << iERROR << "Format error in QM link atom element selection: " 
00930             << current->data
00931             << "\n" << endi;
00932             NAMD_die("Error processing QM information.");
00933         }
00934         
00935         // If there is only one QM-MM bond, we can accept the element only.
00936         if (strVec.size() != 1 && strVec.size() != 3 && qmNumBonds == 1) {
00937             iout << iERROR << "Format error in QM link atom element selection: " 
00938             << current->data
00939             << "\n" << endi;
00940             NAMD_die("Error processing QM information.");
00941         }
00942         
00943         std::stringstream storConv ;
00944         int bondAtom1, bondAtom2;
00945         std::string linkElement ;
00946         
00947         if (strVec.size() == 1) {
00948             linkElement = strVec[0].substr(0,2);
00949         }
00950         else if (strVec.size() == 3) {
00951             
00952             storConv << strVec[0] ;
00953             storConv >> bondAtom1;
00954             if (storConv.fail()) {
00955                 iout << iERROR << "Error parsing link atom element selection: " 
00956                 << current->data
00957                 << "\n" << endi;
00958                 NAMD_die("Error processing QM information.");
00959             }
00960             
00961             storConv.clear() ;
00962             storConv << strVec[1];
00963             storConv >> bondAtom2;
00964             if (storConv.fail()) {
00965                 iout << iERROR << "Error parsing link atom element selection: " 
00966                 << current->data
00967                 << "\n" << endi;
00968                 NAMD_die("Error processing QM information.");
00969             }
00970             
00971             linkElement = strVec[2].substr(0,2);
00972         }
00973         
00974         numParsedLinkElem++;
00975         
00976         if (numParsedLinkElem>qmNumBonds) {
00977             iout << iERROR << "More link atom elements were set than there"
00978             " are QM-MM bonds.\n" << endi;
00979             NAMD_die("Error processing QM information.");
00980         }
00981         
00982         int bondIter;
00983         
00984         if (strVec.size() == 1) {
00985             bondIter = 0;
00986         }
00987         else if (strVec.size() == 3) {
00988             
00989             Bool foundBond = false;
00990             
00991             for (bondIter=0; bondIter<qmNumBonds; bondIter++) {
00992                 
00993                 if ( (qmMMBond[bondIter][0] == bondAtom1 &&
00994                     qmMMBond[bondIter][1] == bondAtom2 ) ||
00995                     (qmMMBond[bondIter][0] == bondAtom2 &&
00996                     qmMMBond[bondIter][1] == bondAtom1 ) ) {
00997                     
00998                     foundBond = true;
00999                     break;
01000                 }
01001             }
01002             
01003             if (! foundBond) {
01004                 iout << iERROR << "Error parsing link atom element selection: " 
01005                 << current->data
01006                 << "\n" << endi;
01007                 iout << iERROR << "Atom pair was not found to be a QM-MM bond.\n"
01008                 << endi;
01009                 NAMD_die("Error processing QM information.");
01010             }
01011         }
01012         
01013         strcpy(qmDummyElement[bondIter],linkElement.c_str());
01014         qmDummyElement[bondIter][2] = '\0';
01015         
01016         iout << iINFO << "Applying user defined link atom element "
01017         << qmDummyElement[bondIter] << " to QM-MM bond "
01018         << bondIter << ": " << qmMMBond[bondIter][0]
01019         << " - " << qmMMBond[bondIter][1]
01020         << "\n" << endi;
01021     }
01022     
01023     
01024     
01027     
01028     
01029     int32 **bondsWithAtomLocal = NULL ;
01030     int32 *byAtomSizeLocal = NULL;
01031     ObjectArena <int32 >* tmpArenaLocal = NULL;
01032     if (qmNumBonds > 0) {
01033         
01034         bondsWithAtomLocal = new int32 *[numAtoms];
01035         byAtomSizeLocal = new int32[numAtoms];
01036         tmpArenaLocal = new ObjectArena<int32>;
01037         
01038         //  Build the bond lists
01039        for (int i=0; i<numAtoms; i++)
01040        {
01041          byAtomSizeLocal[i] = 0;
01042        }
01043        for (int i=0; i<numRealBonds; i++)
01044        {
01045          byAtomSizeLocal[bonds[i].atom1]++;
01046          byAtomSizeLocal[bonds[i].atom2]++;
01047        }
01048        for (int i=0; i<numAtoms; i++)
01049        {
01050          bondsWithAtomLocal[i] = tmpArenaLocal->getNewArray(byAtomSizeLocal[i]+1);
01051          bondsWithAtomLocal[i][byAtomSizeLocal[i]] = -1;
01052          byAtomSizeLocal[i] = 0;
01053        }
01054        for (int i=0; i<numRealBonds; i++)
01055        {
01056          int a1 = bonds[i].atom1;
01057          int a2 = bonds[i].atom2;
01058          bondsWithAtomLocal[a1][byAtomSizeLocal[a1]++] = i;
01059          bondsWithAtomLocal[a2][byAtomSizeLocal[a2]++] = i;
01060        }
01061     }
01062     
01063     // In this loops we try to find other bonds in which the MM atoms from
01064     // QM-MM bonds may be involved. The other MM atoms (which we will call M2 and M3)
01065     // will be involved in charge manipulation. See ComputeQM.C for details.
01066     for (int qmBndIt = 0; qmBndIt < qmNumBonds; qmBndIt++) {
01067         
01068         // The charge targets are accumulated in a temporary vector and then 
01069         // transfered to an array that will be transmited to the ComputeQMMgr object.
01070         std::vector<int> chrgTrgt ;
01071         
01072         int MM1 = qmMMBond[qmBndIt][0], MM2, MM2BondIndx, MM3, MM3BondIndx;
01073         
01074         switch (simParams->qmBondScheme) {
01075             
01076             case QMSCHEMERCD:
01077             
01078             case QMSCHEMECS:
01079             {
01080                 // Selects ALL MM2 atoms.
01081                 for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
01082                     MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
01083                     
01084                     // Checks which of the atoms in the bond structure is the
01085                     // MM2 atom.
01086                     if (bonds[MM2BondIndx].atom1 == MM1)
01087                         MM2 = bonds[MM2BondIndx].atom2;
01088                     else
01089                         MM2 = bonds[MM2BondIndx].atom1;
01090                     
01091                     // In case the bonded atom is a QM atom,
01092                     // skips the index.
01093                     if (qmAtomGroup[MM2] > 0)
01094                         continue;
01095                     
01096                     chrgTrgt.push_back(MM2);
01097                 }
01098                 
01099             } break;
01100             
01101             case QMSCHEMEZ3:
01102             {
01103                 // Selects all MM3 atoms
01104                 for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
01105                     MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
01106                     
01107                     // Checks which of the atoms in the bond structure is the
01108                     // MM2 atom.
01109                     if (bonds[MM2BondIndx].atom1 == MM1)
01110                         MM2 = bonds[MM2BondIndx].atom2;
01111                     else
01112                         MM2 = bonds[MM2BondIndx].atom1;
01113                     
01114                     // In case the bonded atom is a QM atom,
01115                     // skips the index.
01116                     if (qmAtomGroup[MM2] > 0)
01117                         continue;
01118                     
01119                     for (int i=0; i<byAtomSizeLocal[MM2]; i++) {
01120                         MM3BondIndx = bondsWithAtomLocal[MM2][i] ;
01121                         
01122                         // Checks which of the atoms in the bond structure is the
01123                         // MM3 atom.
01124                         if (bonds[MM3BondIndx].atom1 == MM2)
01125                             MM3 = bonds[MM3BondIndx].atom2;
01126                         else
01127                             MM3 = bonds[MM3BondIndx].atom1;
01128                         
01129                         // In case the bonded atom is a QM atom,
01130                         // skips the index.
01131                         // We also keep the search from going back to MM1.
01132                         if (qmAtomGroup[MM3] > 0 || MM3 == MM1)
01133                             continue;
01134                         
01135                         chrgTrgt.push_back(MM3);
01136                     }
01137                     
01138                 }
01139                 
01140             };
01141             
01142             case QMSCHEMEZ2:
01143             {
01144                 // Selects all MM2 atoms
01145                 for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
01146                     MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
01147                     
01148                     // Checks which of the atoms in the bond structure is the
01149                     // MM2 atom.
01150                     if (bonds[MM2BondIndx].atom1 == MM1)
01151                         MM2 = bonds[MM2BondIndx].atom2;
01152                     else
01153                         MM2 = bonds[MM2BondIndx].atom1;
01154                     
01155                     // In case the bonded atom is a QM atom,
01156                     // skips the index.
01157                     if (qmAtomGroup[MM2] > 0)
01158                         continue;
01159                     
01160                     chrgTrgt.push_back(MM2);
01161                 }
01162                 
01163             };
01164             
01165             case QMSCHEMEZ1:
01166             {
01167                 // Selects all MM1 atoms
01168                 chrgTrgt.push_back(MM1);
01169             } break;
01170         }
01171         
01172         
01173         qmMMChargeTarget[qmBndIt] = new int[chrgTrgt.size()] ;
01174         qmMMNumTargs[qmBndIt] =  chrgTrgt.size();
01175         
01176         DebugM(3, "MM-QM bond " << qmBndIt << "; MM atom " 
01177         << qmMMBond[qmBndIt][0] << " conections: \n" );
01178         
01179         for (size_t i=0; i < chrgTrgt.size(); i++) {
01180             qmMMChargeTarget[qmBndIt][i] = chrgTrgt[i];
01181             DebugM(3,"MM Bonded to: " << chrgTrgt[i] << "\n" );
01182         }
01183         
01184         chrgTrgt.clear();
01185     }
01186     
01187     if (bondsWithAtomLocal != NULL)
01188         delete [] bondsWithAtomLocal;  bondsWithAtomLocal = 0;
01189     if (byAtomSizeLocal != NULL)
01190         delete [] byAtomSizeLocal;  byAtomSizeLocal = 0;
01191     if (tmpArenaLocal != NULL)
01192         delete tmpArenaLocal;  tmpArenaLocal = 0;
01193     
01194     
01197     
01198     
01199     if(simParams->qmLSSOn) {
01200         
01201         std::map<Real,int> grpLSSSize ;
01202         std::map<Real,int>::iterator itGrpSize;
01203         
01204         qmLSSTotalNumAtms = 0;
01205         qmLSSResidueSize = 0;
01206         
01207         if (simParams->qmLSSFreq == 0)
01208             qmLSSFreq = simParams->stepsPerCycle ;
01209         else
01210             qmLSSFreq = simParams->qmLSSFreq;
01211             
01212         #ifdef DEBUG_QM
01213         int resSize = -1;
01214         #endif
01215         
01216         std::map<Real, int> grpNumLSSRes;
01217         std::map<Real, int>::iterator itGrpNumRes;
01218         
01219         for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
01220             
01221             if (it->atmIDs.size() != it->size) {
01222                 iout << iERROR << "The number of atoms loaded for residue "
01223                 << it->resID << " does not match the expected for this residue type.\n" 
01224                 << endi;
01225                 NAMD_die("Error parsing data for LSS.");
01226             }
01227             
01228             qmLSSTotalNumAtms += it->size;
01229             
01230             #ifdef DEBUG_QM
01231             if (resSize < 0) resSize = it->size ;
01232             if (resSize > 0 and resSize != it->size) {
01233                 iout << iERROR << "The number of atoms loaded for residue "
01234                 << it->resID << " does not match previously loaded residues.\n" 
01235                 << endi;
01236                 NAMD_die("Error parsing data for LSS.");
01237             }
01238                 
01239 //             DebugM(3,"Residue " << it->resID << ": " << it->segName
01240 //             << " - from " << it->begAtmID << " with size "
01241 //             << it->size << " (QM ID: " << it->qmGrpID
01242 //             << ") has " << it->atmIDs.size() << " atoms: \n" ) ;
01243 //             for (int i=0; i<it->atmIDs.size(); i++) 
01244 //                 DebugM(3, it->atmIDs[i] << "\n" );
01245             #endif
01246             
01247             // Calculating total number of atoms per group
01248             itGrpSize = grpLSSSize.find(it->qmGrpID) ;
01249             if (itGrpSize != grpLSSSize.end())
01250                 itGrpSize->second += it->size;
01251             else
01252                 grpLSSSize.insert(std::pair<Real,int>(it->qmGrpID, it->size));
01253             
01254             // Calculating total number of solvent residues per group
01255             itGrpNumRes = grpNumLSSRes.find(it->qmGrpID) ;
01256             if (itGrpNumRes != grpNumLSSRes.end())
01257                 itGrpNumRes->second += 1;
01258             else
01259                 grpNumLSSRes.insert(std::pair<Real,int>(it->qmGrpID, 1));
01260         }
01261         
01262         qmLSSResidueSize = lssGrpRes.begin()->size ;
01263         
01264         qmLSSSize = new int[qmNumGrps];
01265         
01266         qmLSSIdxs = new int[qmLSSTotalNumAtms];
01267         int lssAtmIndx = 0;
01268         
01269         switch (simParams->qmLSSMode) {
01270         
01271         case QMLSSMODECOM:
01272         {
01273             
01274             qmLSSRefSize = new int[qmNumGrps];
01275             
01276             qmLSSMass = new Mass[qmLSSTotalNumAtms];
01277             
01278             qmLSSRefIDs = new int[totalNumRefAtms];
01279             qmLSSRefMass = new Mass[totalNumRefAtms];
01280             int lssRefIndx = 0;
01281             
01282             for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
01283                 
01284                 itGrpSize = grpNumLSSRes.find(qmGrpID[grpIndx]) ;
01285                 
01286                 if (itGrpSize != grpNumLSSRes.end())
01287                     qmLSSSize[grpIndx] =  itGrpSize->second;
01288                 else
01289                     qmLSSSize[grpIndx] =  0;
01290                 
01291                 for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
01292                     
01293                     if (it->qmGrpID == qmGrpID[grpIndx]) {
01294                         for (int i=0; i<it->atmIDs.size(); i++) {
01295                             qmLSSIdxs[lssAtmIndx] = it->atmIDs[i];
01296                             qmLSSMass[lssAtmIndx] = atoms[it->atmIDs[i]].mass;
01297                             lssAtmIndx++;
01298                         }
01299                     }
01300                 }
01301                 
01302                 DebugM(4, "QM group " << qmGrpID[grpIndx]
01303                 << " has " << qmLSSSize[grpIndx] << " solvent molecules, "
01304                 << "totalling " << grpLSSSize[qmGrpID[grpIndx]]
01305                 << " atoms (check " << lssAtmIndx << ").\n" );
01306                 
01307                 qmLSSRefSize[grpIndx] = lssGrpRefIDs[qmGrpID[grpIndx]].size();
01308                 for(int i=0; i < qmLSSRefSize[grpIndx]; i++) {
01309                     qmLSSRefIDs[lssRefIndx] = lssGrpRefIDs[qmGrpID[grpIndx]][i];
01310                     qmLSSRefMass[lssRefIndx] =  atoms[qmLSSRefIDs[lssRefIndx]].mass;
01311                     lssRefIndx++;
01312                 }
01313                 
01314                 DebugM(4, "QM group " << qmGrpID[grpIndx]
01315                 << " has " << qmLSSRefSize[grpIndx] << " non-solvent atoms for LSS.\n" );
01316             }
01317             
01318         } break ;
01319         
01320         case QMLSSMODEDIST:
01321         {
01322             for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
01323                 
01324                 itGrpSize = grpNumLSSRes.find(qmGrpID[grpIndx]) ;
01325                 
01326                 if (itGrpSize != grpNumLSSRes.end())
01327                     qmLSSSize[grpIndx] =  itGrpSize->second;
01328                 else
01329                     qmLSSSize[grpIndx] =  0;
01330                 
01331                 for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
01332                     
01333                     if (it->qmGrpID == qmGrpID[grpIndx]) {
01334                         for (int i=0; i<it->atmIDs.size(); i++) {
01335                             qmLSSIdxs[lssAtmIndx] = it->atmIDs[i];
01336                             lssAtmIndx++;
01337                         }
01338                     }
01339                 }
01340                 
01341                 DebugM(4, "QM group " << qmGrpID[grpIndx]
01342                 << " has " << qmLSSSize[grpIndx] << " solvent molecules, "
01343                 << "totalling " << grpLSSSize[qmGrpID[grpIndx]]
01344                 << " atoms (check " << lssAtmIndx << ").\n" );
01345             }
01346             
01347         } break ;
01348             
01349         }
01350     }
01351     
01352     
01355     
01356     
01357     PDB *customPCPDB;
01358     
01359     // In case we have a custom and fixed set of point charges for each QM group,
01360     // we process the files containing information.
01361     current = NULL;
01362     if (simParams->qmCustomPCSel) {
01363         current = cfgList->find("QMCustomPCFile");
01364     }
01365     
01366     std::map<Real,std::vector<int> > qmPCVecMap ;
01367     
01368     int numParsedPBDs = 0;
01369     for ( ; current != NULL; current = current->next ) {
01370         
01371         iout << iINFO << "Parsing QM Custom PC file " << current->data << "\n" << endi;
01372         customPCPDB = new PDB(current->data);
01373         
01374         if (customPCPDB->num_atoms() != numAtoms)
01375            NAMD_die("Number of atoms in QM Custom PC PDB file doesn't match coordinate PDB");
01376         
01377         std::vector< int > currPCSel ;
01378         Real currQMID = 0 ;
01379         int currGrpSize = 0 ;
01380         
01381         for (int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
01382             
01383             BigReal beta = customPCPDB->atom(atmInd)->temperaturefactor() ;
01384             BigReal occ = customPCPDB->atom(atmInd)->occupancy() ;
01385             
01386             if ( beta != 0 && occ != 0)
01387                 NAMD_die("An atom cannot be marked as QM and poitn charge simultaneously!");
01388             
01389             // If this is not a QM atom and 
01390             if (occ != 0) {
01391                 currPCSel.push_back(atmInd) ;
01392             }
01393             
01394             if (beta != 0) {
01395                 if (pdbP->atom(atmInd)->temperaturefactor() != beta)
01396                     NAMD_die("QM Group selection is different in reference PDB and Custom-PC PDB!");
01397                 
01398                 if (currQMID == 0) {
01399                     // If this is the first QM atom we find, keep the QM Group ID.
01400                     currQMID = beta;
01401                 }
01402                 else {
01403                     // For all other atoms, check if it is the same group. It must be!!
01404                     if (currQMID != beta)
01405                         NAMD_die("Found two different QM group IDs in this file!");
01406                 }
01407                 
01408                 currGrpSize++;
01409             }
01410             
01411         }
01412         
01413         if (currGrpSize != qmGrpSizeVec[ qmGrpIDMap[currQMID] ])
01414             NAMD_die("Reference PDB and Custom-PC PDB have different QM group sizes!") ;
01415         
01416         qmPCVecMap.insert(std::pair<Real,std::vector<int> >(
01417             currQMID, currPCSel ));
01418         
01419         numParsedPBDs++;
01420         delete customPCPDB;
01421     }
01422     
01423     delete pdbP;
01424     
01425     if (numParsedPBDs != qmNumGrps && simParams->qmCustomPCSel) {
01426         iout << iWARN << "The number of files provided for custom point "
01427         "charges is not equal to the number of QM groups!\n" << endi;
01428     }
01429     
01430     // Initializes an array with the number of Custom Point Charges per 
01431     // QM group.
01432     qmCustPCSizes = new int[qmNumGrps];
01433     for (int i=0; i<qmNumGrps; i++)
01434         qmCustPCSizes[i] = 0;
01435     
01436     qmTotCustPCs = 0;
01437     
01438     // Stores the size of each Custom PC vector in the array.
01439     // We may not have PCs for all QM groups.
01440     for (auto it = qmPCVecMap.begin(); it != qmPCVecMap.end(); it++) {
01441         qmTotCustPCs += (*it).second.size();
01442         int qmIndx = qmGrpIDMap[(*it).first];
01443         qmCustPCSizes[qmIndx] = (*it).second.size();
01444     }
01445     
01446     qmCustomPCIdxs = new int[qmTotCustPCs];
01447     
01448     if (simParams->qmCustomPCSel) {
01449         
01450         int auxIter = 0;
01451         for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
01452             
01453             Real currQMID = qmGrpID[grpIndx];
01454             
01455             iout << iINFO << "Loading " << qmPCVecMap[currQMID].size()
01456             << " custom point charges to QM Group " << grpIndx
01457             << " (ID: " << currQMID << ")\n" << endi;
01458             
01459             for (int pcIndx=0; pcIndx<qmPCVecMap[currQMID].size(); pcIndx++) {
01460                 qmCustomPCIdxs[auxIter] = qmPCVecMap[currQMID][pcIndx] ;
01461                 auxIter++;
01462             }
01463         }
01464     } 
01465     
01466     // Conditional SMD option
01467     qmcSMD = simParams->qmCSMD;
01468     if (qmcSMD) {
01469         read_qm_csdm_file(qmGrpIDMap);
01470     }
01471     
01474     
01475     if (simParams->qmElecEmbed) {
01476         // Modifies Atom charges for Electrostatic Embedding.
01477         // QM atoms cannot have charges in the standard location, to keep
01478         // NAMD from calculating electrostatic interactions between QM and MM atoms.
01479         // We handle electrostatics ourselves in ComputeQM.C and in special
01480         // modifications for PME.
01481         for (int i=0; i<qmNumQMAtoms; i++) {
01482             qmAtmChrg[i] = atoms[qmAtmIndx[i]].charge;
01483             atoms[qmAtmIndx[i]].charge = 0;
01484         }
01485     }
01486     
01487     
01488     if ( simParams->extraBondsOn) {
01489         // Lifted from Molecule::build_extra_bonds
01490         
01491         StringList *file = cfgList->find("extraBondsFile");
01492         
01493         char err_msg[512];
01494         int a1,a2; float k, ref;
01495         
01496         for ( ; file; file = file->next ) {  // loop over files
01497             FILE *f = fopen(file->data,"r");
01498 //             if ( ! f ) {
01499 //               sprintf(err_msg, "UNABLE TO OPEN EXTRA BONDS FILE %s", file->data);
01500 //               NAMD_err(err_msg);
01501 //             } else {
01502 //               iout << iINFO << "READING EXTRA BONDS FILE " << file->data <<"\n"<<endi;
01503 //             }
01504             
01505             while ( 1 ) {
01506               char buffer[512];
01507               int ret_code;
01508               do {
01509                 ret_code = NAMD_read_line(f, buffer);
01510               } while ( (ret_code==0) && (NAMD_blank_string(buffer)) );
01511               if (ret_code!=0) break;
01512 
01513               char type[512];
01514               sscanf(buffer,"%s",type);
01515               
01516               int badline = 0;
01517               if ( ! strncasecmp(type,"bond",4) ) {
01518                 if ( sscanf(buffer, "%s %d %d %f %f %s",
01519                     type, &a1, &a2, &k, &ref, err_msg) != 5 ) badline = 1;
01520                 
01521                 // If an extra bond is defined between QM atoms, we make
01522                 // note so that it wont be deleted when we delete bonded 
01523                 // interactions between QM atoms.
01524                 if( qmAtomGroup[a1] > 0 && qmAtomGroup[a2]) {
01525                     Bond tmp;
01526                     tmp.bond_type = 0;
01527                     tmp.atom1 = a1;  tmp.atom2 = a2;
01528                     qmExtraBonds.add(tmp);
01529                 }
01530                 
01531                 
01532               }else if ( ! strncasecmp(type,"#",1) ) {
01533                 continue;  // comment
01534               }
01535               
01536             }
01537             fclose(f);
01538         }
01539     }
01540     
01541     return;
01542     
01543 #endif // MEM_OPT_VERSION
01544 }

void Molecule::print_atoms ( Parameters params  ) 

Definition at line 5168 of file Molecule.C.

References DebugM, endi(), and Parameters::get_vdw_params().

Referenced by NamdState::loadStructure().

05169 {
05170 #ifdef MEM_OPT_VERSION
05171     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
05172 #else   
05173   register int i;
05174   Real sigma;
05175   Real epsilon;
05176   Real sigma14;
05177   Real epsilon14;
05178 
05179   DebugM(2,"ATOM LIST\n" \
05180       << "******************************************\n" \
05181                   << "NUM  NAME TYPE RES  MASS    CHARGE CHARGE   FEP-CHARGE"  \
05182       << "SIGMA   EPSILON SIGMA14 EPSILON14\n" \
05183         << endi);
05184 
05185   for (i=0; i<numAtoms; i++)
05186   {
05187     params->get_vdw_params(&sigma, &epsilon, &sigma14, &epsilon14, 
05188         atoms[i].vdw_type);
05189 
05190     DebugM(2,i+1 << " " << atomNames[i].atomname  \
05191               << " " << atomNames[i].atomtype << " " \
05192               << atomNames[i].resname  << " " << atoms[i].mass  \
05193         << " " << atoms[i].charge << " " << sigma \
05194         << " " << epsilon << " " << sigma14 \
05195         << " " << epsilon14 << "\n" \
05196         << endi);
05197   }
05198 #endif  
05199 }

void Molecule::print_bonds ( Parameters params  ) 

Definition at line 5211 of file Molecule.C.

References DebugM, endi(), and Parameters::get_bond_params().

Referenced by NamdState::loadStructure().

05212 {
05213 #ifdef MEM_OPT_VERSION
05214     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
05215 #else   
05216   register int i;
05217   Real k;
05218   Real x0;
05219 
05220   DebugM(2,"BOND LIST\n" << "********************************\n" \
05221       << "ATOM1 ATOM2 TYPE1 TYPE2      k        x0" \
05222       << endi);
05223 
05224   for (i=0; i<numBonds; i++)
05225   {
05226     params->get_bond_params(&k, &x0, bonds[i].bond_type);
05227 
05228     DebugM(2,bonds[i].atom1+1 << " " \
05229        << bonds[i].atom2+1 << " "   \
05230        << atomNames[bonds[i].atom1].atomtype << " "  \
05231        << atomNames[bonds[i].atom2].atomtype << " " << k \
05232        << " " << x0 << endi);
05233   }
05234   
05235 #endif  
05236 }

void Molecule::print_exclusions (  ) 

Definition at line 5248 of file Molecule.C.

References DebugM, and endi().

Referenced by NamdState::loadStructure().

05249 {
05250 #ifdef MEM_OPT_VERSION
05251     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
05252 #else
05253   register int i;
05254 
05255   DebugM(2,"EXPLICIT EXCLUSION LIST\n" \
05256       << "********************************\n" \
05257             << "ATOM1 ATOM2 " \
05258       << endi);
05259 
05260   for (i=0; i<numExclusions; i++)
05261   {
05262     DebugM(2,exclusions[i].atom1+1 << "  " \
05263        << exclusions[i].atom2+1 << endi);
05264   }
05265 #endif
05266 }

void Molecule::print_go_params (  ) 

Definition at line 548 of file GoMolecule.C.

References DebugM, go_array, j, MAX_GO_CHAINS, and NumGoChains.

00549 {
00550   int i;
00551   int j;
00552   int index;
00553 
00554   DebugM(3,NumGoChains << " Go PARAMETERS 3\n" \
00555          << "*****************************************" << std::endl);
00556 
00557   for (i=0; i<NumGoChains; i++) {
00558     for (j=0; j<NumGoChains; j++) {
00559       index = (i * MAX_GO_CHAINS) + j;
00560       //  Real epsilon;    // Epsilon
00561       //  Real exp_a;      // First exponent for attractive L-J term
00562       //  Real exp_b;      // Second exponent for attractive L-J term
00563       //  Real sigmaRep;   // Sigma for repulsive term
00564       //  Real epsilonRep; // Epsilon for replusive term
00565       DebugM(3,"Go index=(" << i << "," << j << ") epsilon=" << go_array[index].epsilon \
00566              << " exp_a=" << go_array[index].exp_a << " exp_b=" << go_array[index].exp_b \
00567              << " exp_rep=" << go_array[index].exp_rep << " sigmaRep=" \
00568              << go_array[index].sigmaRep << " epsilonRep=" << go_array[index].epsilonRep \
00569              << " cutoff=" << go_array[index].cutoff << std::endl);
00570     }
00571   }
00572 
00573 }

void Molecule::print_go_sigmas (  ) 

Definition at line 1134 of file GoMolecule.C.

References DebugM, goSigmaIndices, goSigmas, j, numAtoms, and numGoAtoms.

01135 {
01136   int i;  //  Counter
01137   int j;  //  Counter
01138   Real sigma;
01139 
01140   DebugM(3,"GO SIGMA ARRAY\n" \
01141          << "***************************" << std::endl);
01142   DebugM(3, "numGoAtoms: " << numGoAtoms << std::endl);
01143 
01144   if (goSigmaIndices == NULL) {
01145     DebugM(3, "GO SIGMAS HAVE NOT BEEN BUILT" << std::endl);
01146     return;
01147   }
01148 
01149   for (i=0; i<numAtoms; i++) {
01150     for (j=0; j<numAtoms; j++) {
01151       if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 ) {
01152         //DebugM(3, "i: " << i << ", j: " << j << std::endl);
01153         sigma = goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]];
01154         if (sigma > 0.0) {
01155           DebugM(3, "(" << i << "," << j << ") - +" << sigma << " ");
01156         }
01157         else {
01158           DebugM(3, "(" << i << "," << j << ") - " << sigma << " ");
01159         }
01160       } else {
01161         //DebugM(3, "0 ");
01162       }
01163     }
01164     if ( goSigmaIndices[i] != -1 ) {
01165       DebugM(3, "-----------" << std::endl);
01166     }
01167   }
01168   return;
01169 }

void Molecule::put_stir_startTheta ( Real  theta,
int  atomnum 
) const [inline]

Definition at line 1295 of file Molecule.h.

01296   {
01297     stirParams[stirIndexes[atomnum]].startTheta = theta;
01298   }

void Molecule::read_go_file ( char *  fname  ) 

Definition at line 113 of file GoMolecule.C.

References go_val::cutoff, DebugM, endi(), go_val::epsilon, go_val::epsilonRep, go_val::exp_a, go_val::exp_b, go_val::exp_rep, FALSE, go_array, go_indices, iout, iWARN(), j, MAX_GO_CHAINS, MAX_RESTRICTIONS, NAMD_blank_string(), NAMD_die(), NAMD_find_first_word(), NAMD_read_line(), NumGoChains, go_val::restrictions, go_val::sigmaRep, and TRUE.

Referenced by build_go_params().

00115 {
00116 
00117   int i;                   // Counter
00118   int j;                   // Counter
00119   int  par_type=0;         //  What type of parameter are we currently
00120                            //  dealing with? (vide infra)
00121   // JLai -- uncommented
00122   int  skipline;           //  skip this line?
00123   int  skipall = 0;        //  skip rest of file;
00124   char buffer[512];           //  Buffer to store each line of the file
00125   char first_word[512];           //  First word of the current line
00126   int read_count = 0;      //  Count of input parameters on a given line
00127   int chain1 = 0;          //  First chain type for pair interaction
00128   int chain2 = 0;          //  Second chain type for pair interaction
00129   int int1;                //  First parameter int
00130   int int2;                //  Second parameter int
00131   Real r1;                 //  Parameter Real
00132   char in2[512];           //  Second parameter word
00133   GoValue *goValue1 = NULL;    //  current GoValue for loading parameters
00134   GoValue *goValue2 = NULL;    //  other current GoValue for loading parameters
00135   Bool sameGoChain = FALSE;    //  whether the current GoValue is within a chain or between chains
00136   int restrictionCount = 0;    //  number of Go restrictions set for a given chain pair
00137   FILE *pfile;                 //  File descriptor for the parameter file
00138 
00139   /*  Check to make sure that we haven't previously been told     */
00140   /*  that all the files were read                                */
00141   /*if (AllFilesRead)
00142     {
00143     NAMD_die("Tried to read another parameter file after being told that all files were read!");
00144     }*/
00145   
00146   /*  Initialize go_indices  */
00147   for (i=0; i<MAX_GO_CHAINS+1; i++) {
00148     go_indices[i] = -1;
00149   }
00150 
00151   /*  Initialize go_array   */
00152   for (i=0; i<MAX_GO_CHAINS*MAX_GO_CHAINS; i++) {
00153     go_array[i].epsilon = 1.25;
00154     go_array[i].exp_a = 12;
00155     go_array[i].exp_b = 6;
00156     go_array[i].exp_rep = 12;
00157     go_array[i].sigmaRep = 2.25;
00158     go_array[i].epsilonRep = 0.03;
00159     go_array[i].cutoff = 4.0;
00160     for (j=0; j<MAX_RESTRICTIONS; j++) {
00161       go_array[i].restrictions[j] = -1;
00162     }
00163   }
00164 
00165   /*  Try and open the file                                        */
00166   if ( (pfile = fopen(fname, "r")) == NULL)
00167     {
00168       char err_msg[256];
00169       
00170       sprintf(err_msg, "UNABLE TO OPEN GO PARAMETER FILE %s\n", fname);
00171       NAMD_die(err_msg);
00172     }
00173   
00174   /*  Keep reading in lines until we hit the EOF                        */
00175   while (NAMD_read_line(pfile, buffer) != -1)
00176     {
00177       /*  Get the first word of the line                        */
00178       NAMD_find_first_word(buffer, first_word);
00179       skipline=0;
00180       
00181       /*  First, screen out things that we ignore.                   */   
00182       /*  blank lines, lines that start with '!' or '*', lines that  */
00183       /*  start with "END".                                          */
00184       if (!NAMD_blank_string(buffer) &&
00185           (strncmp(first_word, "!", 1) != 0) &&
00186           (strncmp(first_word, "*", 1) != 0) &&
00187           (strncasecmp(first_word, "END", 3) != 0))
00188         {
00189           if ( skipall ) {
00190             iout << iWARN << "SKIPPING PART OF GO PARAMETER FILE AFTER RETURN STATEMENT\n" << endi;
00191             break;
00192           }
00193           /*  Now, determine the apropriate parameter type.   */
00194           if (strncasecmp(first_word, "chaintypes", 10)==0)
00195             {
00196               read_count=sscanf(buffer, "%s %d %d\n", first_word, &int1, &int2);
00197               if (read_count != 3) {
00198                 char err_msg[512];
00199                 sprintf(err_msg, "UNKNOWN PARAMETER IN GO PARAMETER FILE %s\nLINE=*%s*\nread_count=%d, int1=%d, int2=%d", fname, buffer, read_count, int1, int2);
00200                 NAMD_die(err_msg);
00201               }
00202               chain1 = int1;
00203               chain2 = int2;
00204               if (chain1 < 1 || chain1 > MAX_GO_CHAINS ||
00205                   chain2 < 1 || chain2 > MAX_GO_CHAINS) {
00206                 char err_msg[512];
00207                 sprintf(err_msg, "GO PARAMETER FILE: CHAIN INDEX MUST BE [1-%d] %s\nLINE=*%s*\nread_count=%d, int1=%d, int2=%d", MAX_GO_CHAINS, fname, buffer, read_count, int1, int2);
00208                 NAMD_die(err_msg);
00209               }
00210               if (go_indices[chain1] == -1) {
00211                 go_indices[chain1] = NumGoChains;
00212                 NumGoChains++;
00213               }
00214               if (go_indices[chain2] == -1) {
00215                 go_indices[chain2] = NumGoChains;
00216                 NumGoChains++;
00217               }
00218               if (chain1 == chain2) {
00219                 sameGoChain = TRUE;
00220               } else {
00221                 sameGoChain = FALSE;
00222               }
00223               //goValue = &go_array[(chain1 * MAX_GO_CHAINS) + chain2];
00224               goValue1 = &go_array[(chain1*MAX_GO_CHAINS) + chain2];
00225               goValue2 = &go_array[(chain2*MAX_GO_CHAINS) + chain1];
00226 #if CODE_REDUNDANT
00227               goValue1 = &go_array[(go_indices[chain1]*MAX_GO_CHAINS) + go_indices[chain2]];
00228               goValue2 = &go_array[(go_indices[chain2]*MAX_GO_CHAINS) + go_indices[chain1]];
00229 #endif
00230               restrictionCount = 0;    //  restrictionCount applies to each chain pair separately
00231             }
00232           else if (strncasecmp(first_word, "epsilonRep", 10)==0)
00233             {
00234               read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
00235               if (read_count != 2) {}
00236               goValue1->epsilonRep = r1;
00237               if (!sameGoChain) {
00238                 goValue2->epsilonRep = r1;
00239               }
00240             }
00241           else if (strncasecmp(first_word, "epsilon", 7)==0)
00242             {
00243               // Read in epsilon
00244               read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
00245               if (read_count != 2) {}
00246               goValue1->epsilon = r1;
00247               if (!sameGoChain) {
00248                 goValue2->epsilon = r1;
00249               }
00250             }
00251           else if (strncasecmp(first_word, "exp_a", 5)==0)
00252             {
00253               read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
00254               if (read_count != 2) {}
00255               goValue1->exp_a = int1;
00256               if (!sameGoChain) {
00257                 goValue2->exp_a = int1;
00258               }
00259             }
00260           else if (strncasecmp(first_word, "exp_b", 5)==0)
00261             {
00262               read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
00263               if (read_count != 2) {}
00264               goValue1->exp_b = int1;
00265               if (!sameGoChain) {
00266                 goValue2->exp_b = int1;
00267               }
00268             }
00269           else if (strncasecmp(first_word, "exp_rep", 5)==0)
00270             {
00271               read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
00272               if (read_count != 2) {}
00273               goValue1->exp_rep = int1;
00274               if (!sameGoChain) {
00275                 goValue2->exp_rep = int1;
00276               }
00277             }
00278           else if (strncasecmp(first_word, "exp_Rep", 5)==0)
00279             {
00280               read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
00281               if (read_count != 2) {}
00282               goValue1->exp_rep = int1;
00283               if (!sameGoChain) {
00284               goValue2->exp_rep = int1;
00285               }
00286             }
00287           else if (strncasecmp(first_word, "sigmaRep", 8)==0)
00288             {
00289               read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
00290               if (read_count != 2) {}
00291               goValue1->sigmaRep = r1;
00292               if (!sameGoChain) {
00293                 goValue2->sigmaRep = r1;
00294               }
00295             }
00296           else if (strncasecmp(first_word, "cutoff", 6)==0)
00297             {
00298               read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
00299               if (read_count != 2) {}
00300               goValue1->cutoff = r1;
00301               if (!sameGoChain) {
00302                 goValue2->cutoff = r1;
00303               }
00304             }
00305           else if (strncasecmp(first_word, "restriction", 10)==0)
00306             {
00307               read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
00308               if (read_count != 2) {}
00309               if (int1 < 0) {
00310                 DebugM(3, "ERROR: residue restriction value must be nonnegative.  int1=" << int1 << "\n");
00311               }
00312               if (!sameGoChain) {
00313                 //goValue2->restrictions[restrictionCount] = int1;
00314                 DebugM(3, "ERROR: residue restrictions should not be defined between two separate chains.  chain1=" << chain1 << ", chain2=" << chain2 << "\n");
00315               }
00316               else {
00317                 goValue1->restrictions[restrictionCount] = int1;
00318               }
00319               restrictionCount++;
00320             }
00321           else if (strncasecmp(first_word, "return", 4)==0)
00322             {
00323               skipall=8;
00324               skipline=1;
00325             }        
00326           else // if (par_type == 0)
00327             {
00328               /*  This is an unknown paramter.        */
00329               /*  This is BAD                                */
00330               char err_msg[512];
00331               
00332               sprintf(err_msg, "UNKNOWN PARAMETER IN GO PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00333               NAMD_die(err_msg);
00334             }
00335         }
00336       else
00337         {
00338           skipline=1;
00339         }
00340     }
00341   
00342   /*  Close the file  */
00343   fclose(pfile);
00344   
00345   return;
00346 }

void Molecule::read_parm ( Ambertoppar  ) 
void Molecule::receive_GoMolecule ( MIStream msg  ) 

Definition at line 1744 of file GoMolecule.C.

References atomChainTypes, go_val::cutoff, go_val::epsilon, go_val::epsilonRep, go_val::exp_a, go_val::exp_b, go_val::exp_rep, MIStream::get(), go_array, go_indices, goCoordinates, SimParameters::goForcesOn, goIndxLJA, goIndxLJB, SimParameters::goMethod, goNumLJPair, goResidIndices, goResids, goSigmaIndices, goSigmaPairA, goSigmaPairB, goSigmas, goWithinCutoff, j, MAX_GO_CHAINS, MAX_RESTRICTIONS, NAMD_die(), numAtoms, numGoAtoms, NumGoChains, pointerToGoBeg, pointerToGoEnd, go_val::restrictions, and go_val::sigmaRep.

01744                                                {
01745       // Ported by JLai -- Original by JE
01746       // JE - receive Go info
01747       Real *a1, *a2, *a3, *a4;
01748       int *i1, *i2, *i3, *i4;
01749       int maxGoChainsSqr = MAX_GO_CHAINS*MAX_GO_CHAINS;  // JE JLai Go code
01750       msg->get(NumGoChains);
01751       
01752       if (NumGoChains) {
01753         //go_indices = new int[MAX_GO_CHAINS+1];
01754         //go_array = new GoValue[MAX_GO_CHAINS*MAX_GO_CHAINS];
01755         
01756         //      int go_indices[MAX_GO_CHAINS+1];        //  Indices from chainIDs to go_array      
01757         //      GoValue go_array[MAX_GO_CHAINS*MAX_GO_CHAINS];   //  Array of Go params
01758         msg->get(MAX_GO_CHAINS+1,go_indices);
01759         
01760         a1 = new Real[maxGoChainsSqr];
01761         a2 = new Real[maxGoChainsSqr];
01762         a3 = new Real[maxGoChainsSqr];
01763         a4 = new Real[maxGoChainsSqr];
01764         i1 = new int[maxGoChainsSqr];
01765         i2 = new int[maxGoChainsSqr];
01766         i3 = new int[maxGoChainsSqr];
01767         i4 = new int[maxGoChainsSqr*MAX_RESTRICTIONS];
01768         
01769         if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) || 
01770              (i1 == NULL) || (i2 == NULL) || (i3 == NULL) || (i4 == NULL) )
01771           {
01772             NAMD_die("memory allocation failed in Molecule::send_Molecule");
01773           }
01774         
01775         msg->get(maxGoChainsSqr, a1);
01776         msg->get(maxGoChainsSqr, a2);
01777         msg->get(maxGoChainsSqr, a3);
01778         msg->get(maxGoChainsSqr, a4);
01779         msg->get(maxGoChainsSqr, i1);
01780         msg->get(maxGoChainsSqr, i2);
01781         msg->get(maxGoChainsSqr, i3);
01782         msg->get(maxGoChainsSqr*MAX_RESTRICTIONS, i4);
01783         
01784         for (int i=0; i<maxGoChainsSqr; i++) {
01785           go_array[i].epsilon = a1[i];
01786           go_array[i].sigmaRep = a2[i];
01787           go_array[i].epsilonRep = a3[i];
01788           go_array[i].cutoff = a4[i];
01789           go_array[i].exp_a = i1[i];
01790           go_array[i].exp_b = i2[i];
01791           go_array[i].exp_rep = i3[i];
01792           for (int j=0; j<MAX_RESTRICTIONS; j++) {
01793             go_array[i].restrictions[j] = i4[i*MAX_RESTRICTIONS + j];
01794           }
01795         }
01796         
01797         delete [] a1;
01798         delete [] a2;
01799         delete [] a3;
01800         delete [] a4;
01801         delete [] i1;
01802         delete [] i2;
01803         delete [] i3;
01804         delete [] i4;
01805 
01806         //msg->get(MAX_GO_CHAINS*MAX_GO_CHAINS, (char*)go_array);
01807         
01808         /*DebugM(3,"NumGoChains:" << NumGoChains << std::endl);
01809           for (int ii=0; ii<MAX_GO_CHAINS; ii++) {
01810           for (int jj=0; jj<MAX_GO_CHAINS; jj++) {
01811           DebugM(3,"go_array[" << ii << "][" << jj << "]:" << go_array[ii][jj] << std::endl);
01812           }
01813           }
01814           for (int ii=0; ii<MAX_GO_CHAINS+1; ii++) {
01815           DebugM(3,"go_indices[" << ii << "]:" << go_indices[ii] << std::endl);
01816           }*/
01817       }
01818 
01819       if (simParams->goForcesOn) {
01820         switch(simParams->goMethod) {
01821         case 1:
01822           msg->get(numGoAtoms);
01823           //printf("Deleting goSigmaIndiciesA\n");
01824           delete [] goSigmaIndices;
01825           goSigmaIndices = new int32[numAtoms];
01826           //printf("Deleting atomChainTypesA\n");
01827           delete [] atomChainTypes;
01828           atomChainTypes = new int32[numGoAtoms];
01829           //printf("Deleting goSigmasA\n");
01830           delete [] goSigmas;
01831           goSigmas = new Real[numGoAtoms*numGoAtoms];
01832           //printf("Deleting goWithinCutoffA\n"); 
01833           delete [] goWithinCutoff;
01834           goWithinCutoff = new bool[numGoAtoms*numGoAtoms];
01835           msg->get(numAtoms, goSigmaIndices);
01836           msg->get(numGoAtoms, atomChainTypes);
01837           msg->get(numGoAtoms*numGoAtoms, goSigmas);
01838           msg->get(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
01839           break;          
01840         case 2: //GSR
01841           msg->get(numGoAtoms);
01842           delete [] pointerToGoBeg;
01843           pointerToGoBeg = new int[numAtoms];
01844           msg->get(numAtoms,pointerToGoBeg);
01845           delete [] pointerToGoEnd;
01846           pointerToGoEnd = new int[numAtoms];
01847           msg->get(numAtoms,pointerToGoEnd);
01848           delete [] goSigmaIndices;
01849           goSigmaIndices = new int32[numAtoms];
01850           msg->get(numAtoms,goSigmaIndices);
01851           delete [] goResidIndices;
01852           goResidIndices = new int32[numAtoms];
01853           msg->get(numAtoms,goResidIndices);      
01854           delete [] atomChainTypes;
01855           atomChainTypes = new int32[numGoAtoms];
01856           msg->get(numGoAtoms,atomChainTypes);
01857           msg->get(goNumLJPair);
01858           delete [] goIndxLJA;
01859           goIndxLJA = new int[goNumLJPair];
01860           msg->get(goNumLJPair,goIndxLJA);
01861           delete [] goIndxLJB;
01862           goIndxLJB = new int[goNumLJPair];
01863           msg->get(goNumLJPair,goIndxLJB);
01864           delete [] goSigmaPairA;
01865           goSigmaPairA = new double[goNumLJPair];
01866           msg->get(goNumLJPair,goSigmaPairA);
01867           delete [] goSigmaPairB;
01868           goSigmaPairB = new double[goNumLJPair];
01869           msg->get(goNumLJPair,goSigmaPairB);  
01870           break;
01871         case 3:
01872           msg->get(numGoAtoms);
01873           //printf("Deleting goSigmaIndiciesB\n");
01874           delete [] goSigmaIndices;
01875           goSigmaIndices = new int32[numAtoms];
01876           //printf("Deleting atomChainTypesB\n");
01877           delete [] atomChainTypes;
01878           atomChainTypes = new int32[numGoAtoms];
01879           //delete [] goSigmas;
01880           //goSigmas = new Real[numGoAtoms*numGoAtoms];
01881           //delete [] goWithinCutoff;
01882           //goWithinCutoff = new bool[numGoAtoms*numGoAtoms];
01883           //printf("Deleting goCoordinatesB\n");
01884           delete [] goCoordinates;
01885           goCoordinates = new Real[numGoAtoms*3];
01886           //printf("Deleting goResidsB\n");
01887           delete [] goResids;
01888           goResids = new int[numGoAtoms];
01889           msg->get(numAtoms, goSigmaIndices);
01890           msg->get(numGoAtoms, atomChainTypes);
01891           //msg->get(numGoAtoms*numGoAtoms, goSigmas);
01892           //msg->get(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
01893           msg->get(numGoAtoms*3, goCoordinates);
01894           msg->get(numGoAtoms, goResids);
01895           break;
01896         }
01897       }
01898 
01899       delete msg;
01900 
01901 }

void Molecule::receive_Molecule ( MIStream msg  ) 

Definition at line 5623 of file Molecule.C.

References SimParameters::alchOn, atomNamePool, atomSigPool, consForce, consForceIndexes, SimParameters::consForceOn, consTorqueIndexes, SimParameters::consTorqueOn, consTorqueParams, SimParameters::constraintsOn, DebugM, SimParameters::excludeFromPressure, SimParameters::fixedAtomsOn, gA, MIStream::get(), Parameters::get_num_vdw_params(), GridforceGrid::get_total_grids(), ObjectArena< Type >::getNewArray(), giSigma1, giSigma2, gMu1, gMu2, SimParameters::goGroPair, gRepulsive, gromacsPair_type, indxGaussA, indxGaussB, indxLJA, indxLJB, is_drude_psf, is_lonepairs_psf, isBFactorValid, isOccupancyValid, SimParameters::langevinOn, SimParameters::LCPOOn, SimParameters::lesOn, maxHydrogenGroupSize, maxMigrationGroupSize, SimParameters::mgridforceOn, SimParameters::movDragOn, numAcceptors, numCalcAngles, numCalcBonds, numCalcCrossterms, numCalcDihedrals, numCalcExclusions, numCalcFullExclusions, numCalcImpropers, numCalcLJPair, numConsForce, numConsTorque, numConstraints, numDonors, numExPressureAtoms, numFepFinal, numFepInitial, numFixedRigidBonds, numGaussPair, numGridforceGrids, numGridforces, numHydrogenGroups, numMigrationGroups, numMovDrag, numRotDrag, numStirredAtoms, pairC12, pairC6, SimParameters::pairInteractionOn, pointerToGaussBeg, pointerToGaussEnd, pointerToLJBeg, pointerToLJEnd, SimParameters::qmForcesOn, SimParameters::rotDragOn, SimParameters::soluteScalingOn, ss_index, ss_num_vdw_params, ss_vdw_type, SimParameters::stirOn, SimParameters::tCoupleOn, and GridforceGrid::unpack_grid().

Referenced by Node::resendMolecule().

05623                                             {
05624   //  Get the atom information
05625   msg->get(numAtoms);
05626 
05627 #ifdef MEM_OPT_VERSION
05628 //in the memory optimized version, only the atom signatures are recved
05629 //from the master Node. --Chao Mei
05630 
05631   msg->get(massPoolSize);
05632   if(atomMassPool) delete [] atomMassPool;
05633   atomMassPool = new Real[massPoolSize];
05634   msg->get(massPoolSize, atomMassPool);
05635 
05636   msg->get(chargePoolSize);
05637   if(atomChargePool) delete [] atomChargePool;
05638   atomChargePool = new Real[chargePoolSize];
05639   msg->get(chargePoolSize, atomChargePool);
05640 
05641   //get atoms' signatures
05642   msg->get(atomSigPoolSize);
05643   if(atomSigPool) delete [] atomSigPool;
05644   atomSigPool = new AtomSignature[atomSigPoolSize];
05645   for(int i=0; i<atomSigPoolSize; i++)
05646       atomSigPool[i].unpack(msg);
05647 
05648   //get exclusions' signatures
05649   msg->get(exclSigPoolSize);
05650   if(exclSigPool) delete [] exclSigPool;
05651   exclSigPool = new ExclusionSignature[exclSigPoolSize];
05652   for(int i=0; i<exclSigPoolSize; i++)
05653       exclSigPool[i].unpack(msg);
05654  
05655   msg->get(numHydrogenGroups);      
05656   msg->get(maxHydrogenGroupSize);      
05657   msg->get(numMigrationGroups);      
05658   msg->get(maxMigrationGroupSize);      
05659   msg->get(isOccupancyValid);
05660   msg->get(isBFactorValid);
05661 
05662    //get names for atoms
05663   msg->get(atomNamePoolSize);
05664   atomNamePool = new char *[atomNamePoolSize];
05665   for(int i=0; i<atomNamePoolSize;i++) {
05666     int len;
05667     msg->get(len);
05668     atomNamePool[i] = nameArena->getNewArray(len+1);
05669     msg->get(len, atomNamePool[i]);
05670   }
05671   
05672   if(simParams->fixedAtomsOn){
05673     int numFixedAtomsSet;
05674     msg->get(numFixedAtoms);
05675     msg->get(numFixedAtomsSet);
05676     fixedAtomsSet = new AtomSetList(numFixedAtomsSet);
05677     msg->get(numFixedAtomsSet*sizeof(AtomSet), (char *)(fixedAtomsSet->begin()));
05678   } 
05679 
05680   if(simParams->constraintsOn){
05681     int numConstrainedAtomsSet;
05682     msg->get(numConstraints);
05683     msg->get(numConstrainedAtomsSet);
05684     constrainedAtomsSet = new AtomSetList(numConstrainedAtomsSet);
05685     msg->get(numConstrainedAtomsSet*sizeof(AtomSet), (char *)(constrainedAtomsSet->begin()));
05686   } 
05687 
05688 #else
05689   delete [] atoms;
05690   atoms= new Atom[numAtoms];  
05691   msg->get(numAtoms*sizeof(Atom), (char*)atoms);
05692 
05693   //  Get the bond information
05694   msg->get(numRealBonds);
05695   msg->get(numBonds);    
05696   if (numBonds)
05697   {
05698     delete [] bonds;
05699     bonds=new Bond[numBonds]; 
05700     msg->get(numBonds*sizeof(Bond), (char*)bonds);
05701   }  
05702   
05703   //  Get the angle information
05704   msg->get(numAngles);  
05705   if (numAngles)
05706   {
05707     delete [] angles;
05708     angles=new Angle[numAngles];  
05709     msg->get(numAngles*sizeof(Angle), (char*)angles);
05710   }  
05711   
05712   //  Get the dihedral information
05713   msg->get(numDihedrals);    
05714   if (numDihedrals)
05715   {
05716     delete [] dihedrals;
05717     dihedrals=new Dihedral[numDihedrals];  
05718     msg->get(numDihedrals*sizeof(Dihedral), (char*)dihedrals);
05719   }  
05720   
05721   //  Get the improper information
05722   msg->get(numImpropers);
05723   if (numImpropers)
05724   {
05725     delete [] impropers;
05726     impropers=new Improper[numImpropers];  
05727     msg->get(numImpropers*sizeof(Improper), (char*)impropers);
05728   }
05729   
05730   //  Get the crossterm information
05731   msg->get(numCrossterms);
05732   if (numCrossterms)
05733   {
05734     delete [] crossterms;
05735     crossterms=new Crossterm[numCrossterms];  
05736     msg->get(numCrossterms*sizeof(Crossterm), (char*)crossterms);
05737   }
05738   
05739   //  Get the hydrogen bond donors
05740   msg->get(numDonors);  
05741   if (numDonors)
05742   {
05743     delete [] donors;
05744     donors=new Bond[numDonors];  
05745     msg->get(numDonors*sizeof(Bond), (char*)donors);
05746   }
05747   
05748   //  Get the hydrogen bond acceptors
05749   msg->get(numAcceptors);  
05750   if (numAcceptors)
05751   {
05752     delete [] acceptors;
05753     acceptors=new Bond[numAcceptors];  
05754     msg->get(numAcceptors*sizeof(Bond), (char*)acceptors);
05755   }
05756   
05757   //  Get the exclusion information 
05758   msg->get(numExclusions);  
05759   if (numExclusions)
05760   {
05761     delete [] exclusions;
05762     exclusions=new Exclusion[numExclusions];  
05763     msg->get(numExclusions*sizeof(Exclusion), (char*)exclusions);
05764   }
05765         
05766       //  Get the constraint information, if they are active
05767       if (simParams->constraintsOn)
05768       {
05769          msg->get(numConstraints);
05770 
05771          delete [] consIndexes;
05772          consIndexes = new int32[numAtoms];
05773          
05774          msg->get(numAtoms, consIndexes);
05775          
05776          if (numConstraints)
05777          {
05778            delete [] consParams;
05779            consParams = new ConstraintParams[numConstraints];
05780       
05781            msg->get(numConstraints*sizeof(ConstraintParams), (char*)consParams);
05782          }
05783       }
05784 #endif
05785 
05786       /* BEGIN gf */
05787       if (simParams->mgridforceOn)
05788       {
05789          DebugM(3, "Receiving gridforce info\n");
05790          
05791          msg->get(numGridforceGrids);
05792          
05793          DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n");
05794          
05795          delete [] numGridforces;
05796          numGridforces = new int[numGridforceGrids];
05797          
05798          delete [] gridfrcIndexes;      // Should I be deleting elements of these first?
05799          delete [] gridfrcParams;
05800          delete [] gridfrcGrid;
05801          gridfrcIndexes = new int32*[numGridforceGrids];
05802          gridfrcParams = new GridforceParams*[numGridforceGrids];
05803          gridfrcGrid = new GridforceGrid*[numGridforceGrids];
05804          
05805          int grandTotalGrids = 0;
05806          for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
05807              msg->get(numGridforces[gridnum]);
05808              
05809              gridfrcIndexes[gridnum] = new int32[numAtoms];
05810              msg->get(numAtoms, gridfrcIndexes[gridnum]);
05811          
05812              if (numGridforces[gridnum])
05813              {
05814                  gridfrcParams[gridnum] = new GridforceParams[numGridforces[gridnum]];
05815                  msg->get(numGridforces[gridnum]*sizeof(GridforceParams), (char*)gridfrcParams[gridnum]);
05816              }
05817              
05818              gridfrcGrid[gridnum] = GridforceGrid::unpack_grid(gridnum, msg);
05819              
05820              grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
05821          }
05822       }
05823       /* END gf */
05824       
05825       //  Get the stirring information, if stirring is  active
05826       if (simParams->stirOn)
05827       {
05828          msg->get(numStirredAtoms);
05829 
05830          delete [] stirIndexes;
05831          stirIndexes = new int32[numAtoms];
05832          
05833          msg->get(numAtoms, stirIndexes);
05834          
05835          if (numStirredAtoms)
05836          {
05837            delete [] stirParams;
05838            stirParams = new StirParams[numStirredAtoms];
05839       
05840            msg->get(numStirredAtoms*sizeof(StirParams), (char*)stirParams);
05841          }
05842       }
05843       
05844       //  Get the moving drag information, if it is active
05845       if (simParams->movDragOn) {
05846          msg->get(numMovDrag);
05847          delete [] movDragIndexes;
05848          movDragIndexes = new int32[numAtoms];
05849          msg->get(numAtoms, movDragIndexes);
05850          if (numMovDrag)
05851          {
05852            delete [] movDragParams;
05853            movDragParams = new MovDragParams[numMovDrag];
05854            msg->get(numMovDrag*sizeof(MovDragParams), (char*)movDragParams);
05855          }
05856       }
05857       
05858       //  Get the rotating drag information, if it is active
05859       if (simParams->rotDragOn) {
05860          msg->get(numRotDrag);
05861          delete [] rotDragIndexes;
05862          rotDragIndexes = new int32[numAtoms];
05863          msg->get(numAtoms, rotDragIndexes);
05864          if (numRotDrag)
05865          {
05866            delete [] rotDragParams;
05867            rotDragParams = new RotDragParams[numRotDrag];
05868            msg->get(numRotDrag*sizeof(RotDragParams), (char*)rotDragParams);
05869          }
05870       }
05871       
05872       //  Get the "constant" torque information, if it is active
05873       if (simParams->consTorqueOn) {
05874          msg->get(numConsTorque);
05875          delete [] consTorqueIndexes;
05876          consTorqueIndexes = new int32[numAtoms];
05877          msg->get(numAtoms, consTorqueIndexes);
05878          if (numConsTorque)
05879          {
05880            delete [] consTorqueParams;
05881            consTorqueParams = new ConsTorqueParams[numConsTorque];
05882            msg->get(numConsTorque*sizeof(ConsTorqueParams), (char*)consTorqueParams);
05883          }
05884       }
05885       
05886       // Get the constant force information, if it's active
05887       if (simParams->consForceOn)
05888       { msg->get(numConsForce);
05889         delete [] consForceIndexes;
05890         consForceIndexes = new int32[numAtoms];
05891         msg->get(numAtoms, consForceIndexes);
05892         if (numConsForce)
05893         { delete [] consForce;
05894           consForce = new Vector[numConsForce];
05895           msg->get(numConsForce*sizeof(Vector), (char*)consForce);
05896         }
05897       }
05898 
05899       if (simParams->excludeFromPressure) {
05900         exPressureAtomFlags = new int32[numAtoms];
05901         msg->get(numExPressureAtoms);
05902         msg->get(numAtoms, exPressureAtomFlags);
05903       }
05904 
05905 #ifndef MEM_OPT_VERSION
05906       //  Get the langevin parameters, if they are active
05907       if (simParams->langevinOn || simParams->tCoupleOn)
05908       {
05909         delete [] langevinParams;
05910         langevinParams = new Real[numAtoms];
05911 
05912         msg->get(numAtoms, langevinParams);
05913       }
05914 
05915       //  Get the fixed atoms, if they are active
05916       if (simParams->fixedAtomsOn)
05917       {
05918         delete [] fixedAtomFlags;
05919         fixedAtomFlags = new int32[numAtoms];
05920 
05921         msg->get(numFixedAtoms);
05922         msg->get(numAtoms, fixedAtomFlags);
05923         msg->get(numFixedRigidBonds);
05924       }
05925 
05926       if (simParams->qmForcesOn)
05927       {
05928         if( qmAtomGroup != 0)
05929             delete [] qmAtomGroup;
05930         qmAtomGroup = new Real[numAtoms];
05931         
05932         msg->get(numAtoms, qmAtomGroup);
05933         
05934         msg->get(qmNumQMAtoms);
05935         
05936         if( qmAtmChrg != 0)
05937             delete [] qmAtmChrg;
05938         qmAtmChrg = new Real[qmNumQMAtoms];
05939         
05940         msg->get(qmNumQMAtoms, qmAtmChrg);
05941         
05942         if( qmAtmIndx != 0)
05943             delete [] qmAtmIndx;
05944         qmAtmIndx = new int[qmNumQMAtoms];
05945         
05946         msg->get(qmNumQMAtoms, qmAtmIndx);
05947         
05948         msg->get(qmNoPC);
05949         
05950         msg->get(qmNumBonds);
05951         
05952         msg->get(qmMeNumBonds);
05953         
05954         if( qmMeMMindx != 0)
05955             delete [] qmMeMMindx;
05956         qmMeMMindx = new int[qmMeNumBonds];
05957         
05958         msg->get(qmMeNumBonds, qmMeMMindx);
05959         
05960         if( qmMeQMGrp != 0)
05961             delete [] qmMeQMGrp;
05962         qmMeQMGrp = new Real[qmMeNumBonds];
05963         
05964         msg->get(qmMeNumBonds, qmMeQMGrp);
05965         
05966         msg->get(qmPCFreq);
05967         
05968         msg->get(qmNumGrps);
05969         
05970         if( qmGrpID != 0)
05971             delete [] qmGrpID;
05972         qmGrpID = new Real[qmNumGrps];
05973         msg->get(qmNumGrps, qmGrpID);
05974         
05975         if( qmCustPCSizes != 0)
05976             delete [] qmCustPCSizes;
05977         qmCustPCSizes = new int[qmNumGrps];
05978         msg->get(qmNumGrps, qmCustPCSizes);
05979         
05980         msg->get(qmTotCustPCs);
05981         
05982         if( qmCustomPCIdxs != 0)
05983             delete [] qmCustomPCIdxs;
05984         qmCustomPCIdxs = new int[qmTotCustPCs];
05985         msg->get(qmTotCustPCs, qmCustomPCIdxs);
05986       }
05987     
05988 //fepb
05989       //receive fep atom info
05990       if (simParams->alchOn || simParams->lesOn || simParams->pairInteractionOn) {
05991         delete [] fepAtomFlags;
05992         fepAtomFlags = new unsigned char[numAtoms];
05993 
05994         msg->get(numFepInitial);
05995         msg->get(numFepFinal);
05996         msg->get(numAtoms*sizeof(unsigned char), (char*)fepAtomFlags);
05997       }
05998 //fepe
05999 
06000 //soluteScaling
06001       if (simParams->soluteScalingOn) {
06002         delete [] ssAtomFlags;
06003         delete [] ss_vdw_type;
06004         delete [] ss_index;
06005         ssAtomFlags = new unsigned char[numAtoms];
06006         ss_vdw_type = new int [params->get_num_vdw_params()];
06007         ss_index = new int [numAtoms];
06008         msg->get(numAtoms*sizeof(unsigned char), (char*)ssAtomFlags);
06009         msg->get(ss_num_vdw_params);
06010         msg->get(params->get_num_vdw_params()*sizeof(int), (char*)ss_vdw_type);
06011         msg->get(numAtoms*sizeof(int), (char*)ss_index);
06012       }
06013 //soluteScaling
06014 #ifdef OPENATOM_VERSION
06015       // This needs to be refactored into its own version
06016       if (simParams->openatomOn) {
06017         delete [] fepAtomFlags;
06018         fepAtomFlags = new unsigned char[numAtoms];
06019 
06020         msg->get(numFepInitial);
06021         msg->get(numAtoms*sizeof(unsigned char), (char*)fepAtomFlags);
06022 #endif //OPENATOM_VERSION
06023 
06024       // DRUDE: receive data read from PSF
06025       msg->get(is_lonepairs_psf);
06026       if (is_lonepairs_psf) {
06027         msg->get(numLphosts);
06028         delete[] lphosts;
06029         lphosts = new Lphost[numLphosts];
06030         msg->get(numLphosts*sizeof(Lphost), (char*)lphosts);
06031       }
06032       msg->get(is_drude_psf);
06033       if (is_drude_psf) {
06034         delete[] drudeConsts;
06035         drudeConsts = new DrudeConst[numAtoms];
06036         msg->get(numAtoms*sizeof(DrudeConst), (char*)drudeConsts);
06037         msg->get(numAnisos);
06038         delete[] anisos;
06039         anisos = new Aniso[numAnisos];
06040         msg->get(numAnisos*sizeof(Aniso), (char*)anisos);
06041       }
06042       // DRUDE
06043 
06044   //LCPO
06045   if (simParams->LCPOOn) {
06046     delete [] lcpoParamType;
06047     lcpoParamType = new int[numAtoms];
06048     msg->get(numAtoms, (int*)lcpoParamType);
06049   }
06050 
06051   //Receive GromacsPairStuff -- JLai
06052 
06053   if (simParams->goGroPair) {
06054     msg->get(numLJPair);
06055     delete [] indxLJA;
06056     indxLJA = new int[numLJPair];
06057     msg->get(numLJPair,indxLJA);
06058     delete [] indxLJB;
06059     indxLJB = new int[numLJPair];
06060     msg->get(numLJPair,indxLJB);
06061     delete [] pairC6;
06062     pairC6 = new Real[numLJPair];    
06063     msg->get(numLJPair,pairC6);
06064     delete [] pairC12;
06065     pairC12 = new Real[numLJPair];
06066     msg->get(numLJPair,pairC12);
06067     delete [] gromacsPair_type;
06068     gromacsPair_type = new int[numLJPair];
06069     msg->get(numLJPair,gromacsPair_type);
06070     delete [] pointerToLJBeg;
06071     pointerToLJBeg = new int[numAtoms];
06072     msg->get((numAtoms),pointerToLJBeg);
06073     delete [] pointerToLJEnd;
06074     pointerToLJEnd = new int[numAtoms];
06075     msg->get((numAtoms),pointerToLJEnd);
06076     // JLai
06077     delete [] gromacsPair;
06078     gromacsPair = new GromacsPair[numLJPair];
06079     for(int i=0; i < numLJPair; i++) {
06080         gromacsPair[i].atom1 = indxLJA[i];
06081         gromacsPair[i].atom2 = indxLJB[i];
06082         gromacsPair[i].pairC6  = pairC6[i];
06083         gromacsPair[i].pairC12 = pairC12[i];
06084         gromacsPair[i].gromacsPair_type = gromacsPair_type[i];
06085     }
06086     //
06087     msg->get(numGaussPair);
06088     delete [] indxGaussA;
06089     indxGaussA = new int[numGaussPair];
06090     msg->get(numGaussPair,indxGaussA);
06091     delete [] indxGaussB;
06092     indxGaussB = new int[numGaussPair];
06093     msg->get(numGaussPair,indxGaussB);
06094     delete [] gA;
06095     gA = new Real[numGaussPair];
06096     msg->get(numGaussPair,gA);
06097     delete [] gMu1;
06098     gMu1 = new Real[numGaussPair];
06099     msg->get(numGaussPair,gMu1);
06100     delete [] giSigma1;
06101     giSigma1 = new Real[numGaussPair];
06102     msg->get(numGaussPair,giSigma1);
06103     delete [] gMu2;
06104     gMu2 = new Real[numGaussPair];
06105     msg->get(numGaussPair,gMu2);
06106     delete [] giSigma2;
06107     giSigma2 = new Real[numGaussPair];
06108     msg->get(numGaussPair,giSigma2);
06109     delete [] gRepulsive;
06110     gRepulsive = new Real[numGaussPair];
06111     msg->get(numGaussPair,gRepulsive);
06112     delete [] pointerToGaussBeg;
06113     pointerToGaussBeg = new int[numAtoms];
06114     msg->get((numAtoms),pointerToGaussBeg);
06115     delete [] pointerToGaussEnd;
06116     pointerToGaussEnd = new int[numAtoms];
06117     msg->get((numAtoms),pointerToGaussEnd);
06118     //
06119   }
06120 #endif
06121 
06122       //  Now free the message 
06123       delete msg;
06124 
06125 #ifdef MEM_OPT_VERSION
06126 
06127       build_excl_check_signatures();
06128 
06129     //set num{Calc}Tuples(Bonds,...,Impropers) to 0
06130     numBonds = numCalcBonds = 0;
06131     numAngles = numCalcAngles = 0;
06132     numDihedrals = numCalcDihedrals = 0;
06133     numImpropers = numCalcImpropers = 0;
06134     numCrossterms = numCalcCrossterms = 0;
06135     numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;  
06136     // JLai
06137     numLJPair = numCalcLJPair = 0;
06138     // End of JLai
06139 
06140 #else
06141 
06142       //  analyze the data and find the status of each atom
06143       build_atom_status();
06144       build_lists_by_atom();      
06145 
06146       
06147 #endif
06148 }

void Molecule::reloadCharges ( float  charge[],
int  n 
)

Referenced by Node::reloadCharges().

Real Molecule::rigid_bond_length ( int  atomnum  )  const [inline]

Definition at line 1440 of file Molecule.h.

Referenced by WorkDistrib::createAtomLists(), and outputCompressedFile().

01441   {
01442     return(rigidBondLengths[atomnum]);
01443   }

void Molecule::send_GoMolecule ( MOStream msg  ) 

Definition at line 1635 of file GoMolecule.C.

References atomChainTypes, go_val::cutoff, MOStream::end(), go_val::epsilon, go_val::epsilonRep, go_val::exp_a, go_val::exp_b, go_val::exp_rep, go_array, go_indices, goCoordinates, SimParameters::goForcesOn, goIndxLJA, goIndxLJB, SimParameters::goMethod, goNumLJPair, goResidIndices, goResids, goSigmaIndices, goSigmaPairA, goSigmaPairB, goSigmas, goWithinCutoff, j, MAX_GO_CHAINS, MAX_RESTRICTIONS, NAMD_die(), numAtoms, numGoAtoms, NumGoChains, pointerToGoBeg, pointerToGoEnd, MOStream::put(), go_val::restrictions, and go_val::sigmaRep.

01635                                             {
01636   Real *a1, *a2, *a3, *a4;
01637   int *i1, *i2, *i3, *i4;
01638   int maxGoChainsSqr = MAX_GO_CHAINS*MAX_GO_CHAINS;  // JE JLai Go code
01639   msg->put(NumGoChains);
01640   
01641   if (NumGoChains) {
01642     //      int go_indices[MAX_GO_CHAINS+1];        //  Indices from chainIDs to go_array      
01643     //      GoValue go_array[MAX_GO_CHAINS*MAX_GO_CHAINS];   //  Array of Go params
01644     msg->put(MAX_GO_CHAINS+1,go_indices);
01645 
01646     a1 = new Real[maxGoChainsSqr];
01647     a2 = new Real[maxGoChainsSqr];
01648     a3 = new Real[maxGoChainsSqr];
01649     a4 = new Real[maxGoChainsSqr];
01650     i1 = new int[maxGoChainsSqr];
01651     i2 = new int[maxGoChainsSqr];
01652     i3 = new int[maxGoChainsSqr];
01653     i4 = new int[maxGoChainsSqr*MAX_RESTRICTIONS];
01654 
01655     if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) || 
01656          (i1 == NULL) || (i2 == NULL) || (i3 == NULL) || (i4 == NULL) )
01657     {
01658       NAMD_die("memory allocation failed in Molecules::send_Molecules");
01659     }
01660 
01661     for (int i=0; i<maxGoChainsSqr; i++) {
01662       a1[i] = go_array[i].epsilon;
01663       a2[i] = go_array[i].sigmaRep;
01664       a3[i] = go_array[i].epsilonRep;
01665       a4[i] = go_array[i].cutoff;
01666       i1[i] = go_array[i].exp_a;
01667       i2[i] = go_array[i].exp_b;
01668       i3[i] = go_array[i].exp_rep;
01669       for (int j=0; j<MAX_RESTRICTIONS; j++) {
01670         i4[i*MAX_RESTRICTIONS + j] = go_array[i].restrictions[j];
01671       }
01672     }
01673 
01674     msg->put(maxGoChainsSqr, a1);
01675     msg->put(maxGoChainsSqr, a2);
01676     msg->put(maxGoChainsSqr, a3);
01677     msg->put(maxGoChainsSqr, a4);
01678     msg->put(maxGoChainsSqr, i1);
01679     msg->put(maxGoChainsSqr, i2);
01680     msg->put(maxGoChainsSqr, i3);
01681     msg->put(maxGoChainsSqr*MAX_RESTRICTIONS, i4);
01682 
01683     delete [] a1;
01684     delete [] a2;
01685     delete [] a3;
01686     delete [] a4;
01687     delete [] i1;
01688     delete [] i2;
01689     delete [] i3;
01690     delete [] i4;
01691   }
01692 
01693   //Ported JLai
01694   if (simParams->goForcesOn) {
01695     switch(simParams->goMethod) {
01696     case 1:
01697       msg->put(numGoAtoms);
01698       msg->put(numAtoms, goSigmaIndices);
01699       msg->put(numGoAtoms, atomChainTypes);
01700       msg->put(numGoAtoms*numGoAtoms, goSigmas);
01701       msg->put(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
01702       // printf("Molecule.C sending atomChainTypes %d %d \n", numGoAtoms, atomChainTypes);
01703       break;
01704     case 2: //GSS
01705       msg->put(numGoAtoms);
01706       msg->put(numAtoms,pointerToGoBeg);
01707       msg->put(numAtoms,pointerToGoEnd);
01708       msg->put(numAtoms,goSigmaIndices);
01709       msg->put(numAtoms,goResidIndices);
01710       msg->put(numGoAtoms,atomChainTypes);
01711       msg->put(goNumLJPair);
01712       msg->put(goNumLJPair,goIndxLJA);
01713       msg->put(goNumLJPair,goIndxLJB);
01714       msg->put(goNumLJPair,goSigmaPairA);
01715       msg->put(goNumLJPair,goSigmaPairB);
01716       break;
01717     case 3:
01718       msg->put(numGoAtoms);
01719       msg->put(numAtoms, goSigmaIndices);
01720       msg->put(numGoAtoms, atomChainTypes);
01721       //msg->put(numGoAtoms*numGoAtoms, goSigmas);
01722       //msg->put(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
01723       msg->put(numGoAtoms*3, goCoordinates);
01724       msg->put(numGoAtoms, goResids);
01725       break;
01726     }
01727   } 
01728 
01729   msg->end();
01730   delete msg;
01731 }

void Molecule::send_Molecule ( MOStream msg  ) 

Definition at line 5279 of file Molecule.C.

References SimParameters::alchOn, atomNamePool, atomSigPool, consForce, consForceIndexes, SimParameters::consForceOn, consTorqueIndexes, SimParameters::consTorqueOn, consTorqueParams, SimParameters::constraintsOn, DebugM, MOStream::end(), endi(), SimParameters::excludeFromPressure, SimParameters::fixedAtomsOn, gA, Parameters::get_num_vdw_params(), giSigma1, giSigma2, gMu1, gMu2, SimParameters::goGroPair, gRepulsive, gromacsPair_type, indxGaussA, indxGaussB, indxLJA, indxLJB, is_drude_psf, is_lonepairs_psf, isBFactorValid, isOccupancyValid, SimParameters::langevinOn, SimParameters::LCPOOn, SimParameters::lesOn, maxHydrogenGroupSize, maxMigrationGroupSize, SimParameters::mgridforceOn, SimParameters::movDragOn, numAcceptors, numCalcAngles, numCalcBonds, numCalcCrossterms, numCalcDihedrals, numCalcExclusions, numCalcFullExclusions, numCalcImpropers, numCalcLJPair, numConsForce, numConsTorque, numConstraints, numDonors, numExPressureAtoms, numFepFinal, numFepInitial, numFixedRigidBonds, numGaussPair, numGridforceGrids, numGridforces, numHydrogenGroups, numMigrationGroups, numMovDrag, numRotDrag, numStirredAtoms, GridforceGrid::pack_grid(), pairC12, pairC6, SimParameters::pairInteractionOn, pointerToGaussBeg, pointerToGaussEnd, pointerToLJBeg, pointerToLJEnd, MOStream::put(), SimParameters::qmForcesOn, SimParameters::rotDragOn, SimParameters::soluteScalingOn, ss_index, ss_num_vdw_params, ss_vdw_type, SimParameters::stirOn, and SimParameters::tCoupleOn.

Referenced by Node::resendMolecule().

05279                                          {
05280 #ifdef MEM_OPT_VERSION
05281 //in the memory optimized version, only the atom signatures are broadcast
05282 //to other Nodes. --Chao Mei
05283 
05284   msg->put(numAtoms);
05285 
05286   msg->put(massPoolSize);
05287   msg->put(massPoolSize, atomMassPool);
05288 
05289   msg->put(chargePoolSize);
05290   msg->put(chargePoolSize, atomChargePool); 
05291 
05292   //put atoms' signatures
05293   msg->put(atomSigPoolSize);
05294   for(int i=0; i<atomSigPoolSize; i++)
05295       atomSigPool[i].pack(msg);
05296 
05297   //put atom's exclusion signatures
05298   msg->put(exclSigPoolSize);
05299   for(int i=0; i<exclSigPoolSize; i++)
05300       exclSigPool[i].pack(msg);
05301 
05302   msg->put(numHydrogenGroups);      
05303   msg->put(maxHydrogenGroupSize);      
05304   msg->put(numMigrationGroups);      
05305   msg->put(maxMigrationGroupSize);            
05306   msg->put(isOccupancyValid);
05307   msg->put(isBFactorValid);
05308   
05309   //put names for atoms
05310   msg->put(atomNamePoolSize);
05311   for(int i=0; i<atomNamePoolSize;i++) {
05312     int len = strlen(atomNamePool[i]);
05313     msg->put(len);
05314     msg->put(len*sizeof(char), atomNamePool[i]);
05315   } 
05316   
05317   if(simParams->fixedAtomsOn){
05318     int numFixedAtomsSet = fixedAtomsSet->size();
05319     msg->put(numFixedAtoms);
05320     msg->put(numFixedAtomsSet);
05321     msg->put(numFixedAtomsSet*sizeof(AtomSet), (char *)(fixedAtomsSet->begin()));
05322   }
05323 
05324   if (simParams->constraintsOn) {
05325     int numConstrainedAtomsSet = constrainedAtomsSet->size();
05326     msg->put(numConstraints);
05327     msg->put(numConstrainedAtomsSet);
05328     msg->put(numConstrainedAtomsSet*sizeof(AtomSet), (char *)(constrainedAtomsSet->begin()));
05329   }
05330     
05331 #else
05332   msg->put(numAtoms);
05333   msg->put(numAtoms*sizeof(Atom), (char*)atoms);
05334   
05335   //  Send the bond information
05336   msg->put(numRealBonds);
05337   msg->put(numBonds);
05338  
05339   if (numBonds)
05340   {
05341     msg->put(numBonds*sizeof(Bond), (char*)bonds);
05342   }
05343 
05344   //  Send the angle information
05345   msg->put(numAngles);  
05346   if (numAngles)
05347   {
05348     msg->put(numAngles*sizeof(Angle), (char*)angles);
05349   }  
05350 
05351   //  Send the dihedral information
05352   msg->put(numDihedrals);
05353   if (numDihedrals)
05354   {
05355     msg->put(numDihedrals*sizeof(Dihedral), (char*)dihedrals);
05356   }  
05357 
05358   //  Send the improper information
05359   msg->put(numImpropers);  
05360   if (numImpropers)
05361   {
05362     msg->put(numImpropers*sizeof(Improper), (char*)impropers);
05363   }
05364 
05365   //  Send the crossterm information
05366   msg->put(numCrossterms);
05367   if (numCrossterms)
05368   {
05369     msg->put(numCrossterms*sizeof(Crossterm), (char*)crossterms);
05370   }
05371 
05372   // send the hydrogen bond donor information
05373   msg->put(numDonors);
05374   if(numDonors)
05375   {
05376     msg->put(numDonors*sizeof(Bond), (char*)donors);
05377   }
05378 
05379   // send the hydrogen bond acceptor information
05380   msg->put(numAcceptors);
05381   if(numAcceptors)
05382   {
05383     msg->put(numAcceptors*sizeof(Bond), (char*)acceptors);
05384   }
05385 
05386   //  Send the exclusion information  
05387   msg->put(numExclusions);
05388   if (numExclusions)
05389   {
05390     msg->put(numExclusions*sizeof(Exclusion), (char*)exclusions);
05391   }      
05392   //  Send the constraint information, if used
05393   if (simParams->constraintsOn)
05394   {
05395      msg->put(numConstraints);
05396      
05397      msg->put(numAtoms, consIndexes);
05398      
05399      if (numConstraints)
05400      {
05401        msg->put(numConstraints*sizeof(ConstraintParams), (char*)consParams);
05402      }
05403   }
05404 #endif
05405   
05406   /* BEGIN gf */
05407   // Send the gridforce information, if used
05408   if (simParams->mgridforceOn)
05409   {
05410     DebugM(3, "Sending gridforce info\n" << endi);
05411     msg->put(numGridforceGrids);
05412     
05413     for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
05414       msg->put(numGridforces[gridnum]);
05415       msg->put(numAtoms, gridfrcIndexes[gridnum]);
05416       if (numGridforces[gridnum])
05417       {
05418        msg->put(numGridforces[gridnum]*sizeof(GridforceParams), (char*)gridfrcParams[gridnum]);
05419       }
05420       GridforceGrid::pack_grid(gridfrcGrid[gridnum], msg);
05421     }
05422   }
05423   /* END gf */
05424   
05425   //  Send the stirring information, if used
05426   if (simParams->stirOn)
05427   {
05428      //CkPrintf ("DEBUG: putting numStirredAtoms..\n");
05429      msg->put(numStirredAtoms);
05430      //CkPrintf ("DEBUG: putting numAtoms,stirIndexes.. numAtoms=%d\n",numStirredAtoms);
05431      msg->put(numAtoms, stirIndexes);
05432      //CkPrintf ("DEBUG: if numStirredAtoms..\n");
05433      if (numStirredAtoms)
05434      {
05435        //CkPrintf ("DEBUG: big put, with (char*)stirParams\n");
05436        msg->put(numStirredAtoms*sizeof(StirParams), (char*)stirParams);
05437      }
05438   }
05439   
05440   
05441   //  Send the moving drag information, if used
05442   if (simParams->movDragOn) {
05443      msg->put(numMovDrag);
05444      msg->put(numAtoms, movDragIndexes);
05445      if (numMovDrag)
05446      {
05447        msg->put(numMovDrag*sizeof(MovDragParams), (char*)movDragParams);
05448      }
05449   }
05450   
05451   //  Send the rotating drag information, if used
05452   if (simParams->rotDragOn) {
05453      msg->put(numRotDrag);
05454      msg->put(numAtoms, rotDragIndexes);
05455      if (numRotDrag)
05456      {
05457        msg->put(numRotDrag*sizeof(RotDragParams), (char*)rotDragParams);
05458      }
05459   }
05460   
05461   //  Send the "constant" torque information, if used
05462   if (simParams->consTorqueOn) {
05463      msg->put(numConsTorque);
05464      msg->put(numAtoms, consTorqueIndexes);
05465      if (numConsTorque)
05466      {
05467        msg->put(numConsTorque*sizeof(ConsTorqueParams), (char*)consTorqueParams);
05468      }
05469   }
05470   
05471   // Send the constant force information, if used
05472   if (simParams->consForceOn)
05473   { msg->put(numConsForce);
05474     msg->put(numAtoms, consForceIndexes);
05475     if (numConsForce)
05476       msg->put(numConsForce*sizeof(Vector), (char*)consForce);
05477   }
05478   
05479   if (simParams->excludeFromPressure) {
05480     msg->put(numExPressureAtoms);
05481     msg->put(numAtoms, exPressureAtomFlags);
05482   }
05483   
05484 #ifndef MEM_OPT_VERSION
05485   //  Send the langevin parameters, if active
05486   if (simParams->langevinOn || simParams->tCoupleOn)
05487   {
05488     msg->put(numAtoms, langevinParams);
05489   }
05490   
05491   //  Send fixed atoms, if active
05492   if (simParams->fixedAtomsOn)
05493   {
05494     msg->put(numFixedAtoms);
05495     msg->put(numAtoms, fixedAtomFlags);
05496   msg->put(numFixedRigidBonds);
05497   }
05498   
05499   if (simParams->qmForcesOn)
05500   {
05501     msg->put(numAtoms, qmAtomGroup);
05502     msg->put(qmNumQMAtoms);
05503     msg->put(qmNumQMAtoms, qmAtmChrg);
05504     msg->put(qmNumQMAtoms, qmAtmIndx);
05505     msg->put(qmNoPC);
05506     msg->put(qmNumBonds);
05507     msg->put(qmMeNumBonds);
05508     msg->put(qmMeNumBonds, qmMeMMindx);
05509     msg->put(qmMeNumBonds, qmMeQMGrp);
05510     msg->put(qmPCFreq);
05511     msg->put(qmNumGrps);
05512     msg->put(qmNumGrps, qmGrpID);
05513     msg->put(qmNumGrps, qmCustPCSizes);
05514     msg->put(qmTotCustPCs);
05515     msg->put(qmTotCustPCs, qmCustomPCIdxs);
05516   }
05517   
05518   //fepb
05519   // send fep atom info
05520   if (simParams->alchOn || simParams->lesOn || simParams->pairInteractionOn) {
05521     msg->put(numFepInitial);
05522     msg->put(numFepFinal);
05523     msg->put(numAtoms*sizeof(char), (char*)fepAtomFlags);
05524   }
05525   //fepe
05526 
05527   if (simParams->soluteScalingOn) {
05528     msg->put(numAtoms*sizeof(char), (char*)ssAtomFlags);
05529     msg->put(ss_num_vdw_params);
05530     msg->put(params->get_num_vdw_params()*sizeof(int), (char*)ss_vdw_type);
05531     msg->put(numAtoms*sizeof(int), (char*)ss_index);
05532   }
05533 
05534   #ifdef OPENATOM_VERSION
05535   // needs to be refactored into its own openatom version
05536   if (simParams->openatomOn ) {
05537     msg->put(numFepInitial);
05538     msg->put(numAtoms*sizeof(char), (char*)fepAtomFlags);
05539   }
05540   #endif //OPENATOM_VERSION
05541   
05542   // DRUDE: send data read from PSF
05543   msg->put(is_lonepairs_psf);
05544   if (is_lonepairs_psf) {
05545     msg->put(numLphosts);
05546     msg->put(numLphosts*sizeof(Lphost), (char*)lphosts);
05547   }
05548   msg->put(is_drude_psf);
05549   if (is_drude_psf) {
05550     msg->put(numAtoms*sizeof(DrudeConst), (char*)drudeConsts);
05551     msg->put(numAnisos);
05552     msg->put(numAnisos*sizeof(Aniso), (char*)anisos);
05553   }
05554   // DRUDE
05555 
05556   //LCPO
05557   if (simParams->LCPOOn) {
05558     msg->put(numAtoms, (int*)lcpoParamType);
05559   }
05560   
05561   //Send GromacsPairStuff -- JLai
05562   if (simParams->goGroPair) {
05563     msg->put(numLJPair);
05564     msg->put(numLJPair,indxLJA);
05565     msg->put(numLJPair,indxLJB);
05566     msg->put(numLJPair,pairC6);
05567     msg->put(numLJPair,pairC12);
05568     msg->put(numLJPair,gromacsPair_type);
05569     msg->put((numAtoms),pointerToLJBeg);
05570     msg->put((numAtoms),pointerToLJEnd);
05571     msg->put(numGaussPair);
05572     msg->put(numGaussPair,indxGaussA);
05573     msg->put(numGaussPair,indxGaussB);
05574     msg->put(numGaussPair,gA);
05575     msg->put(numGaussPair,gMu1);
05576     msg->put(numGaussPair,giSigma1);
05577     msg->put(numGaussPair,gMu2);
05578     msg->put(numGaussPair,giSigma2);
05579     msg->put(numGaussPair,gRepulsive);
05580     msg->put((numAtoms),pointerToGaussBeg);
05581     msg->put((numAtoms),pointerToGaussEnd);
05582   }
05583 #endif
05584 
05585   // Broadcast the message to the other nodes
05586   msg->end();
05587   delete msg;
05588 
05589 #ifdef MEM_OPT_VERSION
05590 
05591   build_excl_check_signatures();
05592 
05593   //set num{Calc}Tuples(Bonds,...,Impropers) to 0
05594   numBonds = numCalcBonds = 0;
05595   numAngles = numCalcAngles = 0;
05596   numDihedrals = numCalcDihedrals = 0;
05597   numImpropers = numCalcImpropers = 0;
05598   numCrossterms = numCalcCrossterms = 0;
05599   numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;  
05600   // JLai
05601   numLJPair = numCalcLJPair = 0;
05602   // End of JLai
05603 
05604 #else
05605 
05606   //  Now build arrays of indexes into these arrays by atom      
05607   build_lists_by_atom();
05608 
05609 #endif
05610 }

int Molecule::set_gridfrc_grid ( int  gridnum,
GridforceGrid grid 
) [inline]

Definition at line 1272 of file Molecule.h.

References numGridforceGrids.

Referenced by Node::reloadGridforceGrid().

01273   {
01274       if (grid && gridnum >= 0 && gridnum < numGridforceGrids) {
01275           gridfrcGrid[gridnum] = grid;
01276           return 0;
01277       } else {
01278           return -1;
01279       }
01280   }

void Molecule::set_qm_replaceAll ( Bool  newReplaceAll  )  [inline]

Definition at line 809 of file Molecule.h.

Referenced by NamdState::loadStructure().

00809 { qmReplaceAll = newReplaceAll; } ;

void Molecule::setBFactorData ( molfile_atom_t *  atomarray  ) 

Definition at line 3074 of file Molecule.C.

Referenced by Molecule().

03074                                                       {
03075     bfactor = new float[numAtoms];
03076     for(int i=0; i<numAtoms; i++) {
03077         bfactor[i] = atomarray[i].bfactor;
03078     }
03079 }

void Molecule::setOccupancyData ( molfile_atom_t *  atomarray  ) 

Definition at line 3067 of file Molecule.C.

Referenced by Molecule().

03067                                                         {
03068     occupancy = new float[numAtoms];
03069     for(int i=0; i<numAtoms; i++) {
03070         occupancy[i] = atomarray[i].occupancy;
03071     }
03072 }


Friends And Related Function Documentation

friend class AngleElem [friend]

Definition at line 215 of file Molecule.h.

friend class AnisoElem [friend]

Definition at line 219 of file Molecule.h.

friend class BondElem [friend]

Definition at line 214 of file Molecule.h.

friend class CrosstermElem [friend]

Definition at line 220 of file Molecule.h.

friend class DihedralElem [friend]

Definition at line 216 of file Molecule.h.

friend class ExclElem [friend]

Definition at line 213 of file Molecule.h.

friend class GromacsPairElem [friend]

Definition at line 222 of file Molecule.h.

friend class ImproperElem [friend]

Definition at line 217 of file Molecule.h.

friend class TholeElem [friend]

Definition at line 218 of file Molecule.h.

friend class WorkDistrib [friend]

Definition at line 224 of file Molecule.h.


Member Data Documentation

Definition at line 561 of file Molecule.h.

Referenced by NamdState::loadStructure().

Definition at line 562 of file Molecule.h.

Referenced by NamdState::loadStructure().

Definition at line 563 of file Molecule.h.

Referenced by NamdState::loadStructure().

ConsTorqueParams* Molecule::consTorqueParams

Definition at line 608 of file Molecule.h.

Referenced by get_constorque_params(), receive_Molecule(), and send_Molecule().

Definition at line 676 of file Molecule.h.

Referenced by build_go_arrays(), and goInit().

Definition at line 677 of file Molecule.h.

Referenced by build_go_arrays(), and goInit().

Definition at line 668 of file Molecule.h.

Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().

Definition at line 670 of file Molecule.h.

Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().

Definition at line 672 of file Molecule.h.

Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().

Definition at line 669 of file Molecule.h.

Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().

Definition at line 671 of file Molecule.h.

Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().

GoValue Molecule::go_array[MAX_GO_CHAINS *MAX_GO_CHAINS]
int Molecule::go_indices[MAX_GO_CHAINS+1]

Definition at line 1572 of file Molecule.h.

Referenced by read_go_file(), receive_GoMolecule(), and send_GoMolecule().

Definition at line 645 of file Molecule.h.

Referenced by build_go_sigmas2(), receive_GoMolecule(), and send_GoMolecule().

Definition at line 646 of file Molecule.h.

Referenced by build_go_sigmas2(), get_go_force2(), receive_GoMolecule(), and send_GoMolecule().

Definition at line 644 of file Molecule.h.

Referenced by build_go_sigmas2(), receive_GoMolecule(), and send_GoMolecule().

Definition at line 642 of file Molecule.h.

Referenced by build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), and goInit().

Definition at line 637 of file Molecule.h.

Referenced by build_go_sigmas2(), get_go_force2(), receive_GoMolecule(), and send_GoMolecule().

Definition at line 647 of file Molecule.h.

Referenced by build_go_sigmas2(), get_go_force2(), receive_GoMolecule(), and send_GoMolecule().

Definition at line 648 of file Molecule.h.

Referenced by build_go_sigmas2(), get_go_force2(), receive_GoMolecule(), and send_GoMolecule().

Definition at line 673 of file Molecule.h.

Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().

Definition at line 661 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 631 of file Molecule.h.

Referenced by WorkDistrib::createAtomLists(), and outputCompressedFile().

Definition at line 666 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 667 of file Molecule.h.

Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().

Definition at line 657 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 658 of file Molecule.h.

Referenced by get_gro_force2(), receive_Molecule(), and send_Molecule().

Definition at line 459 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 460 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 1452 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 1452 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 567 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Number of anisotropic terms.

Definition at line 575 of file Molecule.h.

Referenced by AnisoElem::getMoleculePointers().

Definition at line 623 of file Molecule.h.

Referenced by Controller::compareChecksums().

Definition at line 654 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 622 of file Molecule.h.

Referenced by Controller::compareChecksums().

Definition at line 603 of file Molecule.h.

Referenced by NamdState::loadStructure(), receive_Molecule(), and send_Molecule().

Definition at line 587 of file Molecule.h.

Referenced by is_atom_constorqued(), receive_Molecule(), and send_Molecule().

Definition at line 566 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Number of Drude particles.

Definition at line 573 of file Molecule.h.

Referenced by NamdState::loadStructure(), and Controller::receivePressure().

Definition at line 590 of file Molecule.h.

Referenced by is_atom_exPressure(), receive_Molecule(), and send_Molecule().

Definition at line 600 of file Molecule.h.

Referenced by NamdState::loadStructure(), receive_Molecule(), and send_Molecule().

Definition at line 665 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 1573 of file Molecule.h.

Referenced by print_go_params(), read_go_file(), receive_GoMolecule(), and send_GoMolecule().

Definition at line 583 of file Molecule.h.

Referenced by build_gridforce_params(), receive_Molecule(), and send_Molecule().

Definition at line 653 of file Molecule.h.

Referenced by build_go_sigmas2(), and GromacsPairElem::getMoleculePointers().

Number of lone pairs.

Definition at line 572 of file Molecule.h.

Referenced by NamdState::loadStructure(), num_deg_freedom(), and Controller::receivePressure().

Number of lone pair host records in PSF.

Definition at line 576 of file Molecule.h.

Definition at line 585 of file Molecule.h.

Referenced by is_atom_movdragged(), receive_Molecule(), and send_Molecule().

Definition at line 627 of file Molecule.h.

Referenced by ParallelIOMgr::bcastHydroBasedCounter(), and NamdState::loadStructure().

Definition at line 629 of file Molecule.h.

Referenced by ParallelIOMgr::bcastHydroBasedCounter(), and NamdState::loadStructure().