Molecule Class Reference

#include <Molecule.h>

List of all members.

Public Member Functions

void compute_LJcorrection ()
BigReal getEnergyTailCorr (const BigReal)
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_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
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 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_virial_1
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
int numTotalExclusions
int numLonepairs
int numDrudeAtoms
int numTholes
int numAnisos
int numLphosts
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
int numCalcExclusions
int 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

Classes

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


Detailed Description

Definition at line 150 of file Molecule.h.


Constructor & Destructor Documentation

Molecule::Molecule ( SimParameters ,
Parameters param 
)

Definition at line 418 of file Molecule.C.

References initialize(), and param.

00419 {
00420   initialize(simParams,param);
00421 }

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

Definition at line 431 of file Molecule.C.

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

00432 {
00433   initialize(simParams,param);
00434 
00435 #ifdef MEM_OPT_VERSION
00436   if(simParams->useCompressedPsf)
00437       read_mol_signatures(filename, param, cfgList);
00438 #else
00439         read_psf_file(filename, param); 
00440  //LCPO
00441   if (simParams->LCPOOn)
00442     assignLCPOTypes( 0 );
00443 #endif      
00444  }

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

Definition at line 453 of file Molecule.C.

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

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

Molecule::Molecule ( SimParameters ,
Parameters ,
Ambertoppar  
)

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

Molecule::~Molecule (  ) 

Definition at line 547 of file Molecule.C.

References atomSigPool, and numAtoms.

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


Member Function Documentation

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

Definition at line 1005 of file Molecule.h.

References atoms.

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

01006   {
01007     #ifdef MEM_OPT_VERSION
01008     return atomChargePool[eachAtomCharge[anum]];
01009     #else
01010     return(atoms[anum].charge);
01011     #endif
01012   }

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

Definition at line 995 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().

00996   {
00997     #ifdef MEM_OPT_VERSION
00998     return atomMassPool[eachAtomMass[anum]];
00999     #else
01000     return(atoms[anum].mass);
01001     #endif
01002   }

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

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 1015 of file Molecule.h.

References atoms.

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

01016   {      
01017       return(atoms[anum].vdw_type);
01018   }

void Molecule::build_constant_forces ( char *   ) 

Referenced by NamdState::loadStructure().

void Molecule::build_constorque_params ( StringList ,
StringList ,
StringList ,
StringList ,
StringList ,
StringList ,
PDB ,
char *   
)

Referenced by NamdState::loadStructure().

void Molecule::build_constraint_params ( StringList ,
StringList ,
StringList ,
PDB ,
char *   
)

Referenced by NamdState::loadStructure().

void Molecule::build_exPressure_atoms ( StringList ,
StringList ,
PDB ,
char *   
)

Referenced by NamdState::loadStructure().

void Molecule::build_extra_bonds ( Parameters parameters,
StringList file 
)

Referenced by NamdState::loadStructure().

void Molecule::build_fep_flags ( StringList ,
StringList ,
PDB ,
char *  ,
const char *   
)

Referenced by NamdState::loadStructure().

void Molecule::build_fixed_atoms ( StringList ,
StringList ,
PDB ,
char *   
)

Referenced by NamdState::loadStructure().

void Molecule::build_go_arrays ( StringList ,
char *   
)

Definition at line 950 of file GoMolecule.C.

References PDB::atom(), atomChainTypes, atoms_1to4(), StringList::data, DebugM, endi(), energyNative, energyNonnative, get_go_cutoff(), 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  ) 

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 ,
char *   
)

Definition at line 577 of file GoMolecule.C.

References PDB::atom(), atomChainTypes, atoms_1to4(), StringList::data, DebugM, endi(), get_go_cutoff(), get_go_exp_b(), 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 ,
char *   
)

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_b(), 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, and pointerToGoEnd.

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 ,
StringList ,
StringList ,
StringList ,
PDB ,
char *   
)

Definition at line 6098 of file Molecule.C.

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

Referenced by NamdState::loadStructure().

06104 {
06105     PDB *kPDB;
06106     register int i;             //  Loop counters
06107     register int j;
06108     register int k;
06109 
06110     DebugM(3,  "Entered build_gridforce_params multi...\n");
06111 //     DebugM(3, "\tgridfrcfile = " << gridfrcfile->data << endi);
06112 //     DebugM(3, "\tgridfrccol = " << gridfrccol->data << endi);
06113     
06114     MGridforceParams* mgridParams = simParams->mgridforcelist.get_first();
06115     numGridforceGrids = 0;
06116     while (mgridParams != NULL) {
06117         numGridforceGrids++;
06118         mgridParams = mgridParams->next;
06119     }
06120     
06121     DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n");
06122     gridfrcIndexes = new int32*[numGridforceGrids];
06123     gridfrcParams = new GridforceParams*[numGridforceGrids];
06124     gridfrcGrid = new GridforceGrid*[numGridforceGrids];
06125     numGridforces = new int[numGridforceGrids];
06126     
06127     int grandTotalGrids = 0;    // including all subgrids
06128     
06129     mgridParams = simParams->mgridforcelist.get_first();
06130     for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
06131         int current_index=0;    //  Index into values used
06132         int kcol = 5;           //  Column to look for force constant in
06133         int qcol = 0;           //  Column for charge (default 0: use electric charge)
06134         Real kval = 0;          //  Force constant value retreived
06135         char filename[129];     //  PDB filename
06136         char potfilename[129];  //  Potential file name
06137         
06138         if (mgridParams == NULL) {
06139             NAMD_die("Problem with mgridParams!");
06140         }
06141         
06142         // Now load values from mgridforcelist object
06143         if (mgridParams->gridforceFile == NULL)
06144         {
06145             if ( ! initial_pdb ) NAMD_die("Initial PDB file unavailable, gridforceFile required.");
06146             kPDB = initial_pdb;
06147         }
06148         else
06149         {
06150             DebugM(4, "mgridParams->gridforceFile = " << mgridParams->gridforceFile << "\n" << endi);
06151             
06152             if ( (cwd == NULL) || (mgridParams->gridforceFile[0] == '/') )
06153             {
06154                 strcpy(filename, mgridParams->gridforceFile);
06155             }
06156             else
06157             {
06158                 strcpy(filename, cwd);
06159                 strcat(filename, mgridParams->gridforceFile);
06160             }
06161         
06162             kPDB = new PDB(filename);
06163             if ( kPDB == NULL )
06164             {
06165                 NAMD_die("Memory allocation failed in Molecule::build_gridforce_params");
06166             }
06167            
06168             if (kPDB->num_atoms() != numAtoms)
06169             {
06170                 NAMD_die("Number of atoms in grid force PDB doesn't match coordinate PDB");
06171             }
06172         }
06173 
06174         //  Get the column that the force constant is going to be in.  It
06175         //  can be in any of the 5 floating point fields in the PDB, according
06176         //  to what the user wants.  The allowable fields are X, Y, Z, O, or
06177         //  B which correspond to the 1st, 2nd, ... 5th floating point fields.
06178         //  The default is the 5th field, which is beta (temperature factor)
06179         if (mgridParams->gridforceCol == NULL)
06180         {
06181             kcol = 5;
06182         }
06183         else
06184         {
06185             if (strcasecmp(mgridParams->gridforceCol, "X") == 0)
06186             {
06187                 kcol=1;
06188             }
06189             else if (strcasecmp(mgridParams->gridforceCol, "Y") == 0)
06190             {
06191                 kcol=2;
06192             }
06193             else if (strcasecmp(mgridParams->gridforceCol, "Z") == 0)
06194             {
06195                 kcol=3;
06196             }
06197             else if (strcasecmp(mgridParams->gridforceCol, "O") == 0)
06198             {
06199                 kcol=4;
06200             }
06201             else if (strcasecmp(mgridParams->gridforceCol, "B") == 0)
06202             {
06203                 kcol=5;
06204             }
06205             else
06206             {
06207                 NAMD_die("gridforcecol must have value of X, Y, Z, O, or B");
06208             }
06209         }
06210     
06211         //  Get the column that the charge is going to be in.
06212         if (mgridParams->gridforceQcol == NULL)
06213         {
06214             qcol = 0;   // Default: don't read charge from file, use electric charge
06215         }
06216         else
06217         {
06218             if (strcasecmp(mgridParams->gridforceQcol, "X") == 0)
06219             {
06220                 qcol=1;
06221             }
06222             else if (strcasecmp(mgridParams->gridforceQcol, "Y") == 0)
06223             {
06224                 qcol=2;
06225             }
06226             else if (strcasecmp(mgridParams->gridforceQcol, "Z") == 0)
06227             {
06228                 qcol=3;
06229             }
06230             else if (strcasecmp(mgridParams->gridforceQcol, "O") == 0)
06231             {
06232                 qcol=4;
06233             }
06234             else if (strcasecmp(mgridParams->gridforceQcol, "B") == 0)
06235             {
06236                 qcol=5;
06237             }
06238             else
06239             {
06240                 NAMD_die("gridforcechargecol must have value of X, Y, Z, O, or B");
06241             }
06242         }
06243     
06244         if (kcol == qcol) {
06245             NAMD_die("gridforcecol and gridforcechargecol cannot have same value");
06246         }
06247 
06248     
06249         //  Allocate an array that will store an index into the constraint
06250         //  parameters for each atom.  If the atom is not constrained, its
06251         //  value will be set to -1 in this array.
06252         gridfrcIndexes[gridnum] = new int32[numAtoms];
06253        
06254         if (gridfrcIndexes[gridnum] == NULL)
06255         {
06256             NAMD_die("memory allocation failed in Molecule::build_gridforce_params()");
06257         }
06258         
06259         //  Loop through all the atoms and find out which ones are constrained
06260         for (i=0; i<numAtoms; i++)
06261         {
06262             //  Get the k value based on where we were told to find it
06263             switch (kcol)
06264             {
06265             case 1:
06266                 kval = (kPDB->atom(i))->xcoor();
06267                 break;
06268             case 2:
06269                 kval = (kPDB->atom(i))->ycoor();
06270                 break;
06271             case 3:
06272                 kval = (kPDB->atom(i))->zcoor();
06273                 break;
06274             case 4:
06275                 kval = (kPDB->atom(i))->occupancy();
06276                 break;
06277             case 5:
06278                 kval = (kPDB->atom(i))->temperaturefactor();
06279                 break;
06280             }
06281            
06282             if (kval > 0.0)
06283             {
06284                 //  This atom is constrained
06285                 gridfrcIndexes[gridnum][i] = current_index;
06286                 current_index++;
06287             }
06288             else
06289             {
06290                 //  This atom is not constrained
06291                 gridfrcIndexes[gridnum][i] = -1;
06292             }
06293         }
06294     
06295         if (current_index == 0)
06296         {
06297             //  Constraints were turned on, but there weren't really any constrained
06298             iout << iWARN << "NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" << endi;
06299         }
06300         else
06301         {
06302             //  Allocate an array to hold the constraint parameters
06303             gridfrcParams[gridnum] = new GridforceParams[current_index];
06304             if (gridfrcParams[gridnum] == NULL)
06305             {
06306                 NAMD_die("memory allocation failed in Molecule::build_gridforce_params");
06307             }
06308         }
06309     
06310         numGridforces[gridnum] = current_index;
06311 
06312         //  Loop through all the atoms and assign the parameters for those
06313         //  that are constrained
06314         for (i=0; i<numAtoms; i++)
06315         {
06316             if (gridfrcIndexes[gridnum][i] != -1)
06317             {
06318                 //  This atom has grid force, so get the k value again
06319                 switch (kcol)
06320                 {
06321                 case 1:
06322                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->xcoor();
06323                     break;
06324                 case 2:
06325                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->ycoor();
06326                     break;
06327                 case 3:
06328                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->zcoor();
06329                     break;
06330                 case 4:
06331                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->occupancy();
06332                     break;
06333                 case 5:
06334                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->temperaturefactor();
06335                     break;
06336                 }
06337             
06338                 //  Also get charge column
06339                 switch (qcol)
06340                 {
06341                 case 0:
06342 #ifdef MEM_OPT_VERSION
06343                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]];
06344 #else
06345                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge;
06346 #endif
06347                     break;
06348                 case 1:
06349                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->xcoor();
06350                     break;
06351                 case 2:
06352                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->ycoor();
06353                     break;
06354                 case 3:
06355                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->zcoor();
06356                     break;
06357                 case 4:
06358                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->occupancy();
06359                     break;
06360                 case 5:
06361                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->temperaturefactor();
06362                     break;
06363                 }
06364             }
06365         }
06366        
06367         //  If we had to create new PDB objects, delete them now
06368         if (mgridParams->gridforceFile != NULL)
06369         {
06370             delete kPDB;
06371         }
06372     
06373         //  Now we fill in our grid information
06374     
06375         // Open potential file
06376         if ( (cwd == NULL) || (mgridParams->gridforceVfile[0] == '/') )
06377         {
06378             strcpy(potfilename, mgridParams->gridforceVfile);
06379         }
06380         else
06381         {
06382             strcpy(potfilename, cwd);
06383             strcat(potfilename, mgridParams->gridforceVfile);
06384         }
06385     
06386 //        iout << iINFO << "Allocating grid " << gridnum
06387 //             << "\n" << endi;
06388         
06389         DebugM(3, "allocating GridforceGrid(" << gridnum << ")\n");
06390         gridfrcGrid[gridnum] = GridforceGrid::new_grid(gridnum, potfilename, simParams, mgridParams);
06391         
06392         grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
06393         DebugM(4, "grandTotalGrids = " << grandTotalGrids << "\n" << endi);
06394         
06395         // Finally, get next mgridParams pointer
06396         mgridParams = mgridParams->next;
06397     }
06398 }

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 
)

Referenced by NamdState::loadStructure().

void Molecule::build_movdrag_params ( StringList ,
StringList ,
StringList ,
PDB ,
char *   
)

Referenced by NamdState::loadStructure().

void Molecule::build_rotdrag_params ( StringList ,
StringList ,
StringList ,
StringList ,
StringList ,
StringList ,
PDB ,
char *   
)

Referenced by NamdState::loadStructure().

void Molecule::build_stirred_atoms ( StringList ,
StringList ,
PDB ,
char *   
)

Referenced by NamdState::loadStructure().

int Molecule::checkexcl ( int  atom1,
int  atom2 
) const

void Molecule::compute_LJcorrection (  ) 

Referenced by NamdState::loadStructure().

void Molecule::delete_alch_bonded ( void   ) 

Referenced by NamdState::loadStructure().

void Molecule::delete_qm_bonded (  ) 

Definition at line 1547 of file MoleculeQM.C.

References bond::atom1, bond::atom2, DebugM, SimParameters::extraBondsOn, iERROR(), iout, and NAMD_die().

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 991 of file Molecule.h.

Referenced by NamdState::loadStructure().

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

void Molecule::freeOccupancyData (  )  [inline]

Definition at line 987 of file Molecule.h.

Referenced by NamdState::loadStructure().

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

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

Definition at line 1062 of file Molecule.h.

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

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

Definition at line 1025 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

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

Definition at line 1104 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01105       { 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 1073 of file Molecule.h.

References atomTypePool, and NAMD_die().

01074   {
01075     if (atomNames == NULL)
01076     {
01077       NAMD_die("Tried to find atom type on node other than node 0");
01078     }
01079 
01080     #ifdef MEM_OPT_VERSION    
01081     return atomTypePool[atomNames[anum].atomtypeIdx];
01082     #else
01083     return(atomNames[anum].atomtype);
01084     #endif
01085   }

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

Definition at line 1022 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

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

Definition at line 1102 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01103       { return bondsByAtom[anum]; } 

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

Definition at line 980 of file Molecule.h.

Referenced by wrap_coor_int().

00980 { return cluster[anum]; }

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

Definition at line 981 of file Molecule.h.

Referenced by wrap_coor_int().

00981 { return clusterSize[anum]; }

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

Definition at line 1219 of file Molecule.h.

Referenced by ComputeRestraints::doForce().

01220   {
01221     k = consParams[consIndexes[atomnum]].k;
01222     refPos = consParams[consIndexes[atomnum]].refPos;
01223   }

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

Definition at line 1293 of file Molecule.h.

References consTorqueIndexes, and consTorqueParams.

Referenced by ComputeConsTorque::doForce().

01295   {
01296     v = consTorqueParams[consTorqueIndexes[atomnum]].v;
01297     a = consTorqueParams[consTorqueIndexes[atomnum]].a;
01298     p = consTorqueParams[consTorqueIndexes[atomnum]].p;
01299   }

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

Definition at line 1034 of file Molecule.h.

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

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

Definition at line 1110 of file Molecule.h.

01111       { return crosstermsByAtom[anum]; }  

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

Definition at line 854 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00854 {return cSMDcoffs;} ;

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

Definition at line 850 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00850 {return cSMDindex;} ;

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

Definition at line 849 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00849 { return cSMDindxLen;} ;

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

Definition at line 852 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00852 {return cSMDKs;} ;

int Molecule::get_cSMDnumInst (  )  [inline]

Definition at line 848 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00848 { return cSMDnumInst;} ;

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

Definition at line 851 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00851 {return cSMDpairs;} ;

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

Definition at line 853 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00853 {return cSMDVels;} ;

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

Definition at line 1031 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

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

Definition at line 1106 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01107       { return dihedralsByAtom[anum]; }

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

Definition at line 1059 of file Molecule.h.

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

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

Definition at line 1130 of file Molecule.h.

Referenced by dumpbench(), and SELF().

01130                                                                {      
01131       return &all_exclusions[anum];             
01132   }

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

Definition at line 1069 of file Molecule.h.

References exclusions.

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

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

Definition at line 1112 of file Molecule.h.

References exclusionsByAtom.

01113       { return exclusionsByAtom[anum]; }

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

Definition at line 1335 of file Molecule.h.

References if(), NAMD_die(), and simParams.

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

01336   {
01337     int typeSum = 0;
01338     for ( int i=0; i < order; ++i ) {
01339       typeSum += (fepAtomFlags[atomID[i]] == 2 ? -1 : fepAtomFlags[atomID[i]]);
01340     }
01341     // Increase the cutoff if scaling purely alchemical bonds. 
01342     // This really only applies when order = 2.
01343     if ( simParams->alchBondDecouple ) order++;
01344 
01345     if ( typeSum == 0 ) return 0; // Most interactions get caught here.
01346     else if ( 0 < typeSum && typeSum < order ) return 1;
01347     else if ( 0 > typeSum && typeSum > -order ) return 2;
01348 
01349     if ( simParams->alchBondDecouple ) {
01350       // Alchemify should always keep this from bombing, but just in case...
01351       NAMD_die("Unexpected alchemical bonded interaction!");
01352     }
01353   }

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

Definition at line 1302 of file Molecule.h.

Referenced by WorkDistrib::createAtomLists(), and ComputeHomeTuples< AnisoElem, aniso, aniso >::loadTuples().

01303   {
01304     return(fepAtomFlags[anum]);
01305   }

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

Definition at line 1114 of file Molecule.h.

Referenced by ComputeNonbondedCUDA::build_exclusions(), and SELF().

01115       { return fullExclusionsByAtom[anum]; }

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

Definition at line 1551 of file Molecule.h.

References go_val::cutoff, 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().

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

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

Definition at line 1560 of file Molecule.h.

References go_val::epsilon, and MAX_GO_CHAINS.

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

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

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

Definition at line 1554 of file Molecule.h.

References go_val::epsilonRep, and MAX_GO_CHAINS.

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

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

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

Definition at line 1563 of file Molecule.h.

References go_val::exp_a, and MAX_GO_CHAINS.

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

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

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

Definition at line 1566 of file Molecule.h.

References go_val::exp_b, and MAX_GO_CHAINS.

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

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

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

Definition at line 1569 of file Molecule.h.

References go_val::exp_rep, and MAX_GO_CHAINS.

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

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

BigReal Molecule::get_go_force ( BigReal  ,
int  ,
int  ,
BigReal ,
BigReal  
) 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(), goForce, 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  ,
BigReal  ,
BigReal  ,
int  ,
int  ,
BigReal ,
BigReal  
) 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  ,
int  ,
int  ,
BigReal ,
BigReal  
) 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, goForce, 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*(expRep/(r*r)) * 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 1557 of file Molecule.h.

References MAX_GO_CHAINS, and go_val::sigmaRep.

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

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

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

Definition at line 1233 of file Molecule.h.

References numGridforceGrids.

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

01234   {
01235       GridforceGrid *result = NULL;
01236       if (gridnum >= 0 && gridnum < numGridforceGrids) {
01237           result = gridfrcGrid[gridnum];
01238       }
01239       return result;
01240   }

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

Definition at line 1227 of file Molecule.h.

Referenced by ComputeGridForce::do_calc().

01228   {
01229       k = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].k;
01230       q = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].q;
01231   }

BigReal Molecule::get_gro_force ( BigReal  ,
BigReal  ,
BigReal  ,
int  ,
int   
) const

BigReal Molecule::get_gro_force2 ( BigReal  ,
BigReal  ,
BigReal  ,
int  ,
int  ,
BigReal ,
BigReal  
) 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 1028 of file Molecule.h.

Referenced by dumpbench().

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

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

Definition at line 1108 of file Molecule.h.

Referenced by dumpbench().

01109       { return impropersByAtom[anum]; }  

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

Definition at line 1038 of file Molecule.h.

01038                                        {
01039     // don't call unless simParams->drudeOn == TRUE
01040     // otherwise lphostIndexes array doesn't exist!
01041     int index = lphostIndexes[atomid];
01042     return (index != -1 ? &(lphosts[index]) : NULL);
01043   }

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

Definition at line 1116 of file Molecule.h.

Referenced by SELF().

01117       { return modExclusionsByAtom[anum]; }

int Molecule::get_mother_atom ( int   ) 

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

Definition at line 1278 of file Molecule.h.

Referenced by Sequencer::addMovDragToPosition().

01279   {
01280     v = movDragParams[movDragIndexes[atomnum]].v;
01281   }

Bool Molecule::get_noPC (  )  [inline]

Definition at line 836 of file Molecule.h.

Referenced by ComputeQM::initialize().

00836 { return qmNoPC; } ;

int Molecule::get_numQMAtoms (  )  [inline]

Definition at line 805 of file Molecule.h.

Referenced by ComputePme::doQMWork(), ComputeQM::initialize(), ComputeQMMgr::procQMRes(), HomePatch::qmSwapAtoms(), ComputeQMMgr::recvPartQM(), and ComputeQM::saveResults().

00805 {return qmNumQMAtoms; } ;

Real* Molecule::get_qmAtmChrg (  )  [inline]

Definition at line 802 of file Molecule.h.

Referenced by ComputePme::doQMWork(), ComputeQM::initialize(), HomePatch::qmSwapAtoms(), and ComputeQM::saveResults().

00802 {return qmAtmChrg; } ;

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

Definition at line 803 of file Molecule.h.

Referenced by ComputePme::doQMWork(), ComputeQM::initialize(), HomePatch::qmSwapAtoms(), and ComputeQM::saveResults().

00803 {return qmAtmIndx; } ;

Real Molecule::get_qmAtomGroup ( int  indx  )  const [inline]

Definition at line 800 of file Molecule.h.

00800 {return qmAtomGroup[indx]; } ;

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

Definition at line 799 of file Molecule.h.

Referenced by ComputePme::doQMWork(), ComputeQM::initialize(), HomePatch::qmSwapAtoms(), ComputeQM::saveResults(), and SELF().

00799 {return qmAtomGroup; } ;

Bool Molecule::get_qmcSMD (  )  [inline]

Definition at line 847 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00847 { return qmcSMD;} ;

int* Molecule::get_qmCustomPCIdxs (  )  [inline]

Definition at line 845 of file Molecule.h.

Referenced by ComputeQM::initialize().

00845 { return qmCustomPCIdxs; } ;

int* Molecule::get_qmCustPCSizes (  )  [inline]

Definition at line 844 of file Molecule.h.

Referenced by ComputeQM::initialize().

00844 { return qmCustPCSizes; } ;

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

Definition at line 814 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00814 { return qmDummyBondVal; } ;

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

Definition at line 819 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00819 { return qmDummyElement; } ;

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

Definition at line 806 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00806 {return qmElementArray;} ;

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

Definition at line 817 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00817 { return qmGrpBonds; } ;

Real* Molecule::get_qmGrpChrg (  )  [inline]

Definition at line 810 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00810 {return qmGrpChrg; } ;

Real* Molecule::get_qmGrpID (  )  [inline]

Definition at line 809 of file Molecule.h.

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

00809 { return qmGrpID; } ;

Real* Molecule::get_qmGrpMult (  )  [inline]

Definition at line 811 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00811 {return qmGrpMult; } ;

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

Definition at line 815 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00815 { return qmGrpNumBonds; } ;

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

Definition at line 808 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00808 {return qmGrpSizes; } ;

int Molecule::get_qmLSSFreq (  )  [inline]

Definition at line 827 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00827 { return qmLSSFreq; } ;

int* Molecule::get_qmLSSIdxs (  )  [inline]

Definition at line 825 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00825 { return qmLSSIdxs;} ;

Mass* Molecule::get_qmLSSMass (  )  [inline]

Definition at line 826 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00826 { return qmLSSMass; } ;

int* Molecule::get_qmLSSRefIDs (  )  [inline]

Definition at line 829 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00829 { return qmLSSRefIDs ; } ;

Mass* Molecule::get_qmLSSRefMass (  )  [inline]

Definition at line 831 of file Molecule.h.

00831 {return qmLSSRefMass; } ;

int* Molecule::get_qmLSSRefSize (  )  [inline]

Definition at line 830 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00830 { return qmLSSRefSize ; } ;

int Molecule::get_qmLSSResSize (  )  [inline]

Definition at line 828 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00828 { return qmLSSResidueSize; } ;

int* Molecule::get_qmLSSSize (  )  [inline]

Definition at line 824 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00824 { return qmLSSSize;} ;

int* Molecule::get_qmMeMMindx (  )  [inline]

Definition at line 838 of file Molecule.h.

Referenced by ComputeQM::initialize().

00838 { return qmMeMMindx; } ;

int Molecule::get_qmMeNumBonds (  )  [inline]

Definition at line 837 of file Molecule.h.

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

00837 { return qmMeNumBonds; } ;

Real* Molecule::get_qmMeQMGrp (  )  [inline]

Definition at line 839 of file Molecule.h.

Referenced by ComputeQM::initialize().

00839 { return qmMeQMGrp; } ;

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

Definition at line 816 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00816 { return qmMMBond; } ;

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

Definition at line 818 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00818 { return qmMMBondedIndx; } ;

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

Definition at line 821 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00821 { return qmMMChargeTarget; } ;

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

Definition at line 822 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00822 { return qmMMNumTargs; } ;

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

Definition at line 832 of file Molecule.h.

00832 { return qmClassicSolv;} ;

int Molecule::get_qmNumBonds (  )  [inline]

Definition at line 813 of file Molecule.h.

00813 { return qmNumBonds; } ;

int Molecule::get_qmNumGrps (  )  const [inline]

Definition at line 807 of file Molecule.h.

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

00807 {return qmNumGrps; } ;

int Molecule::get_qmPCFreq (  )  [inline]

Definition at line 841 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00841 { return qmPCFreq; } ;

Bool Molecule::get_qmReplaceAll (  )  [inline]

Definition at line 834 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00834 {return qmReplaceAll; } ;

int Molecule::get_qmTotCustPCs (  )  [inline]

Definition at line 843 of file Molecule.h.

Referenced by ComputeQM::initialize().

00843 { 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(), and GlobalMasterEasy::getNumAtoms().

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 1284 of file Molecule.h.

Referenced by Sequencer::addRotDragToPosition().

01286   {
01287     v = rotDragParams[rotDragIndexes[atomnum]].v;
01288     a = rotDragParams[rotDragIndexes[atomnum]].a;
01289     p = rotDragParams[rotDragIndexes[atomnum]].p;
01290   }

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

Definition at line 1259 of file Molecule.h.

01260   {
01261     refPos = stirParams[stirIndexes[atomnum]].refPos;
01262   }

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

Definition at line 1271 of file Molecule.h.

Referenced by ComputeStir::doForce().

01272   {
01273     return stirParams[stirIndexes[atomnum]].startTheta;
01274   }

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

Definition at line 1065 of file Molecule.h.

01065 {return acceptors;}

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

Definition at line 1048 of file Molecule.h.

Referenced by buildAngleData().

01048 {return angles;}

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

Definition at line 1047 of file Molecule.h.

Referenced by buildBondData().

01047 {return bonds;}

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

Definition at line 1051 of file Molecule.h.

Referenced by buildCrosstermData().

01051 {return crossterms;}

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

Definition at line 1050 of file Molecule.h.

Referenced by buildDihedralData().

01050 {return dihedrals;}

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

Definition at line 1064 of file Molecule.h.

01064 {return donors;}

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

Definition at line 1049 of file Molecule.h.

Referenced by buildImproperData().

01049 {return impropers;}

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

Definition at line 1055 of file Molecule.h.

01055 { return lphosts; }

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

Definition at line 469 of file Molecule.h.

References drude_constants::alpha.

Referenced by SELF().

00469                                     {
00470     return drudeConsts[i].alpha;
00471   }

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

Definition at line 478 of file Molecule.h.

Referenced by buildAtomData().

00478 { return atomNames; }

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

Definition at line 477 of file Molecule.h.

References atoms.

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

00477 { return atoms; }

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

Definition at line 481 of file Molecule.h.

Referenced by buildAtomData().

00481 { return atomSegResids; }

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

Definition at line 989 of file Molecule.h.

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

00989 { return (const float *)bfactor; }

BigReal Molecule::getEnergyTailCorr ( const   BigReal  ) 

Referenced by Controller::printEnergies().

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

Definition at line 465 of file Molecule.h.

Referenced by HomePatch::setLcpoType().

00465                                  {
00466     return lcpoParamType;
00467   }

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

Definition at line 985 of file Molecule.h.

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

00985 { return (const float *)occupancy; }

BigReal Molecule::getVirialTailCorr ( const   BigReal  ) 

Referenced by Controller::calcPressure().

Bool Molecule::go_restricted ( int  ,
int  ,
int   
)

Definition at line 525 of file GoMolecule.C.

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

Referenced by 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 1203 of file Molecule.h.

References consTorqueIndexes, FALSE, and numConsTorque.

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

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

Definition at line 1154 of file Molecule.h.

References FALSE, and numConstraints.

Referenced by ComputeRestraints::doForce().

01155   {
01156     if (numConstraints)
01157     {
01158       //  Check the index to see if it is constrained
01159       return(consIndexes[atomnum] != -1);
01160     }
01161     else
01162     {
01163       //  No constraints at all, so just return FALSE
01164       return(FALSE);
01165     }
01166   }

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

Definition at line 1397 of file Molecule.h.

Referenced by Sequencer::langevinPiston().

01398   {
01399     return (numExPressureAtoms && exPressureAtomFlags[atomnum]);
01400   }

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

Definition at line 1357 of file Molecule.h.

Referenced by WorkDistrib::createAtomLists().

01358   {
01359     return (numFixedAtoms && fixedAtomFlags[atomnum]);
01360   }

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

Definition at line 1138 of file Molecule.h.

References FALSE, and numGridforceGrids.

Referenced by ComputeGridForce::do_calc().

01139   {
01140       if (numGridforceGrids)
01141       {
01142           return(gridfrcIndexes[gridnum][atomnum] != -1);
01143       }
01144       else
01145       {
01146           return(FALSE);
01147       }
01148   }

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

Definition at line 1171 of file Molecule.h.

References FALSE, and numMovDrag.

Referenced by Sequencer::addMovDragToPosition().

01172   {
01173     if (numMovDrag)
01174     {
01175       //  Check the index to see if it is constrained
01176       return(movDragIndexes[atomnum] != -1);
01177     }
01178     else
01179     {
01180       //  No constraints at all, so just return FALSE
01181       return(FALSE);
01182     }
01183   }

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

Definition at line 1187 of file Molecule.h.

References FALSE, and numRotDrag.

Referenced by Sequencer::addRotDragToPosition().

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

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

Definition at line 1378 of file Molecule.h.

References FALSE.

Referenced by ComputeStir::doForce().

01379   {
01380     if (numStirredAtoms)
01381     {
01382       //  Check the index to see if it is constrained
01383       return(stirIndexes[atomnum] != -1);
01384     }
01385     else
01386     {
01387       //  No constraints at all, so just return FALSE
01388       return(FALSE);
01389     }
01390   }

Bool Molecule::is_drude ( int   ) 

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

Definition at line 1393 of file Molecule.h.

01394   {
01395     return (numFixedAtoms && (fixedAtomFlags[atomnum] == -1));
01396   }

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   ) 

Referenced by outputCompressedFile(), and wrap_coor_int().

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

Definition at line 1253 of file Molecule.h.

Referenced by WorkDistrib::createAtomLists().

01254   {
01255     return(langevinParams ? langevinParams[atomnum] : 0.);
01256   }

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

Definition at line 510 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().

00510                                                          {
00511     // local variables prefixed by s_
00512     int64_t s_NumDegFreedom = 3 * (int64_t) numAtoms;
00513     int64_t s_NumFixedAtoms = num_fixed_atoms();
00514     if (s_NumFixedAtoms) s_NumDegFreedom -= 3 * s_NumFixedAtoms;
00515     if (numLonepairs) s_NumDegFreedom -= 3 * numLonepairs;
00516     if ( ! (s_NumFixedAtoms || numConstraints
00517           || simParams->comMove || simParams->langevinOn) ) {
00518       s_NumDegFreedom -= 3;
00519     }
00520     if ( ! isInitialReport && simParams->pairInteractionOn) {
00521       //
00522       // DJH: a kludge?  We want to initially report system degrees of freedom
00523       //
00524       // this doesn't attempt to deal with fixed atoms or constraints
00525       s_NumDegFreedom = 3 * numFepInitial;
00526     }
00527     int s_NumFixedRigidBonds = 
00528       (simParams->fixedAtomsOn ? numFixedRigidBonds : 0);
00529     if (simParams->watmodel == WAT_TIP4) {
00530       // numLonepairs is subtracted here because all lonepairs have a rigid bond
00531       // to oxygen, but all of the LP degrees of freedom are dealt with above
00532       s_NumDegFreedom -= (numRigidBonds - s_NumFixedRigidBonds - numLonepairs);
00533     }
00534     else {
00535       // Non-polarized systems don't have LPs.
00536       // For Drude model, bonds that attach LPs are not counted as rigid;
00537       // LPs have already been subtracted from degrees of freedom.
00538       s_NumDegFreedom -= (numRigidBonds - s_NumFixedRigidBonds);
00539     }
00540     return s_NumDegFreedom;
00541   }

int Molecule::num_fixed_atoms (  )  const [inline]

Definition at line 484 of file Molecule.h.

References numFixedAtoms, and simParams.

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

00484                               {
00485     // local variables prefixed by s_
00486     int s_NumFixedAtoms = (simParams->fixedAtomsOn ? numFixedAtoms : 0);
00487     return s_NumFixedAtoms;  // value is "turned on" SimParameters
00488   }

int Molecule::num_fixed_groups (  )  const [inline]

Definition at line 490 of file Molecule.h.

References num_fixed_atoms(), and numFixedGroups.

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

00490                                {
00491     // local variables prefixed by s_
00492     int s_NumFixedAtoms = num_fixed_atoms();
00493     int s_NumFixedGroups = (s_NumFixedAtoms ? numFixedGroups : 0);
00494     return s_NumFixedGroups;  // value is "turned on" SimParameters
00495   }

int64_t Molecule::num_group_deg_freedom (  )  const [inline]

Definition at line 497 of file Molecule.h.

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

Referenced by Controller::receivePressure().

00497                                         {
00498     // local variables prefixed by s_
00499     int64_t s_NumGroupDegFreedom = 3 * (int64_t) numHydrogenGroups;
00500     int64_t s_NumFixedAtoms = num_fixed_atoms();
00501     int s_NumFixedGroups = num_fixed_groups();
00502     if (s_NumFixedGroups) s_NumGroupDegFreedom -= 3 * s_NumFixedGroups;
00503     if ( ! (s_NumFixedAtoms || numConstraints
00504           || simParams->comMove || simParams->langevinOn) ) {
00505       s_NumGroupDegFreedom -= 3;
00506     }
00507     return s_NumGroupDegFreedom;
00508   }

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

Definition at line 109 of file MoleculeQM.C.

References PDB::atom(), StringList::data, DebugM, endi(), ConfigList::find(), iERROR(), if(), iINFO(), iout, NAMD_die(), StringList::next, PDB::num_atoms(), numAtoms, PDBAtom::occupancy(), SimParameters::qmBondOn, SimParameters::qmColumn, SimParameters::qmLSSMode, QMLSSMODECOM, SimParameters::qmLSSOn, SimParameters::qmNoPC, SimParameters::qmPCSelFreq, eabffunc::split(), 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     
01467     if (simParams->qmCSMD) {
01468         qmcSMD = simParams->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  ) 

Definition at line 5116 of file Molecule.C.

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

Referenced by NamdState::loadStructure().

05117 {
05118 #ifdef MEM_OPT_VERSION
05119     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
05120 #else   
05121   register int i;
05122   Real sigma;
05123   Real epsilon;
05124   Real sigma14;
05125   Real epsilon14;
05126 
05127   DebugM(2,"ATOM LIST\n" \
05128       << "******************************************\n" \
05129                   << "NUM  NAME TYPE RES  MASS    CHARGE CHARGE   FEP-CHARGE"  \
05130       << "SIGMA   EPSILON SIGMA14 EPSILON14\n" \
05131         << endi);
05132 
05133   for (i=0; i<numAtoms; i++)
05134   {
05135     params->get_vdw_params(&sigma, &epsilon, &sigma14, &epsilon14, 
05136         atoms[i].vdw_type);
05137 
05138     DebugM(2,i+1 << " " << atomNames[i].atomname  \
05139               << " " << atomNames[i].atomtype << " " \
05140               << atomNames[i].resname  << " " << atoms[i].mass  \
05141         << " " << atoms[i].charge << " " << sigma \
05142         << " " << epsilon << " " << sigma14 \
05143         << " " << epsilon14 << "\n" \
05144         << endi);
05145   }
05146 #endif  
05147 }

void Molecule::print_bonds ( Parameters  ) 

Definition at line 5159 of file Molecule.C.

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

Referenced by NamdState::loadStructure().

05160 {
05161 #ifdef MEM_OPT_VERSION
05162     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
05163 #else   
05164   register int i;
05165   Real k;
05166   Real x0;
05167 
05168   DebugM(2,"BOND LIST\n" << "********************************\n" \
05169       << "ATOM1 ATOM2 TYPE1 TYPE2      k        x0" \
05170       << endi);
05171 
05172   for (i=0; i<numBonds; i++)
05173   {
05174     params->get_bond_params(&k, &x0, bonds[i].bond_type);
05175 
05176     DebugM(2,bonds[i].atom1+1 << " " \
05177        << bonds[i].atom2+1 << " "   \
05178        << atomNames[bonds[i].atom1].atomtype << " "  \
05179        << atomNames[bonds[i].atom2].atomtype << " " << k \
05180        << " " << x0 << endi);
05181   }
05182   
05183 #endif  
05184 }

void Molecule::print_exclusions (  ) 

Definition at line 5196 of file Molecule.C.

References DebugM, and endi().

Referenced by NamdState::loadStructure().

05197 {
05198 #ifdef MEM_OPT_VERSION
05199     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
05200 #else
05201   register int i;
05202 
05203   DebugM(2,"EXPLICIT EXCLUSION LIST\n" \
05204       << "********************************\n" \
05205             << "ATOM1 ATOM2 " \
05206       << endi);
05207 
05208   for (i=0; i<numExclusions; i++)
05209   {
05210     DebugM(2,exclusions[i].atom1+1 << "  " \
05211        << exclusions[i].atom2+1 << endi);
05212   }
05213 #endif
05214 }

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 1265 of file Molecule.h.

01266   {
01267     stirParams[stirIndexes[atomnum]].startTheta = theta;
01268   }

void Molecule::read_go_file ( char *   ) 

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  ) 

Definition at line 1744 of file GoMolecule.C.

References atomChainTypes, go_val::cutoff, 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, 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  ) 

Definition at line 5564 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(), 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::stirOn, SimParameters::tCoupleOn, and GridforceGrid::unpack_grid().

Referenced by Node::resendMolecule().

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

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

Referenced by Node::reloadCharges().

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

Definition at line 1403 of file Molecule.h.

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

01404   {
01405     return(rigidBondLengths[atomnum]);
01406   }

void Molecule::send_GoMolecule ( MOStream  ) 

Definition at line 1635 of file GoMolecule.C.

References atomChainTypes, go_val::cutoff, MOStream::end(), 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(), 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  ) 

Definition at line 5227 of file Molecule.C.

References SimParameters::alchOn, atomNamePool, atomSigPool, consForce, consForceIndexes, SimParameters::consForceOn, consTorqueIndexes, SimParameters::consTorqueOn, consTorqueParams, SimParameters::constraintsOn, DebugM, endi(), SimParameters::excludeFromPressure, SimParameters::fixedAtomsOn, gA, 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::stirOn, and SimParameters::tCoupleOn.

Referenced by Node::resendMolecule().

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

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

Definition at line 1242 of file Molecule.h.

References numGridforceGrids.

Referenced by Node::reloadGridforceGrid().

01243   {
01244       if (grid && gridnum >= 0 && gridnum < numGridforceGrids) {
01245           gridfrcGrid[gridnum] = grid;
01246           return 0;
01247       } else {
01248           return -1;
01249       }
01250   }

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

Definition at line 797 of file Molecule.h.

Referenced by NamdState::loadStructure().

00797 { qmReplaceAll = newReplaceAll; } ;

void Molecule::setBFactorData ( molfile_atom_t *  atomarray  ) 

Definition at line 3024 of file Molecule.C.

Referenced by Molecule().

03024                                                       {
03025     bfactor = new float[numAtoms];
03026     for(int i=0; i<numAtoms; i++) {
03027         bfactor[i] = atomarray[i].bfactor;
03028     }
03029 }

void Molecule::setOccupancyData ( molfile_atom_t *  atomarray  ) 

Definition at line 3017 of file Molecule.C.

Referenced by Molecule().

03017                                                         {
03018     occupancy = new float[numAtoms];
03019     for(int i=0; i<numAtoms; i++) {
03020         occupancy[i] = atomarray[i].occupancy;
03021     }
03022 }


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

int Molecule::alchDroppedAngles

Definition at line 550 of file Molecule.h.

Referenced by NamdState::loadStructure().

int Molecule::alchDroppedDihedrals

Definition at line 551 of file Molecule.h.

Referenced by NamdState::loadStructure().

int Molecule::alchDroppedImpropers

Definition at line 552 of file Molecule.h.

Referenced by NamdState::loadStructure().

int32* Molecule::atomChainTypes

Definition at line 623 of file Molecule.h.

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

Vector* Molecule::consForce

Definition at line 593 of file Molecule.h.

Referenced by ComputeConsForce::doForce(), receive_Molecule(), ComputeMgr::recvComputeConsForceMsg(), and send_Molecule().

int32* Molecule::consForceIndexes

Definition at line 592 of file Molecule.h.

Referenced by ComputeConsForce::doForce(), receive_Molecule(), ComputeMgr::recvComputeConsForceMsg(), and send_Molecule().

int32* Molecule::consTorqueIndexes

Definition at line 595 of file Molecule.h.

Referenced by ComputeConsTorque::doForce(), get_constorque_params(), is_atom_constorqued(), receive_Molecule(), and send_Molecule().

ConsTorqueParams* Molecule::consTorqueParams

Definition at line 596 of file Molecule.h.

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

BigReal Molecule::energyNative

Definition at line 664 of file Molecule.h.

Referenced by build_go_arrays(), and goInit().

BigReal Molecule::energyNonnative

Definition at line 665 of file Molecule.h.

Referenced by build_go_arrays(), and goInit().

Real* Molecule::gA

Definition at line 656 of file Molecule.h.

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

Real* Molecule::giSigma1

Definition at line 658 of file Molecule.h.

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

Real* Molecule::giSigma2

Definition at line 660 of file Molecule.h.

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

Real* Molecule::gMu1

Definition at line 657 of file Molecule.h.

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

Real* Molecule::gMu2

Definition at line 659 of file Molecule.h.

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

GoValue Molecule::go_array[MAX_GO_CHAINS *MAX_GO_CHAINS]

Definition at line 1534 of file Molecule.h.

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

int Molecule::go_indices[MAX_GO_CHAINS+1]

Definition at line 1535 of file Molecule.h.

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

Real* Molecule::goCoordinates

Definition at line 628 of file Molecule.h.

Referenced by build_go_arrays(), get_go_force_new(), goInit(), receive_GoMolecule(), and send_GoMolecule().

int* Molecule::goIndxLJA

Definition at line 633 of file Molecule.h.

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

int* Molecule::goIndxLJB

Definition at line 634 of file Molecule.h.

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

int Molecule::goNumLJPair

Definition at line 632 of file Molecule.h.

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

PDB* Molecule::goPDB

Definition at line 630 of file Molecule.h.

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

int32* Molecule::goResidIndices

Definition at line 625 of file Molecule.h.

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

int* Molecule::goResids

Definition at line 629 of file Molecule.h.

Referenced by build_go_arrays(), get_go_force_new(), goInit(), receive_GoMolecule(), and send_GoMolecule().

int32* Molecule::goSigmaIndices

Definition at line 624 of file Molecule.h.

Referenced by build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), get_go_force(), get_go_force_new(), goInit(), print_go_sigmas(), receive_GoMolecule(), and send_GoMolecule().

double* Molecule::goSigmaPairA

Definition at line 635 of file Molecule.h.

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

double* Molecule::goSigmaPairB

Definition at line 636 of file Molecule.h.

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

Real* Molecule::goSigmas

Definition at line 626 of file Molecule.h.

Referenced by build_go_sigmas(), get_go_force(), goInit(), print_go_sigmas(), receive_GoMolecule(), and send_GoMolecule().

bool* Molecule::goWithinCutoff

Definition at line 627 of file Molecule.h.

Referenced by build_go_sigmas(), get_go_force(), goInit(), receive_GoMolecule(), and send_GoMolecule().

Real* Molecule::gRepulsive

Definition at line 661 of file Molecule.h.

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

int* Molecule::gromacsPair_type

Definition at line 649 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

HydrogenGroup Molecule::hydrogenGroup

Definition at line 619 of file Molecule.h.

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

int* Molecule::indxGaussA

Definition at line 654 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

int* Molecule::indxGaussB

Definition at line 655 of file Molecule.h.

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

int* Molecule::indxLJA

Definition at line 645 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

int* Molecule::indxLJB

Definition at line 646 of file Molecule.h.

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

int Molecule::is_drude_psf

Definition at line 448 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

int Molecule::is_lonepairs_psf

Definition at line 449 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

int Molecule::isBFactorValid

Definition at line 1415 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

int Molecule::isOccupancyValid

Definition at line 1415 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

int Molecule::maxHydrogenGroupSize

Definition at line 580 of file Molecule.h.

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

int Molecule::maxMigrationGroupSize

Definition at line 582 of file Molecule.h.

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

int Molecule::numAcceptors

Definition at line 556 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

int Molecule::numAngles

Definition at line 547 of file Molecule.h.

Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ParallelIOMgr::bcastMolInfo(), buildAngleData(), dumpbench(), AngleElem::getMoleculePointers(), NamdState::loadStructure(), Molecule(), and ParallelIOMgr::recvMolInfo().

int Molecule::numAnisos

Definition at line 564 of file Molecule.h.

Referenced by AnisoElem::getMoleculePointers().

int Molecule::numAtoms

Definition at line 543 of file Molecule.h.

Referenced by GlobalMasterFreeEnergy::addForce(), GlobalMasterEasy::addForce(), WorkDistrib::assignNodeToPatch(), ParallelIOMgr::bcastHydroBasedCounter(), build12Excls(), build13Excls(), build14Excls(), ComputeNonbondedCUDA::build_exclusions(), build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), buildAtomData(), buildBondData(), buildExclusions(), colvarproxy_namd::calculate(), Controller::compareChecksums(), ComputeMgr::createComputes(), dumpbench(), GlobalMasterFreeEnergy::getMass(), GlobalMasterEasy::getMass(), GlobalMasterSymmetry::GlobalMasterSymmetry(), GlobalMasterTMD::GlobalMasterTMD(), colvarproxy_namd::init_atom_group(), integrateAllAtomSigs(), loadMolInfo(), NamdState::loadStructure(), Molecule(), num_deg_freedom(), outputCompressedFile(), WorkDistrib::patchMapInit(), prepare_qm(), print_go_sigmas(), receive_GoMolecule(), Controller::receivePressure(), ParallelIOMgr::recvAtomsCntPerPatch(), ComputeMgr::recvComputeConsForceMsg(), ComputeMsmSerialMgr::recvCoord(), ComputeGBISserMgr::recvCoord(), ComputeFmmSerialMgr::recvCoord(), ComputeExtMgr::recvCoord(), ComputeQMMgr::recvPartQM(), Node::reloadCharges(), GlobalMasterFreeEnergy::requestAtom(), GlobalMasterEasy::requestAtom(), Node::resendMolecule(), Node::resendMolecule2(), SELF(), send_GoMolecule(), Node::startup(), Tcl_centerOfMass(), Tcl_centerOfNumber(), Tcl_loadCoords(), Tcl_radiusOfGyration(), wrap_coor_int(), and ~Molecule().

int Molecule::numBonds

Definition at line 546 of file Molecule.h.

Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ParallelIOMgr::bcastMolInfo(), buildBondData(), dumpbench(), BondElem::getMoleculePointers(), NamdState::loadStructure(), Molecule(), and ParallelIOMgr::recvMolInfo().

int Molecule::numCalcAngles

Definition at line 602 of file Molecule.h.

Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), dumpbench(), receive_Molecule(), ParallelIOMgr::recvMolInfo(), and send_Molecule().

int Molecule::numCalcAnisos

Definition at line 611 of file Molecule.h.

Referenced by Controller::compareChecksums().

int Molecule::numCalcBonds

Definition at line 601 of file Molecule.h.

Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), dumpbench(), receive_Molecule(), ParallelIOMgr::recvMolInfo(), and send_Molecule().

int Molecule::numCalcCrossterms

Definition at line 605 of file Molecule.h.

Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), receive_Molecule(), ParallelIOMgr::recvMolInfo(), and send_Molecule().

int Molecule::numCalcDihedrals

Definition at line 603 of file Molecule.h.

Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), dumpbench(), receive_Molecule(), ParallelIOMgr::recvMolInfo(), and send_Molecule().

int Molecule::numCalcExclusions

Definition at line 606 of file Molecule.h.

Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), dumpbench(), receive_Molecule(), ParallelIOMgr::recvMolInfo(), and send_Molecule().

int Molecule::numCalcFullExclusions

Definition at line 607 of file Molecule.h.

Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), receive_Molecule(), ParallelIOMgr::recvMolInfo(), and send_Molecule().

int Molecule::numCalcImpropers

Definition at line 604 of file Molecule.h.

Referenced by ParallelIOMgr::bcastMolInfo(), Controller::compareChecksums(), dumpbench(), receive_Molecule(), ParallelIOMgr::recvMolInfo(), and send_Molecule().

int Molecule::numCalcLJPair

Definition at line 642 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

int Molecule::numCalcTholes

Definition at line 610 of file Molecule.h.

Referenced by Controller::compareChecksums().

int Molecule::numConsForce

Definition at line 591 of file Molecule.h.

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

int Molecule::numConsTorque

Definition at line 575 of file Molecule.h.

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

int Molecule::numConstraints

Definition at line 568 of file Molecule.h.

Referenced by is_atom_constrained(), NamdState::loadStructure(), num_deg_freedom(), num_group_deg_freedom(), receive_Molecule(), Controller::receivePressure(), and send_Molecule().

int Molecule::numCrossterms

Definition at line 554 of file Molecule.h.

Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ParallelIOMgr::bcastMolInfo(), buildCrosstermData(), CrosstermElem::getMoleculePointers(), NamdState::loadStructure(), Molecule(), and ParallelIOMgr::recvMolInfo().

int Molecule::numDihedrals

Definition at line 548 of file Molecule.h.

Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ParallelIOMgr::bcastMolInfo(), buildDihedralData(), dumpbench(), DihedralElem::getMoleculePointers(), NamdState::loadStructure(), Molecule(), and ParallelIOMgr::recvMolInfo().

int Molecule::numDonors

Definition at line 555 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

int Molecule::numDrudeAtoms

Definition at line 562 of file Molecule.h.

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

int Molecule::numExclusions

Definition at line 557 of file Molecule.h.

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

int Molecule::numExPressureAtoms

Definition at line 578 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

int Molecule::numFepFinal

Definition at line 588 of file Molecule.h.

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

int Molecule::numFepInitial

Definition at line 587 of file Molecule.h.

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

int Molecule::numFixedAtoms

Definition at line 576 of file Molecule.h.

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

int Molecule::numFixedGroups

Definition at line 583 of file Molecule.h.

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

int Molecule::numFixedRigidBonds

Definition at line 585 of file Molecule.h.

Referenced by ParallelIOMgr::bcastHydroBasedCounter(), NamdState::loadStructure(), num_deg_freedom(), receive_Molecule(), Controller::receivePressure(), ParallelIOMgr::recvHydroBasedCounter(), and send_Molecule().

int Molecule::numGaussPair

Definition at line 653 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

int Molecule::numGoAtoms

Definition at line 622 of file Molecule.h.

Referenced by build_go_arrays(), build_go_sigmas(), build_go_sigmas2(), get_go_force(), goInit(), print_go_sigmas(), receive_GoMolecule(), and send_GoMolecule().

int Molecule::NumGoChains

Definition at line 1536 of file Molecule.h.

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

int Molecule::numGridforceGrids

Definition at line 570 of file Molecule.h.

Referenced by build_gridforce_params(), ComputeGridForce::doForce(), get_gridfrc_grid(), is_atom_gridforced(), NamdState::loadStructure(), receive_Molecule(), send_Molecule(), and set_gridfrc_grid().

int* Molecule::numGridforces

Definition at line 571 of file Molecule.h.

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

int Molecule::numHydrogenGroups

Definition at line 579 of file Molecule.h.

Referenced by ParallelIOMgr::bcastHydroBasedCounter(), NamdState::loadStructure(), num_group_deg_freedom(), outputCompressedFile(), receive_Molecule(), Controller::receivePressure(), and send_Molecule().

int Molecule::numImpropers

Definition at line 553 of file Molecule.h.

Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ParallelIOMgr::bcastMolInfo(), buildImproperData(), dumpbench(), ImproperElem::getMoleculePointers(), NamdState::loadStructure(), Molecule(), and ParallelIOMgr::recvMolInfo().

int Molecule::numLJPair

Definition at line 641 of file Molecule.h.

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

int Molecule::numLonepairs

Definition at line 561 of file Molecule.h.

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

int Molecule::numLphosts

Definition at line 565 of file Molecule.h.

int Molecule::numMigrationGroups

Definition at line 581 of file Molecule.h.

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

int Molecule::numMovDrag

Definition at line 573 of file Molecule.h.

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

int Molecule::numMultipleDihedrals

Definition at line 615 of file Molecule.h.

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

int Molecule::numMultipleImpropers

Definition at line 617 of file Molecule.h.

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

int Molecule::numPair

Definition at line 640 of file Molecule.h.

int Molecule::numRealBonds

Definition at line 545 of file Molecule.h.

Referenced by buildBondData(), and Molecule().

int Molecule::numRigidBonds

Definition at line 584 of file Molecule.h.

Referenced by ParallelIOMgr::bcastHydroBasedCounter(), ParallelIOMgr::bcastMolInfo(), ParallelIOMgr::integrateMigratedAtoms(), NamdState::loadStructure(), num_deg_freedom(), Controller::receivePressure(), and ParallelIOMgr::recvMolInfo().

int Molecule::numRotDrag

Definition at line 574 of file Molecule.h.

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

int Molecule::numStirredAtoms

Definition at line 577 of file Molecule.h.

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

int Molecule::numTholes

Definition at line 563 of file Molecule.h.

Referenced by TholeElem::getMoleculePointers().

int Molecule::numTotalExclusions

Definition at line 558 of file Molecule.h.

Referenced by ParallelIOMgr::bcastMolInfo(), and ParallelIOMgr::recvMolInfo().

Real* Molecule::pairC12

Definition at line 648 of file Molecule.h.

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

Real* Molecule::pairC6

Definition at line 647 of file Molecule.h.

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

int* Molecule::pointerToGaussBeg

Definition at line 651 of file Molecule.h.

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

int* Molecule::pointerToGaussEnd

Definition at line 652 of file Molecule.h.

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

int* Molecule::pointerToGoBeg

Definition at line 637 of file Molecule.h.

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

int* Molecule::pointerToGoEnd

Definition at line 638 of file Molecule.h.

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

int* Molecule::pointerToLJBeg

Definition at line 643 of file Molecule.h.

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

int* Molecule::pointerToLJEnd

Definition at line 644 of file Molecule.h.

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

Real Molecule::r_ohc

Definition at line 454 of file Molecule.h.

Real Molecule::r_om

Definition at line 453 of file Molecule.h.

int Molecule::suspiciousAlchBonds

Definition at line 549 of file Molecule.h.

Referenced by NamdState::loadStructure().

BigReal Molecule::tail_corr_dUdl_1

Definition at line 459 of file Molecule.h.

Referenced by Controller::printEnergies().

BigReal Molecule::tail_corr_ener

Definition at line 457 of file Molecule.h.

Referenced by Controller::printEnergies(), and Controller::rescaleaccelMD().

BigReal Molecule::tail_corr_virial

Definition at line 458 of file Molecule.h.

BigReal Molecule::tail_corr_virial_1

Definition at line 460 of file Molecule.h.


The documentation for this class was generated from the following files:
Generated on Sat Sep 23 01:17:20 2017 for NAMD by  doxygen 1.4.7