Molecule Class Reference

#include <Molecule.h>

List of all members.

Classes

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

Public Member Functions

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

Public Attributes

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

Friends

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

Detailed Description

Definition at line 150 of file Molecule.h.


Constructor & Destructor Documentation

Molecule::Molecule ( SimParameters simParams,
Parameters param 
)

Definition at line 429 of file Molecule.C.

References initialize().

00430 {
00431   initialize(simParams,param);
00432 }

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

Definition at line 442 of file Molecule.C.

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

00443 {
00444   initialize(simParams,param);
00445 
00446 #ifdef MEM_OPT_VERSION
00447   if(simParams->useCompressedPsf)
00448       read_mol_signatures(filename, param, cfgList);
00449 #else
00450         read_psf_file(filename, param); 
00451  //LCPO
00452   if (simParams->LCPOOn)
00453     assignLCPOTypes( 0 );
00454 #endif      
00455  }

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

Definition at line 464 of file Molecule.C.

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

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

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

Definition at line 558 of file Molecule.C.

References alch_unpert_angles, alch_unpert_bonds, alch_unpert_dihedrals, atomSigPool, numAtoms, ss_index, and ss_vdw_type.

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


Member Function Documentation

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

Definition at line 1048 of file Molecule.h.

References atoms.

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

01049   {
01050     #ifdef MEM_OPT_VERSION
01051     return atomChargePool[eachAtomCharge[anum]];
01052     #else
01053     return(atoms[anum].charge);
01054     #endif
01055   }

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

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

01039   {
01040     #ifdef MEM_OPT_VERSION
01041     return atomMassPool[eachAtomMass[anum]];
01042     #else
01043     return(atoms[anum].mass);
01044     #endif
01045   }

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

Definition at line 1537 of file GoMolecule.C.

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

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

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

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

Definition at line 1058 of file Molecule.h.

References atoms.

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

01059   {      
01060       return(atoms[anum].vdw_type);
01061   }

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

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

Definition at line 950 of file GoMolecule.C.

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

Referenced by NamdState::loadStructure().

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

void Molecule::build_go_params ( StringList g  ) 

Definition at line 79 of file GoMolecule.C.

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

Referenced by NamdState::loadStructure().

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

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

Definition at line 577 of file GoMolecule.C.

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

Referenced by NamdState::loadStructure().

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

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

Definition at line 747 of file GoMolecule.C.

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

Referenced by NamdState::loadStructure().

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

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

Definition at line 6373 of file Molecule.C.

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

Referenced by NamdState::loadStructure().

06379 {
06380     PDB *kPDB;
06381     register int i;             //  Loop counters
06382     register int j;
06383     register int k;
06384 
06385     DebugM(3,  "Entered build_gridforce_params multi...\n");
06386 //     DebugM(3, "\tgridfrcfile = " << gridfrcfile->data << endi);
06387 //     DebugM(3, "\tgridfrccol = " << gridfrccol->data << endi);
06388     
06389     MGridforceParams* mgridParams = simParams->mgridforcelist.get_first();
06390     numGridforceGrids = 0;
06391     while (mgridParams != NULL) {
06392         numGridforceGrids++;
06393         mgridParams = mgridParams->next;
06394     }
06395     
06396     DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n");
06397     gridfrcIndexes = new int32*[numGridforceGrids];
06398     gridfrcParams = new GridforceParams*[numGridforceGrids];
06399     gridfrcGrid = new GridforceGrid*[numGridforceGrids];
06400     numGridforces = new int[numGridforceGrids];
06401     
06402     int grandTotalGrids = 0;    // including all subgrids
06403     
06404     mgridParams = simParams->mgridforcelist.get_first();
06405     for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
06406         int current_index=0;    //  Index into values used
06407         int kcol = 5;           //  Column to look for force constant in
06408         int qcol = 0;           //  Column for charge (default 0: use electric charge)
06409         Real kval = 0;          //  Force constant value retreived
06410         char filename[129];     //  PDB filename
06411         char potfilename[129];  //  Potential file name
06412         
06413         if (mgridParams == NULL) {
06414             NAMD_die("Problem with mgridParams!");
06415         }
06416         
06417         // Now load values from mgridforcelist object
06418         if (mgridParams->gridforceFile == NULL)
06419         {
06420             if ( ! initial_pdb ) NAMD_die("Initial PDB file unavailable, gridforceFile required.");
06421             kPDB = initial_pdb;
06422         }
06423         else
06424         {
06425             DebugM(4, "mgridParams->gridforceFile = " << mgridParams->gridforceFile << "\n" << endi);
06426             
06427             if ( (cwd == NULL) || (mgridParams->gridforceFile[0] == '/') )
06428             {
06429                 strcpy(filename, mgridParams->gridforceFile);
06430             }
06431             else
06432             {
06433                 strcpy(filename, cwd);
06434                 strcat(filename, mgridParams->gridforceFile);
06435             }
06436         
06437             kPDB = new PDB(filename);
06438             if ( kPDB == NULL )
06439             {
06440                 NAMD_die("Memory allocation failed in Molecule::build_gridforce_params");
06441             }
06442            
06443             if (kPDB->num_atoms() != numAtoms)
06444             {
06445                 NAMD_die("Number of atoms in grid force PDB doesn't match coordinate PDB");
06446             }
06447         }
06448 
06449         //  Get the column that the force constant is going to be in.  It
06450         //  can be in any of the 5 floating point fields in the PDB, according
06451         //  to what the user wants.  The allowable fields are X, Y, Z, O, or
06452         //  B which correspond to the 1st, 2nd, ... 5th floating point fields.
06453         //  The default is the 5th field, which is beta (temperature factor)
06454         if (mgridParams->gridforceCol == NULL)
06455         {
06456             kcol = 5;
06457         }
06458         else
06459         {
06460             if (strcasecmp(mgridParams->gridforceCol, "X") == 0)
06461             {
06462                 kcol=1;
06463             }
06464             else if (strcasecmp(mgridParams->gridforceCol, "Y") == 0)
06465             {
06466                 kcol=2;
06467             }
06468             else if (strcasecmp(mgridParams->gridforceCol, "Z") == 0)
06469             {
06470                 kcol=3;
06471             }
06472             else if (strcasecmp(mgridParams->gridforceCol, "O") == 0)
06473             {
06474                 kcol=4;
06475             }
06476             else if (strcasecmp(mgridParams->gridforceCol, "B") == 0)
06477             {
06478                 kcol=5;
06479             }
06480             else
06481             {
06482                 NAMD_die("gridforcecol must have value of X, Y, Z, O, or B");
06483             }
06484         }
06485     
06486         //  Get the column that the charge is going to be in.
06487         if (mgridParams->gridforceQcol == NULL)
06488         {
06489             qcol = 0;   // Default: don't read charge from file, use electric charge
06490         }
06491         else
06492         {
06493             if (strcasecmp(mgridParams->gridforceQcol, "X") == 0)
06494             {
06495                 qcol=1;
06496             }
06497             else if (strcasecmp(mgridParams->gridforceQcol, "Y") == 0)
06498             {
06499                 qcol=2;
06500             }
06501             else if (strcasecmp(mgridParams->gridforceQcol, "Z") == 0)
06502             {
06503                 qcol=3;
06504             }
06505             else if (strcasecmp(mgridParams->gridforceQcol, "O") == 0)
06506             {
06507                 qcol=4;
06508             }
06509             else if (strcasecmp(mgridParams->gridforceQcol, "B") == 0)
06510             {
06511                 qcol=5;
06512             }
06513             else
06514             {
06515                 NAMD_die("gridforcechargecol must have value of X, Y, Z, O, or B");
06516             }
06517         }
06518     
06519         if (kcol == qcol) {
06520             NAMD_die("gridforcecol and gridforcechargecol cannot have same value");
06521         }
06522 
06523     
06524         //  Allocate an array that will store an index into the constraint
06525         //  parameters for each atom.  If the atom is not constrained, its
06526         //  value will be set to -1 in this array.
06527         gridfrcIndexes[gridnum] = new int32[numAtoms];
06528        
06529         if (gridfrcIndexes[gridnum] == NULL)
06530         {
06531             NAMD_die("memory allocation failed in Molecule::build_gridforce_params()");
06532         }
06533         
06534         //  Loop through all the atoms and find out which ones are constrained
06535         for (i=0; i<numAtoms; i++)
06536         {
06537             //  Get the k value based on where we were told to find it
06538             switch (kcol)
06539             {
06540             case 1:
06541                 kval = (kPDB->atom(i))->xcoor();
06542                 break;
06543             case 2:
06544                 kval = (kPDB->atom(i))->ycoor();
06545                 break;
06546             case 3:
06547                 kval = (kPDB->atom(i))->zcoor();
06548                 break;
06549             case 4:
06550                 kval = (kPDB->atom(i))->occupancy();
06551                 break;
06552             case 5:
06553                 kval = (kPDB->atom(i))->temperaturefactor();
06554                 break;
06555             }
06556            
06557             if (kval > 0.0)
06558             {
06559                 //  This atom is constrained
06560                 gridfrcIndexes[gridnum][i] = current_index;
06561                 current_index++;
06562             }
06563             else
06564             {
06565                 //  This atom is not constrained
06566                 gridfrcIndexes[gridnum][i] = -1;
06567             }
06568         }
06569     
06570         if (current_index == 0)
06571         {
06572             //  Constraints were turned on, but there weren't really any constrained
06573             iout << iWARN << "NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" << endi;
06574         }
06575         else
06576         {
06577             //  Allocate an array to hold the constraint parameters
06578             gridfrcParams[gridnum] = new GridforceParams[current_index];
06579             if (gridfrcParams[gridnum] == NULL)
06580             {
06581                 NAMD_die("memory allocation failed in Molecule::build_gridforce_params");
06582             }
06583         }
06584     
06585         numGridforces[gridnum] = current_index;
06586 
06587         //  Loop through all the atoms and assign the parameters for those
06588         //  that are constrained
06589         for (i=0; i<numAtoms; i++)
06590         {
06591             if (gridfrcIndexes[gridnum][i] != -1)
06592             {
06593                 //  This atom has grid force, so get the k value again
06594                 switch (kcol)
06595                 {
06596                 case 1:
06597                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->xcoor();
06598                     break;
06599                 case 2:
06600                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->ycoor();
06601                     break;
06602                 case 3:
06603                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->zcoor();
06604                     break;
06605                 case 4:
06606                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->occupancy();
06607                     break;
06608                 case 5:
06609                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->temperaturefactor();
06610                     break;
06611                 }
06612             
06613                 //  Also get charge column
06614                 switch (qcol)
06615                 {
06616                 case 0:
06617 #ifdef MEM_OPT_VERSION
06618                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]];
06619 #else
06620                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge;
06621 #endif
06622                     break;
06623                 case 1:
06624                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->xcoor();
06625                     break;
06626                 case 2:
06627                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->ycoor();
06628                     break;
06629                 case 3:
06630                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->zcoor();
06631                     break;
06632                 case 4:
06633                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->occupancy();
06634                     break;
06635                 case 5:
06636                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->temperaturefactor();
06637                     break;
06638                 }
06639             }
06640         }
06641        
06642         //  If we had to create new PDB objects, delete them now
06643         if (mgridParams->gridforceFile != NULL)
06644         {
06645             delete kPDB;
06646         }
06647     
06648         //  Now we fill in our grid information
06649     
06650         // Open potential file
06651         if ( (cwd == NULL) || (mgridParams->gridforceVfile[0] == '/') )
06652         {
06653             strcpy(potfilename, mgridParams->gridforceVfile);
06654         }
06655         else
06656         {
06657             strcpy(potfilename, cwd);
06658             strcat(potfilename, mgridParams->gridforceVfile);
06659         }
06660     
06661 //        iout << iINFO << "Allocating grid " << gridnum
06662 //             << "\n" << endi;
06663         
06664         DebugM(3, "allocating GridforceGrid(" << gridnum << ")\n");
06665         gridfrcGrid[gridnum] = GridforceGrid::new_grid(gridnum, potfilename, simParams, mgridParams);
06666         
06667         grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
06668         DebugM(4, "grandTotalGrids = " << grandTotalGrids << "\n" << endi);
06669         
06670         // Finally, get next mgridParams pointer
06671         mgridParams = mgridParams->next;
06672     }
06673 }

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

Build the flags needed for solute scaling.

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

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

Referenced by NamdState::loadStructure().

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

Definition at line 1547 of file MoleculeQM.C.

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

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

void Molecule::freeBFactorData (  )  [inline]

Definition at line 1034 of file Molecule.h.

Referenced by NamdState::loadStructure().

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

void Molecule::freeOccupancyData (  )  [inline]

Definition at line 1030 of file Molecule.h.

Referenced by NamdState::loadStructure().

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

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

Definition at line 1105 of file Molecule.h.

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

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

Definition at line 1068 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

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

Definition at line 1147 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01148       { 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 1116 of file Molecule.h.

References atomTypePool, and NAMD_die().

01117   {
01118     if (atomNames == NULL)
01119     {
01120       NAMD_die("Tried to find atom type on node other than node 0");
01121     }
01122 
01123     #ifdef MEM_OPT_VERSION    
01124     return atomTypePool[atomNames[anum].atomtypeIdx];
01125     #else
01126     return(atomNames[anum].atomtype);
01127     #endif
01128   }

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

Definition at line 1065 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

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

Definition at line 1145 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01146       { return bondsByAtom[anum]; } 

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

Definition at line 1023 of file Molecule.h.

Referenced by wrap_coor_int().

01023 { return cluster[anum]; }

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

Definition at line 1024 of file Molecule.h.

Referenced by wrap_coor_int().

01024 { return clusterSize[anum]; }

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

Definition at line 1262 of file Molecule.h.

Referenced by ComputeRestraints::doForce().

01263   {
01264     k = consParams[consIndexes[atomnum]].k;
01265     refPos = consParams[consIndexes[atomnum]].refPos;
01266   }

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

Definition at line 1336 of file Molecule.h.

References consTorqueIndexes, and consTorqueParams.

Referenced by ComputeConsTorque::doForce().

01338   {
01339     v = consTorqueParams[consTorqueIndexes[atomnum]].v;
01340     a = consTorqueParams[consTorqueIndexes[atomnum]].a;
01341     p = consTorqueParams[consTorqueIndexes[atomnum]].p;
01342   }

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

Definition at line 1077 of file Molecule.h.

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

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

Definition at line 1153 of file Molecule.h.

01154       { return crosstermsByAtom[anum]; }  

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

Definition at line 874 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00874 {return cSMDcoffs;} ;

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

Definition at line 870 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00870 {return cSMDindex;} ;

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

Definition at line 869 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00869 { return cSMDindxLen;} ;

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

Definition at line 872 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00872 {return cSMDKs;} ;

int Molecule::get_cSMDnumInst (  )  [inline]

Definition at line 868 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00868 { return cSMDnumInst;} ;

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

Definition at line 871 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00871 {return cSMDpairs;} ;

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

Definition at line 873 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00873 {return cSMDVels;} ;

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

Definition at line 1074 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

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

Definition at line 1149 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01150       { return dihedralsByAtom[anum]; }

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

Definition at line 1102 of file Molecule.h.

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

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

Definition at line 1173 of file Molecule.h.

Referenced by dumpbench().

01173                                                                {      
01174       return &all_exclusions[anum];             
01175   }

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

Definition at line 1112 of file Molecule.h.

References exclusions.

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

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

Definition at line 1155 of file Molecule.h.

References exclusionsByAtom.

01156       { return exclusionsByAtom[anum]; }

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

Definition at line 1385 of file Molecule.h.

References NAMD_die(), and simParams.

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

01386   {
01387     int typeSum = 0;
01388     for ( int i=0; i < order; ++i ) {
01389       typeSum += (fepAtomFlags[atomID[i]] == 2 ? -1 : fepAtomFlags[atomID[i]]);
01390     }
01391     // Increase the cutoff if scaling purely alchemical bonds. 
01392     // This really only applies when order = 2.
01393     if ( simParams->alchBondDecouple ) order++;
01394 
01395     // The majority of interactions are type 0, so catch those first.
01396     if ( typeSum == 0 || abs(typeSum) == order ) return 0;
01397     else if ( 0 < typeSum && typeSum < order ) return 1;
01398     else if ( -order < typeSum && typeSum < 0 ) return 2;
01399 
01400     // Alchemify should always keep this from bombing, but just in case...
01401     NAMD_die("Unexpected alchemical bonded interaction!");
01402     return 0;
01403   }

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

Definition at line 1157 of file Molecule.h.

Referenced by ComputeNonbondedCUDA::build_exclusions().

01158       { return fullExclusionsByAtom[anum]; }

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

Definition at line 1601 of file Molecule.h.

References go_val::cutoff, go_array, and MAX_GO_CHAINS.

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

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

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

Definition at line 1610 of file Molecule.h.

References go_val::epsilon, go_array, and MAX_GO_CHAINS.

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

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

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

Definition at line 1604 of file Molecule.h.

References go_val::epsilonRep, go_array, and MAX_GO_CHAINS.

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

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

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

Definition at line 1613 of file Molecule.h.

References go_val::exp_a, go_array, and MAX_GO_CHAINS.

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

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

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

Definition at line 1616 of file Molecule.h.

References go_val::exp_b, go_array, and MAX_GO_CHAINS.

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

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

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

Definition at line 1619 of file Molecule.h.

References go_val::exp_rep, go_array, and MAX_GO_CHAINS.

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

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

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

Definition at line 1260 of file GoMolecule.C.

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

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

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

Definition at line 1456 of file GoMolecule.C.

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

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

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

Definition at line 1334 of file GoMolecule.C.

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

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

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

Definition at line 1607 of file Molecule.h.

References go_array, MAX_GO_CHAINS, and go_val::sigmaRep.

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

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

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

Definition at line 1276 of file Molecule.h.

References numGridforceGrids.

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

01277   {
01278       GridforceGrid *result = NULL;
01279       if (gridnum >= 0 && gridnum < numGridforceGrids) {
01280           result = gridfrcGrid[gridnum];
01281       }
01282       return result;
01283   }

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

Definition at line 1270 of file Molecule.h.

Referenced by ComputeGridForce::do_calc().

01271   {
01272       k = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].k;
01273       q = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].q;
01274   }

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

Definition at line 1173 of file GoMolecule.C.

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

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

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

Definition at line 1071 of file Molecule.h.

Referenced by dumpbench().

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

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

Definition at line 1151 of file Molecule.h.

Referenced by dumpbench().

01152       { return impropersByAtom[anum]; }  

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

Definition at line 1081 of file Molecule.h.

01081                                        {
01082     // don't call unless simParams->drudeOn == TRUE
01083     // otherwise lphostIndexes array doesn't exist!
01084     int index = lphostIndexes[atomid];
01085     return (index != -1 ? &(lphosts[index]) : NULL);
01086   }

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

Definition at line 1159 of file Molecule.h.

01160       { return modExclusionsByAtom[anum]; }

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

Definition at line 1321 of file Molecule.h.

Referenced by Sequencer::addMovDragToPosition().

01322   {
01323     v = movDragParams[movDragIndexes[atomnum]].v;
01324   }

Bool Molecule::get_noPC (  )  [inline]

Definition at line 856 of file Molecule.h.

Referenced by ComputeQM::initialize().

00856 { return qmNoPC; } ;

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

Definition at line 820 of file Molecule.h.

00820 {return qmAtomGroup[indx]; } ;

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

Definition at line 867 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00867 { return qmcSMD;} ;

int* Molecule::get_qmCustomPCIdxs (  )  [inline]

Definition at line 865 of file Molecule.h.

Referenced by ComputeQM::initialize().

00865 { return qmCustomPCIdxs; } ;

int* Molecule::get_qmCustPCSizes (  )  [inline]

Definition at line 864 of file Molecule.h.

Referenced by ComputeQM::initialize().

00864 { return qmCustPCSizes; } ;

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

Definition at line 834 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00834 { return qmDummyBondVal; } ;

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

Definition at line 839 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00839 { return qmDummyElement; } ;

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

Definition at line 826 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00826 {return qmElementArray;} ;

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

Definition at line 837 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00837 { return qmGrpBonds; } ;

Real* Molecule::get_qmGrpChrg (  )  [inline]

Definition at line 830 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00830 {return qmGrpChrg; } ;

Real* Molecule::get_qmGrpID (  )  [inline]

Definition at line 829 of file Molecule.h.

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

00829 { return qmGrpID; } ;

Real* Molecule::get_qmGrpMult (  )  [inline]

Definition at line 831 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00831 {return qmGrpMult; } ;

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

Definition at line 835 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00835 { return qmGrpNumBonds; } ;

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

Definition at line 828 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00828 {return qmGrpSizes; } ;

int Molecule::get_qmLSSFreq (  )  [inline]

Definition at line 847 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00847 { return qmLSSFreq; } ;

int* Molecule::get_qmLSSIdxs (  )  [inline]

Definition at line 845 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00845 { return qmLSSIdxs;} ;

Mass* Molecule::get_qmLSSMass (  )  [inline]

Definition at line 846 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00846 { return qmLSSMass; } ;

int* Molecule::get_qmLSSRefIDs (  )  [inline]

Definition at line 849 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00849 { return qmLSSRefIDs ; } ;

Mass* Molecule::get_qmLSSRefMass (  )  [inline]

Definition at line 851 of file Molecule.h.

00851 {return qmLSSRefMass; } ;

int* Molecule::get_qmLSSRefSize (  )  [inline]

Definition at line 850 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00850 { return qmLSSRefSize ; } ;

int Molecule::get_qmLSSResSize (  )  [inline]

Definition at line 848 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00848 { return qmLSSResidueSize; } ;

int* Molecule::get_qmLSSSize (  )  [inline]

Definition at line 844 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00844 { return qmLSSSize;} ;

int* Molecule::get_qmMeMMindx (  )  [inline]

Definition at line 858 of file Molecule.h.

Referenced by ComputeQM::initialize().

00858 { return qmMeMMindx; } ;

int Molecule::get_qmMeNumBonds (  )  [inline]

Definition at line 857 of file Molecule.h.

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

00857 { return qmMeNumBonds; } ;

Real* Molecule::get_qmMeQMGrp (  )  [inline]

Definition at line 859 of file Molecule.h.

Referenced by ComputeQM::initialize().

00859 { return qmMeQMGrp; } ;

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

Definition at line 836 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00836 { return qmMMBond; } ;

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

Definition at line 838 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00838 { return qmMMBondedIndx; } ;

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

Definition at line 841 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00841 { return qmMMChargeTarget; } ;

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

Definition at line 842 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00842 { return qmMMNumTargs; } ;

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

Definition at line 852 of file Molecule.h.

00852 { return qmClassicSolv;} ;

int Molecule::get_qmNumBonds (  )  [inline]

Definition at line 833 of file Molecule.h.

00833 { return qmNumBonds; } ;

int Molecule::get_qmNumGrps (  )  const [inline]

Definition at line 827 of file Molecule.h.

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

00827 {return qmNumGrps; } ;

int Molecule::get_qmPCFreq (  )  [inline]

Definition at line 861 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00861 { return qmPCFreq; } ;

Bool Molecule::get_qmReplaceAll (  )  [inline]

Definition at line 854 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00854 {return qmReplaceAll; } ;

int Molecule::get_qmTotCustPCs (  )  [inline]

Definition at line 863 of file Molecule.h.

Referenced by ComputeQM::initialize().

00863 { return qmTotCustPCs; } ;

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

Definition at line 144 of file Molecule.C.

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

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

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

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

Definition at line 1327 of file Molecule.h.

Referenced by Sequencer::addRotDragToPosition().

01329   {
01330     v = rotDragParams[rotDragIndexes[atomnum]].v;
01331     a = rotDragParams[rotDragIndexes[atomnum]].a;
01332     p = rotDragParams[rotDragIndexes[atomnum]].p;
01333   }

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

Definition at line 1351 of file Molecule.h.

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

01352         {
01353                 return(ssAtomFlags[anum]);
01354         }

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

Definition at line 1302 of file Molecule.h.

01303   {
01304     refPos = stirParams[stirIndexes[atomnum]].refPos;
01305   }

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

Definition at line 1314 of file Molecule.h.

Referenced by ComputeStir::doForce().

01315   {
01316     return stirParams[stirIndexes[atomnum]].startTheta;
01317   }

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

Definition at line 1108 of file Molecule.h.

01108 {return acceptors;}

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

Definition at line 1091 of file Molecule.h.

Referenced by buildAngleData().

01091 {return angles;}

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

Definition at line 1090 of file Molecule.h.

Referenced by buildBondData().

01090 {return bonds;}

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

Definition at line 1094 of file Molecule.h.

Referenced by buildCrosstermData().

01094 {return crossterms;}

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

Definition at line 1093 of file Molecule.h.

Referenced by buildDihedralData().

01093 {return dihedrals;}

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

Definition at line 1107 of file Molecule.h.

01107 {return donors;}

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

Definition at line 1092 of file Molecule.h.

Referenced by buildImproperData().

01092 {return impropers;}

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

Definition at line 1098 of file Molecule.h.

01098 { return lphosts; }

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

Definition at line 482 of file Molecule.h.

References drude_constants::alpha.

00482                                     {
00483     return drudeConsts[i].alpha;
00484   }

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

Definition at line 491 of file Molecule.h.

Referenced by buildAtomData().

00491 { return atomNames; }

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

Definition at line 490 of file Molecule.h.

References atoms.

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

00490 { return atoms; }

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

Definition at line 494 of file Molecule.h.

Referenced by buildAtomData().

00494 { return atomSegResids; }

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

Definition at line 1032 of file Molecule.h.

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

01032 { return (const float *)bfactor; }

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

Definition at line 478 of file Molecule.h.

Referenced by HomePatch::setLcpoType().

00478                                  {
00479     return lcpoParamType;
00480   }

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

Definition at line 1028 of file Molecule.h.

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

01028 { return (const float *)occupancy; }

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

Definition at line 525 of file GoMolecule.C.

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

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

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

void Molecule::goInit (  ) 

Definition at line 54 of file GoMolecule.C.

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

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

void Molecule::initialize (  ) 

Referenced by Molecule().

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

Definition at line 1246 of file Molecule.h.

References consTorqueIndexes, FALSE, and numConsTorque.

01247   {
01248     if (numConsTorque)
01249     {
01250       //  Check the index to see if it is constrained
01251       return(consTorqueIndexes[atomnum] != -1);
01252     }
01253     else
01254     {
01255       //  No constraints at all, so just return FALSE
01256       return(FALSE);
01257     }
01258   }

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

Definition at line 1197 of file Molecule.h.

References FALSE, and numConstraints.

Referenced by ComputeRestraints::doForce().

01198   {
01199     if (numConstraints)
01200     {
01201       //  Check the index to see if it is constrained
01202       return(consIndexes[atomnum] != -1);
01203     }
01204     else
01205     {
01206       //  No constraints at all, so just return FALSE
01207       return(FALSE);
01208     }
01209   }

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

Definition at line 1447 of file Molecule.h.

References numExPressureAtoms.

Referenced by Sequencer::langevinPiston().

01448   {
01449     return (numExPressureAtoms && exPressureAtomFlags[atomnum]);
01450   }

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

Definition at line 1407 of file Molecule.h.

References numFixedAtoms.

Referenced by WorkDistrib::createAtomLists().

01408   {
01409     return (numFixedAtoms && fixedAtomFlags[atomnum]);
01410   }

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

Definition at line 1181 of file Molecule.h.

References FALSE, and numGridforceGrids.

Referenced by ComputeGridForce::do_calc().

01182   {
01183       if (numGridforceGrids)
01184       {
01185           return(gridfrcIndexes[gridnum][atomnum] != -1);
01186       }
01187       else
01188       {
01189           return(FALSE);
01190       }
01191   }

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

Definition at line 1214 of file Molecule.h.

References FALSE, and numMovDrag.

Referenced by Sequencer::addMovDragToPosition().

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

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

Definition at line 1230 of file Molecule.h.

References FALSE, and numRotDrag.

Referenced by Sequencer::addRotDragToPosition().

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

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

Definition at line 1428 of file Molecule.h.

References FALSE, and numStirredAtoms.

Referenced by ComputeStir::doForce().

01429   {
01430     if (numStirredAtoms)
01431     {
01432       //  Check the index to see if it is constrained
01433       return(stirIndexes[atomnum] != -1);
01434     }
01435     else
01436     {
01437       //  No constraints at all, so just return FALSE
01438       return(FALSE);
01439     }
01440   }

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

Definition at line 1443 of file Molecule.h.

References numFixedAtoms.

01444   {
01445     return (numFixedAtoms && (fixedAtomFlags[atomnum] == -1));
01446   }

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

Definition at line 1296 of file Molecule.h.

Referenced by WorkDistrib::createAtomLists().

01297   {
01298     return(langevinParams ? langevinParams[atomnum] : 0.);
01299   }

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

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

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

int Molecule::num_fixed_atoms (  )  const [inline]

Definition at line 497 of file Molecule.h.

References numFixedAtoms, and simParams.

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

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

int Molecule::num_fixed_groups (  )  const [inline]

Definition at line 503 of file Molecule.h.

References num_fixed_atoms(), and numFixedGroups.

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

00503                                {
00504     // local variables prefixed by s_
00505     int s_NumFixedAtoms = num_fixed_atoms();
00506     int s_NumFixedGroups = (s_NumFixedAtoms ? numFixedGroups : 0);
00507     return s_NumFixedGroups;  // value is "turned on" SimParameters
00508   }

int64_t Molecule::num_group_deg_freedom (  )  const [inline]

Definition at line 510 of file Molecule.h.

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

Referenced by Controller::receivePressure().

00510                                         {
00511     // local variables prefixed by s_
00512     int64_t s_NumGroupDegFreedom = 3 * (int64_t) numHydrogenGroups;
00513     int64_t s_NumFixedAtoms = num_fixed_atoms();
00514     int s_NumFixedGroups = num_fixed_groups();
00515     if (s_NumFixedGroups) s_NumGroupDegFreedom -= 3 * s_NumFixedGroups;
00516     if ( ! (s_NumFixedAtoms || numConstraints
00517           || simParams->comMove || simParams->langevinOn) ) {
00518       s_NumGroupDegFreedom -= 3;
00519     }
00520     return s_NumGroupDegFreedom;
00521   }

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

Prepares Live Solvent Selection

Data gathering from PDB to find QM atom and bond info

Multiplicity of each QM group

Charge of each QM group

Populate arrays that are used throughout the the calculations.

Overides Link Atom element with user selection.

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

Live Solvent Selection

Custom Point Charge selection

Topology preparation

Definition at line 109 of file MoleculeQM.C.

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

Referenced by NamdState::loadStructure().

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

void Molecule::print_atoms ( Parameters params  ) 

Definition at line 5337 of file Molecule.C.

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

Referenced by NamdState::loadStructure().

05338 {
05339 #ifdef MEM_OPT_VERSION
05340     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
05341 #else   
05342   register int i;
05343   Real sigma;
05344   Real epsilon;
05345   Real sigma14;
05346   Real epsilon14;
05347 
05348   DebugM(2,"ATOM LIST\n" \
05349       << "******************************************\n" \
05350                   << "NUM  NAME TYPE RES  MASS    CHARGE CHARGE   FEP-CHARGE"  \
05351       << "SIGMA   EPSILON SIGMA14 EPSILON14\n" \
05352         << endi);
05353 
05354   for (i=0; i<numAtoms; i++)
05355   {
05356     params->get_vdw_params(&sigma, &epsilon, &sigma14, &epsilon14, 
05357         atoms[i].vdw_type);
05358 
05359     DebugM(2,i+1 << " " << atomNames[i].atomname  \
05360               << " " << atomNames[i].atomtype << " " \
05361               << atomNames[i].resname  << " " << atoms[i].mass  \
05362         << " " << atoms[i].charge << " " << sigma \
05363         << " " << epsilon << " " << sigma14 \
05364         << " " << epsilon14 << "\n" \
05365         << endi);
05366   }
05367 #endif  
05368 }

void Molecule::print_bonds ( Parameters params  ) 

Definition at line 5380 of file Molecule.C.

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

Referenced by NamdState::loadStructure().

05381 {
05382 #ifdef MEM_OPT_VERSION
05383     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
05384 #else   
05385   register int i;
05386   Real k;
05387   Real x0;
05388 
05389   DebugM(2,"BOND LIST\n" << "********************************\n" \
05390       << "ATOM1 ATOM2 TYPE1 TYPE2      k        x0" \
05391       << endi);
05392 
05393   for (i=0; i<numBonds; i++)
05394   {
05395     params->get_bond_params(&k, &x0, bonds[i].bond_type);
05396 
05397     DebugM(2,bonds[i].atom1+1 << " " \
05398        << bonds[i].atom2+1 << " "   \
05399        << atomNames[bonds[i].atom1].atomtype << " "  \
05400        << atomNames[bonds[i].atom2].atomtype << " " << k \
05401        << " " << x0 << endi);
05402   }
05403   
05404 #endif  
05405 }

void Molecule::print_exclusions (  ) 

Definition at line 5417 of file Molecule.C.

References DebugM, and endi().

Referenced by NamdState::loadStructure().

05418 {
05419 #ifdef MEM_OPT_VERSION
05420     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
05421 #else
05422   register int i;
05423 
05424   DebugM(2,"EXPLICIT EXCLUSION LIST\n" \
05425       << "********************************\n" \
05426             << "ATOM1 ATOM2 " \
05427       << endi);
05428 
05429   for (i=0; i<numExclusions; i++)
05430   {
05431     DebugM(2,exclusions[i].atom1+1 << "  " \
05432        << exclusions[i].atom2+1 << endi);
05433   }
05434 #endif
05435 }

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

01309   {
01310     stirParams[stirIndexes[atomnum]].startTheta = theta;
01311   }

void Molecule::read_alch_unpert_angles ( FILE *  fd  ) 

Definition at line 1783 of file Molecule.C.

References alch_unpert_angles, angle::atom1, angle::atom2, angle::atom3, NAMD_die(), NAMD_read_int(), and num_alch_unpert_Angles.

01783                                                {
01784   int atom_nums[3];  //  Atom numbers for the three atoms
01785   register int j;    //  Loop counter
01786   int num_read=0;    //  Number of angles read so far
01787 
01788   alch_unpert_angles=new Angle[num_alch_unpert_Angles];
01789 
01790   if (alch_unpert_angles == NULL) {
01791     NAMD_die("memory allocation failed in Molecule::read_alch_unpert_angles");
01792   }
01793 
01794   while (num_read < num_alch_unpert_Angles) {
01795     for (j=0; j<3; j++) {
01796       atom_nums[j]=NAMD_read_int(fd, "ANGLES")-1;
01797 
01798       if (atom_nums[j] >= numAtoms) {
01799         char err_msg[128];
01800         sprintf(err_msg, "ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN ALCH UNPERT PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
01801         NAMD_die(err_msg);
01802       }
01803     }
01804 
01805     alch_unpert_angles[num_read].atom1=atom_nums[0];
01806     alch_unpert_angles[num_read].atom2=atom_nums[1];
01807     alch_unpert_angles[num_read].atom3=atom_nums[2];
01808 
01809     ++num_read;
01810   }
01811   return;
01812 }

void Molecule::read_alch_unpert_bonds ( FILE *  fd  ) 

Definition at line 1660 of file Molecule.C.

References alch_unpert_bonds, bond::atom1, bond::atom2, NAMD_die(), NAMD_read_int(), and num_alch_unpert_Bonds.

01660                                               {
01661   int atom_nums[2];  // Atom indexes for the bonded atoms
01662   register int j;      // Loop counter
01663   int num_read=0;    // Number of bonds read so far
01664 
01665   alch_unpert_bonds=new Bond[num_alch_unpert_Bonds];
01666 
01667   if (alch_unpert_bonds == NULL) {
01668     NAMD_die("memory allocations failed in Molecule::read_alch_unpert_bonds");
01669   }
01670 
01671   while (num_read < num_alch_unpert_Bonds) {
01672     for (j=0; j<2; j++) {
01673       atom_nums[j]=NAMD_read_int(fd, "BONDS")-1;
01674 
01675       if (atom_nums[j] >= numAtoms) {
01676         char err_msg[128];
01677 
01678         sprintf(err_msg, "BOND INDEX %d GREATER THAN NATOM %d IN BOND # %d IN ALCH PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
01679         NAMD_die(err_msg);
01680       }
01681     }
01682 
01683     Bond *b = &(alch_unpert_bonds[num_read]);
01684     b->atom1=atom_nums[0];
01685     b->atom2=atom_nums[1];
01686 
01687     ++num_read;
01688   }
01689   return;
01690 }

void Molecule::read_alch_unpert_dihedrals ( FILE *  fd  ) 

Definition at line 1931 of file Molecule.C.

References alch_unpert_dihedrals, dihedral::atom1, dihedral::atom2, dihedral::atom3, dihedral::atom4, NAMD_die(), NAMD_read_int(), and num_alch_unpert_Dihedrals.

01931                                                   {
01932   int atom_nums[4];  // The 4 atom indexes
01933   int num_read=0;    // number of dihedrals read so far
01934   register int j;    // loop counter
01935 
01936   alch_unpert_dihedrals = new Dihedral[num_alch_unpert_Dihedrals];
01937 
01938   if (alch_unpert_dihedrals == NULL) {
01939     NAMD_die("memory allocation failed in Molecule::read_alch_unpert_dihedrals");
01940   }
01941 
01942   while (num_read < num_alch_unpert_Dihedrals) {
01943     for (j=0; j<4; j++) {
01944       atom_nums[j]=NAMD_read_int(fd, "DIHEDRALS")-1;
01945 
01946       if (atom_nums[j] >= numAtoms) {
01947         char err_msg[128];
01948 
01949         sprintf(err_msg, "DIHEDRALS INDEX %d GREATER THAN NATOM %d IN DIHEDRALS # %d IN ALCH UNPERT PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
01950         NAMD_die(err_msg);
01951       }
01952     }
01953 
01954     alch_unpert_dihedrals[num_read].atom1=atom_nums[0];
01955     alch_unpert_dihedrals[num_read].atom2=atom_nums[1];
01956     alch_unpert_dihedrals[num_read].atom3=atom_nums[2];
01957     alch_unpert_dihedrals[num_read].atom4=atom_nums[3];
01958 
01959     num_read++;
01960   }
01961   return;
01962 }

void Molecule::read_go_file ( char *  fname  ) 

Definition at line 113 of file GoMolecule.C.

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

Referenced by build_go_params().

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

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

Definition at line 1744 of file GoMolecule.C.

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

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

void Molecule::receive_Molecule ( MIStream msg  ) 

Definition at line 5806 of file Molecule.C.

References alch_unpert_angles, alch_unpert_bonds, alch_unpert_dihedrals, SimParameters::alchOn, atomNamePool, atomSigPool, consForce, consForceIndexes, SimParameters::consForceOn, consTorqueIndexes, SimParameters::consTorqueOn, consTorqueParams, SimParameters::constraintsOn, DebugM, SimParameters::excludeFromPressure, SimParameters::fixedAtomsOn, gA, MIStream::get(), Parameters::get_num_vdw_params(), GridforceGrid::get_total_grids(), ObjectArena< Type >::getNewArray(), giSigma1, giSigma2, gMu1, gMu2, SimParameters::goGroPair, gRepulsive, gromacsPair_type, indxGaussA, indxGaussB, indxLJA, indxLJB, is_drude_psf, is_lonepairs_psf, isBFactorValid, isOccupancyValid, SimParameters::langevinOn, SimParameters::LCPOOn, SimParameters::lesOn, maxHydrogenGroupSize, maxMigrationGroupSize, SimParameters::mgridforceOn, SimParameters::movDragOn, num_alch_unpert_Angles, num_alch_unpert_Bonds, num_alch_unpert_Dihedrals, 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, numZeroMassAtoms, pairC12, pairC6, SimParameters::pairInteractionOn, pointerToGaussBeg, pointerToGaussEnd, pointerToLJBeg, pointerToLJEnd, SimParameters::qmForcesOn, SimParameters::rotDragOn, SimParameters::sdScaling, SimParameters::soluteScalingOn, ss_index, ss_num_vdw_params, ss_vdw_type, SimParameters::stirOn, SimParameters::tCoupleOn, and GridforceGrid::unpack_grid().

Referenced by Node::resendMolecule().

05806                                             {
05807   //  Get the atom information
05808   msg->get(numAtoms);
05809 
05810 #ifdef MEM_OPT_VERSION
05811 //in the memory optimized version, only the atom signatures are recved
05812 //from the master Node. --Chao Mei
05813 
05814   msg->get(massPoolSize);
05815   if(atomMassPool) delete [] atomMassPool;
05816   atomMassPool = new Real[massPoolSize];
05817   msg->get(massPoolSize, atomMassPool);
05818 
05819   msg->get(chargePoolSize);
05820   if(atomChargePool) delete [] atomChargePool;
05821   atomChargePool = new Real[chargePoolSize];
05822   msg->get(chargePoolSize, atomChargePool);
05823 
05824   //get atoms' signatures
05825   msg->get(atomSigPoolSize);
05826   if(atomSigPool) delete [] atomSigPool;
05827   atomSigPool = new AtomSignature[atomSigPoolSize];
05828   for(int i=0; i<atomSigPoolSize; i++)
05829       atomSigPool[i].unpack(msg);
05830 
05831   //get exclusions' signatures
05832   msg->get(exclSigPoolSize);
05833   if(exclSigPool) delete [] exclSigPool;
05834   exclSigPool = new ExclusionSignature[exclSigPoolSize];
05835   for(int i=0; i<exclSigPoolSize; i++)
05836       exclSigPool[i].unpack(msg);
05837  
05838   msg->get(numHydrogenGroups);      
05839   msg->get(maxHydrogenGroupSize);      
05840   msg->get(numMigrationGroups);      
05841   msg->get(maxMigrationGroupSize);      
05842   msg->get(isOccupancyValid);
05843   msg->get(isBFactorValid);
05844 
05845    //get names for atoms
05846   msg->get(atomNamePoolSize);
05847   atomNamePool = new char *[atomNamePoolSize];
05848   for(int i=0; i<atomNamePoolSize;i++) {
05849     int len;
05850     msg->get(len);
05851     atomNamePool[i] = nameArena->getNewArray(len+1);
05852     msg->get(len, atomNamePool[i]);
05853   }
05854   
05855   if(simParams->fixedAtomsOn){
05856     int numFixedAtomsSet;
05857     msg->get(numFixedAtoms);
05858     msg->get(numFixedAtomsSet);
05859     fixedAtomsSet = new AtomSetList(numFixedAtomsSet);
05860     msg->get(numFixedAtomsSet*sizeof(AtomSet), (char *)(fixedAtomsSet->begin()));
05861   } 
05862 
05863   if(simParams->constraintsOn){
05864     int numConstrainedAtomsSet;
05865     msg->get(numConstraints);
05866     msg->get(numConstrainedAtomsSet);
05867     constrainedAtomsSet = new AtomSetList(numConstrainedAtomsSet);
05868     msg->get(numConstrainedAtomsSet*sizeof(AtomSet), (char *)(constrainedAtomsSet->begin()));
05869   } 
05870 
05871 #else
05872   delete [] atoms;
05873   atoms= new Atom[numAtoms];  
05874   msg->get(numAtoms*sizeof(Atom), (char*)atoms);
05875 
05876   //  Get the bond information
05877   msg->get(numRealBonds);
05878   msg->get(numBonds);    
05879   if (numBonds)
05880   {
05881     delete [] bonds;
05882     bonds=new Bond[numBonds]; 
05883     msg->get(numBonds*sizeof(Bond), (char*)bonds);
05884   }  
05885 
05886   //  Get the angle information
05887   msg->get(numAngles);  
05888   if (numAngles)
05889   {
05890     delete [] angles;
05891     angles=new Angle[numAngles];  
05892     msg->get(numAngles*sizeof(Angle), (char*)angles);
05893   }  
05894 
05895   //  Get the dihedral information
05896   msg->get(numDihedrals);    
05897   if (numDihedrals)
05898   {
05899     delete [] dihedrals;
05900     dihedrals=new Dihedral[numDihedrals];  
05901     msg->get(numDihedrals*sizeof(Dihedral), (char*)dihedrals);
05902   }  
05903 
05904   if (simParams->sdScaling) {
05905     msg->get(num_alch_unpert_Bonds);
05906     alch_unpert_bonds=new Bond[num_alch_unpert_Bonds];
05907     msg->get(num_alch_unpert_Bonds*sizeof(Bond), (char*)alch_unpert_bonds);
05908 
05909     msg->get(num_alch_unpert_Angles);
05910     alch_unpert_angles=new Angle[num_alch_unpert_Angles];
05911     msg->get(num_alch_unpert_Angles*sizeof(Angle), (char*)alch_unpert_angles);
05912 
05913     msg->get(num_alch_unpert_Dihedrals);
05914     alch_unpert_dihedrals=new Dihedral[num_alch_unpert_Dihedrals];
05915     msg->get(num_alch_unpert_Dihedrals*sizeof(Dihedral), (char*)alch_unpert_dihedrals);
05916   }
05917   
05918   //  Get the improper information
05919   msg->get(numImpropers);
05920   if (numImpropers)
05921   {
05922     delete [] impropers;
05923     impropers=new Improper[numImpropers];  
05924     msg->get(numImpropers*sizeof(Improper), (char*)impropers);
05925   }
05926   
05927   //  Get the crossterm information
05928   msg->get(numCrossterms);
05929   if (numCrossterms)
05930   {
05931     delete [] crossterms;
05932     crossterms=new Crossterm[numCrossterms];  
05933     msg->get(numCrossterms*sizeof(Crossterm), (char*)crossterms);
05934   }
05935   
05936   //  Get the hydrogen bond donors
05937   msg->get(numDonors);  
05938   if (numDonors)
05939   {
05940     delete [] donors;
05941     donors=new Bond[numDonors];  
05942     msg->get(numDonors*sizeof(Bond), (char*)donors);
05943   }
05944   
05945   //  Get the hydrogen bond acceptors
05946   msg->get(numAcceptors);  
05947   if (numAcceptors)
05948   {
05949     delete [] acceptors;
05950     acceptors=new Bond[numAcceptors];  
05951     msg->get(numAcceptors*sizeof(Bond), (char*)acceptors);
05952   }
05953   
05954   //  Get the exclusion information 
05955   msg->get(numExclusions);  
05956   if (numExclusions)
05957   {
05958     delete [] exclusions;
05959     exclusions=new Exclusion[numExclusions];  
05960     msg->get(numExclusions*sizeof(Exclusion), (char*)exclusions);
05961   }
05962         
05963       //  Get the constraint information, if they are active
05964       if (simParams->constraintsOn)
05965       {
05966          msg->get(numConstraints);
05967 
05968          delete [] consIndexes;
05969          consIndexes = new int32[numAtoms];
05970          
05971          msg->get(numAtoms, consIndexes);
05972          
05973          if (numConstraints)
05974          {
05975            delete [] consParams;
05976            consParams = new ConstraintParams[numConstraints];
05977       
05978            msg->get(numConstraints*sizeof(ConstraintParams), (char*)consParams);
05979          }
05980       }
05981 #endif
05982 
05983       /* BEGIN gf */
05984       if (simParams->mgridforceOn)
05985       {
05986          DebugM(3, "Receiving gridforce info\n");
05987          
05988          msg->get(numGridforceGrids);
05989          
05990          DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n");
05991          
05992          delete [] numGridforces;
05993          numGridforces = new int[numGridforceGrids];
05994          
05995          delete [] gridfrcIndexes;      // Should I be deleting elements of these first?
05996          delete [] gridfrcParams;
05997          delete [] gridfrcGrid;
05998          gridfrcIndexes = new int32*[numGridforceGrids];
05999          gridfrcParams = new GridforceParams*[numGridforceGrids];
06000          gridfrcGrid = new GridforceGrid*[numGridforceGrids];
06001          
06002          int grandTotalGrids = 0;
06003          for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
06004              msg->get(numGridforces[gridnum]);
06005              
06006              gridfrcIndexes[gridnum] = new int32[numAtoms];
06007              msg->get(numAtoms, gridfrcIndexes[gridnum]);
06008          
06009              if (numGridforces[gridnum])
06010              {
06011                  gridfrcParams[gridnum] = new GridforceParams[numGridforces[gridnum]];
06012                  msg->get(numGridforces[gridnum]*sizeof(GridforceParams), (char*)gridfrcParams[gridnum]);
06013              }
06014              
06015              gridfrcGrid[gridnum] = GridforceGrid::unpack_grid(gridnum, msg);
06016              
06017              grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
06018          }
06019       }
06020       /* END gf */
06021       
06022       //  Get the stirring information, if stirring is  active
06023       if (simParams->stirOn)
06024       {
06025          msg->get(numStirredAtoms);
06026 
06027          delete [] stirIndexes;
06028          stirIndexes = new int32[numAtoms];
06029          
06030          msg->get(numAtoms, stirIndexes);
06031          
06032          if (numStirredAtoms)
06033          {
06034            delete [] stirParams;
06035            stirParams = new StirParams[numStirredAtoms];
06036       
06037            msg->get(numStirredAtoms*sizeof(StirParams), (char*)stirParams);
06038          }
06039       }
06040       
06041       //  Get the moving drag information, if it is active
06042       if (simParams->movDragOn) {
06043          msg->get(numMovDrag);
06044          delete [] movDragIndexes;
06045          movDragIndexes = new int32[numAtoms];
06046          msg->get(numAtoms, movDragIndexes);
06047          if (numMovDrag)
06048          {
06049            delete [] movDragParams;
06050            movDragParams = new MovDragParams[numMovDrag];
06051            msg->get(numMovDrag*sizeof(MovDragParams), (char*)movDragParams);
06052          }
06053       }
06054       
06055       //  Get the rotating drag information, if it is active
06056       if (simParams->rotDragOn) {
06057          msg->get(numRotDrag);
06058          delete [] rotDragIndexes;
06059          rotDragIndexes = new int32[numAtoms];
06060          msg->get(numAtoms, rotDragIndexes);
06061          if (numRotDrag)
06062          {
06063            delete [] rotDragParams;
06064            rotDragParams = new RotDragParams[numRotDrag];
06065            msg->get(numRotDrag*sizeof(RotDragParams), (char*)rotDragParams);
06066          }
06067       }
06068       
06069       //  Get the "constant" torque information, if it is active
06070       if (simParams->consTorqueOn) {
06071          msg->get(numConsTorque);
06072          delete [] consTorqueIndexes;
06073          consTorqueIndexes = new int32[numAtoms];
06074          msg->get(numAtoms, consTorqueIndexes);
06075          if (numConsTorque)
06076          {
06077            delete [] consTorqueParams;
06078            consTorqueParams = new ConsTorqueParams[numConsTorque];
06079            msg->get(numConsTorque*sizeof(ConsTorqueParams), (char*)consTorqueParams);
06080          }
06081       }
06082       
06083       // Get the constant force information, if it's active
06084       if (simParams->consForceOn)
06085       { msg->get(numConsForce);
06086         delete [] consForceIndexes;
06087         consForceIndexes = new int32[numAtoms];
06088         msg->get(numAtoms, consForceIndexes);
06089         if (numConsForce)
06090         { delete [] consForce;
06091           consForce = new Vector[numConsForce];
06092           msg->get(numConsForce*sizeof(Vector), (char*)consForce);
06093         }
06094       }
06095 
06096       if (simParams->excludeFromPressure) {
06097         exPressureAtomFlags = new int32[numAtoms];
06098         msg->get(numExPressureAtoms);
06099         msg->get(numAtoms, exPressureAtomFlags);
06100       }
06101 
06102 #ifndef MEM_OPT_VERSION
06103       //  Get the langevin parameters, if they are active
06104       if (simParams->langevinOn || simParams->tCoupleOn)
06105       {
06106         delete [] langevinParams;
06107         langevinParams = new Real[numAtoms];
06108 
06109         msg->get(numAtoms, langevinParams);
06110       }
06111 
06112       //  Get the fixed atoms, if they are active
06113       if (simParams->fixedAtomsOn)
06114       {
06115         delete [] fixedAtomFlags;
06116         fixedAtomFlags = new int32[numAtoms];
06117 
06118         msg->get(numFixedAtoms);
06119         msg->get(numAtoms, fixedAtomFlags);
06120         msg->get(numFixedRigidBonds);
06121       }
06122 
06123       if (simParams->qmForcesOn)
06124       {
06125         if( qmAtomGroup != 0)
06126             delete [] qmAtomGroup;
06127         qmAtomGroup = new Real[numAtoms];
06128         
06129         msg->get(numAtoms, qmAtomGroup);
06130         
06131         msg->get(qmNumQMAtoms);
06132         
06133         if( qmAtmChrg != 0)
06134             delete [] qmAtmChrg;
06135         qmAtmChrg = new Real[qmNumQMAtoms];
06136         
06137         msg->get(qmNumQMAtoms, qmAtmChrg);
06138         
06139         if( qmAtmIndx != 0)
06140             delete [] qmAtmIndx;
06141         qmAtmIndx = new int[qmNumQMAtoms];
06142         
06143         msg->get(qmNumQMAtoms, qmAtmIndx);
06144         
06145         msg->get(qmNoPC);
06146         
06147         msg->get(qmNumBonds);
06148         
06149         msg->get(qmMeNumBonds);
06150         
06151         if( qmMeMMindx != 0)
06152             delete [] qmMeMMindx;
06153         qmMeMMindx = new int[qmMeNumBonds];
06154         
06155         msg->get(qmMeNumBonds, qmMeMMindx);
06156         
06157         if( qmMeQMGrp != 0)
06158             delete [] qmMeQMGrp;
06159         qmMeQMGrp = new Real[qmMeNumBonds];
06160         
06161         msg->get(qmMeNumBonds, qmMeQMGrp);
06162         
06163         msg->get(qmPCFreq);
06164         
06165         msg->get(qmNumGrps);
06166         
06167         if( qmGrpID != 0)
06168             delete [] qmGrpID;
06169         qmGrpID = new Real[qmNumGrps];
06170         msg->get(qmNumGrps, qmGrpID);
06171         
06172         if( qmCustPCSizes != 0)
06173             delete [] qmCustPCSizes;
06174         qmCustPCSizes = new int[qmNumGrps];
06175         msg->get(qmNumGrps, qmCustPCSizes);
06176         
06177         msg->get(qmTotCustPCs);
06178         
06179         if( qmCustomPCIdxs != 0)
06180             delete [] qmCustomPCIdxs;
06181         qmCustomPCIdxs = new int[qmTotCustPCs];
06182         msg->get(qmTotCustPCs, qmCustomPCIdxs);
06183       }
06184     
06185 //fepb
06186       //receive fep atom info
06187       if (simParams->alchOn || simParams->lesOn || simParams->pairInteractionOn) {
06188         delete [] fepAtomFlags;
06189         fepAtomFlags = new unsigned char[numAtoms];
06190 
06191         msg->get(numFepInitial);
06192         msg->get(numFepFinal);
06193         msg->get(numAtoms*sizeof(unsigned char), (char*)fepAtomFlags);
06194       }
06195 //fepe
06196 
06197 //soluteScaling
06198       if (simParams->soluteScalingOn) {
06199         delete [] ssAtomFlags;
06200         delete [] ss_vdw_type;
06201         delete [] ss_index;
06202         ssAtomFlags = new unsigned char[numAtoms];
06203         ss_vdw_type = new int [params->get_num_vdw_params()];
06204         ss_index = new int [numAtoms];
06205         msg->get(numAtoms*sizeof(unsigned char), (char*)ssAtomFlags);
06206         msg->get(ss_num_vdw_params);
06207         msg->get(params->get_num_vdw_params()*sizeof(int), (char*)ss_vdw_type);
06208         msg->get(numAtoms*sizeof(int), (char*)ss_index);
06209       }
06210 //soluteScaling
06211 #ifdef OPENATOM_VERSION
06212       // This needs to be refactored into its own version
06213       if (simParams->openatomOn) {
06214         delete [] fepAtomFlags;
06215         fepAtomFlags = new unsigned char[numAtoms];
06216 
06217         msg->get(numFepInitial);
06218         msg->get(numAtoms*sizeof(unsigned char), (char*)fepAtomFlags);
06219 #endif //OPENATOM_VERSION
06220 
06221       // DRUDE: receive data read from PSF
06222       msg->get(is_lonepairs_psf);
06223       if (is_lonepairs_psf) {
06224         msg->get(numLphosts);
06225         delete[] lphosts;
06226         lphosts = new Lphost[numLphosts];
06227         msg->get(numLphosts*sizeof(Lphost), (char*)lphosts);
06228       }
06229       msg->get(is_drude_psf);
06230       if (is_drude_psf) {
06231         delete[] drudeConsts;
06232         drudeConsts = new DrudeConst[numAtoms];
06233         msg->get(numAtoms*sizeof(DrudeConst), (char*)drudeConsts);
06234         msg->get(numTholes);
06235         delete[] tholes;
06236         tholes = new Thole[numTholes];
06237         msg->get(numTholes*sizeof(Thole), (char*)tholes);
06238         msg->get(numAnisos);
06239         delete[] anisos;
06240         anisos = new Aniso[numAnisos];
06241         msg->get(numAnisos*sizeof(Aniso), (char*)anisos);
06242       }
06243       msg->get(numZeroMassAtoms);
06244       // DRUDE
06245 
06246   //LCPO
06247   if (simParams->LCPOOn) {
06248     delete [] lcpoParamType;
06249     lcpoParamType = new int[numAtoms];
06250     msg->get(numAtoms, (int*)lcpoParamType);
06251   }
06252 
06253   //Receive GromacsPairStuff -- JLai
06254 
06255   if (simParams->goGroPair) {
06256     msg->get(numLJPair);
06257     delete [] indxLJA;
06258     indxLJA = new int[numLJPair];
06259     msg->get(numLJPair,indxLJA);
06260     delete [] indxLJB;
06261     indxLJB = new int[numLJPair];
06262     msg->get(numLJPair,indxLJB);
06263     delete [] pairC6;
06264     pairC6 = new Real[numLJPair];    
06265     msg->get(numLJPair,pairC6);
06266     delete [] pairC12;
06267     pairC12 = new Real[numLJPair];
06268     msg->get(numLJPair,pairC12);
06269     delete [] gromacsPair_type;
06270     gromacsPair_type = new int[numLJPair];
06271     msg->get(numLJPair,gromacsPair_type);
06272     delete [] pointerToLJBeg;
06273     pointerToLJBeg = new int[numAtoms];
06274     msg->get((numAtoms),pointerToLJBeg);
06275     delete [] pointerToLJEnd;
06276     pointerToLJEnd = new int[numAtoms];
06277     msg->get((numAtoms),pointerToLJEnd);
06278     // JLai
06279     delete [] gromacsPair;
06280     gromacsPair = new GromacsPair[numLJPair];
06281     for(int i=0; i < numLJPair; i++) {
06282         gromacsPair[i].atom1 = indxLJA[i];
06283         gromacsPair[i].atom2 = indxLJB[i];
06284         gromacsPair[i].pairC6  = pairC6[i];
06285         gromacsPair[i].pairC12 = pairC12[i];
06286         gromacsPair[i].gromacsPair_type = gromacsPair_type[i];
06287     }
06288     //
06289     msg->get(numGaussPair);
06290     delete [] indxGaussA;
06291     indxGaussA = new int[numGaussPair];
06292     msg->get(numGaussPair,indxGaussA);
06293     delete [] indxGaussB;
06294     indxGaussB = new int[numGaussPair];
06295     msg->get(numGaussPair,indxGaussB);
06296     delete [] gA;
06297     gA = new Real[numGaussPair];
06298     msg->get(numGaussPair,gA);
06299     delete [] gMu1;
06300     gMu1 = new Real[numGaussPair];
06301     msg->get(numGaussPair,gMu1);
06302     delete [] giSigma1;
06303     giSigma1 = new Real[numGaussPair];
06304     msg->get(numGaussPair,giSigma1);
06305     delete [] gMu2;
06306     gMu2 = new Real[numGaussPair];
06307     msg->get(numGaussPair,gMu2);
06308     delete [] giSigma2;
06309     giSigma2 = new Real[numGaussPair];
06310     msg->get(numGaussPair,giSigma2);
06311     delete [] gRepulsive;
06312     gRepulsive = new Real[numGaussPair];
06313     msg->get(numGaussPair,gRepulsive);
06314     delete [] pointerToGaussBeg;
06315     pointerToGaussBeg = new int[numAtoms];
06316     msg->get((numAtoms),pointerToGaussBeg);
06317     delete [] pointerToGaussEnd;
06318     pointerToGaussEnd = new int[numAtoms];
06319     msg->get((numAtoms),pointerToGaussEnd);
06320     //
06321   }
06322 #endif
06323 
06324       //  Now free the message 
06325       delete msg;
06326 
06327 #ifdef MEM_OPT_VERSION
06328 
06329       build_excl_check_signatures();
06330 
06331     //set num{Calc}Tuples(Bonds,...,Impropers) to 0
06332     numBonds = numCalcBonds = 0;
06333     numAngles = numCalcAngles = 0;
06334     numDihedrals = numCalcDihedrals = 0;
06335     numImpropers = numCalcImpropers = 0;
06336     numCrossterms = numCalcCrossterms = 0;
06337     numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;  
06338     // JLai
06339     numLJPair = numCalcLJPair = 0;
06340     // End of JLai
06341 
06342 #else
06343 
06344       //  analyze the data and find the status of each atom
06345       build_atom_status();
06346       build_lists_by_atom();      
06347 
06348       
06349 #endif
06350 }

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

Referenced by Node::reloadCharges().

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

Definition at line 1453 of file Molecule.h.

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

01454   {
01455     return(rigidBondLengths[atomnum]);
01456   }

void Molecule::send_GoMolecule ( MOStream msg  ) 

Definition at line 1635 of file GoMolecule.C.

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

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

void Molecule::send_Molecule ( MOStream msg  ) 

Definition at line 5448 of file Molecule.C.

References alch_unpert_angles, alch_unpert_bonds, alch_unpert_dihedrals, SimParameters::alchOn, atomNamePool, atomSigPool, consForce, consForceIndexes, SimParameters::consForceOn, consTorqueIndexes, SimParameters::consTorqueOn, consTorqueParams, SimParameters::constraintsOn, DebugM, MOStream::end(), endi(), SimParameters::excludeFromPressure, SimParameters::fixedAtomsOn, gA, Parameters::get_num_vdw_params(), giSigma1, giSigma2, gMu1, gMu2, SimParameters::goGroPair, gRepulsive, gromacsPair_type, indxGaussA, indxGaussB, indxLJA, indxLJB, is_drude_psf, is_lonepairs_psf, isBFactorValid, isOccupancyValid, SimParameters::langevinOn, SimParameters::LCPOOn, SimParameters::lesOn, maxHydrogenGroupSize, maxMigrationGroupSize, SimParameters::mgridforceOn, SimParameters::movDragOn, num_alch_unpert_Angles, num_alch_unpert_Bonds, num_alch_unpert_Dihedrals, 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, numZeroMassAtoms, GridforceGrid::pack_grid(), pairC12, pairC6, SimParameters::pairInteractionOn, pointerToGaussBeg, pointerToGaussEnd, pointerToLJBeg, pointerToLJEnd, MOStream::put(), SimParameters::qmForcesOn, SimParameters::rotDragOn, SimParameters::sdScaling, SimParameters::soluteScalingOn, ss_index, ss_num_vdw_params, ss_vdw_type, SimParameters::stirOn, and SimParameters::tCoupleOn.

Referenced by Node::resendMolecule().

05448                                          {
05449 #ifdef MEM_OPT_VERSION
05450 //in the memory optimized version, only the atom signatures are broadcast
05451 //to other Nodes. --Chao Mei
05452 
05453   msg->put(numAtoms);
05454 
05455   msg->put(massPoolSize);
05456   msg->put(massPoolSize, atomMassPool);
05457 
05458   msg->put(chargePoolSize);
05459   msg->put(chargePoolSize, atomChargePool); 
05460 
05461   //put atoms' signatures
05462   msg->put(atomSigPoolSize);
05463   for(int i=0; i<atomSigPoolSize; i++)
05464       atomSigPool[i].pack(msg);
05465 
05466   //put atom's exclusion signatures
05467   msg->put(exclSigPoolSize);
05468   for(int i=0; i<exclSigPoolSize; i++)
05469       exclSigPool[i].pack(msg);
05470 
05471   msg->put(numHydrogenGroups);      
05472   msg->put(maxHydrogenGroupSize);      
05473   msg->put(numMigrationGroups);      
05474   msg->put(maxMigrationGroupSize);            
05475   msg->put(isOccupancyValid);
05476   msg->put(isBFactorValid);
05477   
05478   //put names for atoms
05479   msg->put(atomNamePoolSize);
05480   for(int i=0; i<atomNamePoolSize;i++) {
05481     int len = strlen(atomNamePool[i]);
05482     msg->put(len);
05483     msg->put(len*sizeof(char), atomNamePool[i]);
05484   } 
05485   
05486   if(simParams->fixedAtomsOn){
05487     int numFixedAtomsSet = fixedAtomsSet->size();
05488     msg->put(numFixedAtoms);
05489     msg->put(numFixedAtomsSet);
05490     msg->put(numFixedAtomsSet*sizeof(AtomSet), (char *)(fixedAtomsSet->begin()));
05491   }
05492 
05493   if (simParams->constraintsOn) {
05494     int numConstrainedAtomsSet = constrainedAtomsSet->size();
05495     msg->put(numConstraints);
05496     msg->put(numConstrainedAtomsSet);
05497     msg->put(numConstrainedAtomsSet*sizeof(AtomSet), (char *)(constrainedAtomsSet->begin()));
05498   }
05499     
05500 #else
05501   msg->put(numAtoms);
05502   msg->put(numAtoms*sizeof(Atom), (char*)atoms);
05503   
05504   //  Send the bond information
05505   msg->put(numRealBonds);
05506   msg->put(numBonds);
05507  
05508   if (numBonds)
05509   {
05510     msg->put(numBonds*sizeof(Bond), (char*)bonds);
05511   }
05512 
05513   //  Send the angle information
05514   msg->put(numAngles);  
05515   if (numAngles)
05516   {
05517     msg->put(numAngles*sizeof(Angle), (char*)angles);
05518   }  
05519 
05520   //  Send the dihedral information
05521   msg->put(numDihedrals);
05522   if (numDihedrals)
05523   {
05524     msg->put(numDihedrals*sizeof(Dihedral), (char*)dihedrals);
05525   }  
05526 
05527   if (simParams->sdScaling) {
05528     msg->put(num_alch_unpert_Bonds);
05529     msg->put(num_alch_unpert_Bonds*sizeof(Bond), (char*)alch_unpert_bonds);
05530 
05531     msg->put(num_alch_unpert_Angles);
05532     msg->put(num_alch_unpert_Angles*sizeof(Angle), (char*)alch_unpert_angles);
05533 
05534     msg->put(num_alch_unpert_Dihedrals);
05535     msg->put(num_alch_unpert_Dihedrals*sizeof(Dihedral), (char*)alch_unpert_dihedrals);
05536   }
05537 
05538   //  Send the improper information
05539   msg->put(numImpropers);  
05540   if (numImpropers)
05541   {
05542     msg->put(numImpropers*sizeof(Improper), (char*)impropers);
05543   }
05544 
05545   //  Send the crossterm information
05546   msg->put(numCrossterms);
05547   if (numCrossterms)
05548   {
05549     msg->put(numCrossterms*sizeof(Crossterm), (char*)crossterms);
05550   }
05551 
05552   // send the hydrogen bond donor information
05553   msg->put(numDonors);
05554   if(numDonors)
05555   {
05556     msg->put(numDonors*sizeof(Bond), (char*)donors);
05557   }
05558 
05559   // send the hydrogen bond acceptor information
05560   msg->put(numAcceptors);
05561   if(numAcceptors)
05562   {
05563     msg->put(numAcceptors*sizeof(Bond), (char*)acceptors);
05564   }
05565 
05566   //  Send the exclusion information  
05567   msg->put(numExclusions);
05568   if (numExclusions)
05569   {
05570     msg->put(numExclusions*sizeof(Exclusion), (char*)exclusions);
05571   }      
05572   //  Send the constraint information, if used
05573   if (simParams->constraintsOn)
05574   {
05575      msg->put(numConstraints);
05576      
05577      msg->put(numAtoms, consIndexes);
05578      
05579      if (numConstraints)
05580      {
05581        msg->put(numConstraints*sizeof(ConstraintParams), (char*)consParams);
05582      }
05583   }
05584 #endif
05585   
05586   /* BEGIN gf */
05587   // Send the gridforce information, if used
05588   if (simParams->mgridforceOn)
05589   {
05590     DebugM(3, "Sending gridforce info\n" << endi);
05591     msg->put(numGridforceGrids);
05592     
05593     for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
05594       msg->put(numGridforces[gridnum]);
05595       msg->put(numAtoms, gridfrcIndexes[gridnum]);
05596       if (numGridforces[gridnum])
05597       {
05598        msg->put(numGridforces[gridnum]*sizeof(GridforceParams), (char*)gridfrcParams[gridnum]);
05599       }
05600       GridforceGrid::pack_grid(gridfrcGrid[gridnum], msg);
05601     }
05602   }
05603   /* END gf */
05604   
05605   //  Send the stirring information, if used
05606   if (simParams->stirOn)
05607   {
05608      //CkPrintf ("DEBUG: putting numStirredAtoms..\n");
05609      msg->put(numStirredAtoms);
05610      //CkPrintf ("DEBUG: putting numAtoms,stirIndexes.. numAtoms=%d\n",numStirredAtoms);
05611      msg->put(numAtoms, stirIndexes);
05612      //CkPrintf ("DEBUG: if numStirredAtoms..\n");
05613      if (numStirredAtoms)
05614      {
05615        //CkPrintf ("DEBUG: big put, with (char*)stirParams\n");
05616        msg->put(numStirredAtoms*sizeof(StirParams), (char*)stirParams);
05617      }
05618   }
05619   
05620   
05621   //  Send the moving drag information, if used
05622   if (simParams->movDragOn) {
05623      msg->put(numMovDrag);
05624      msg->put(numAtoms, movDragIndexes);
05625      if (numMovDrag)
05626      {
05627        msg->put(numMovDrag*sizeof(MovDragParams), (char*)movDragParams);
05628      }
05629   }
05630   
05631   //  Send the rotating drag information, if used
05632   if (simParams->rotDragOn) {
05633      msg->put(numRotDrag);
05634      msg->put(numAtoms, rotDragIndexes);
05635      if (numRotDrag)
05636      {
05637        msg->put(numRotDrag*sizeof(RotDragParams), (char*)rotDragParams);
05638      }
05639   }
05640   
05641   //  Send the "constant" torque information, if used
05642   if (simParams->consTorqueOn) {
05643      msg->put(numConsTorque);
05644      msg->put(numAtoms, consTorqueIndexes);
05645      if (numConsTorque)
05646      {
05647        msg->put(numConsTorque*sizeof(ConsTorqueParams), (char*)consTorqueParams);
05648      }
05649   }
05650   
05651   // Send the constant force information, if used
05652   if (simParams->consForceOn)
05653   { msg->put(numConsForce);
05654     msg->put(numAtoms, consForceIndexes);
05655     if (numConsForce)
05656       msg->put(numConsForce*sizeof(Vector), (char*)consForce);
05657   }
05658   
05659   if (simParams->excludeFromPressure) {
05660     msg->put(numExPressureAtoms);
05661     msg->put(numAtoms, exPressureAtomFlags);
05662   }
05663   
05664 #ifndef MEM_OPT_VERSION
05665   //  Send the langevin parameters, if active
05666   if (simParams->langevinOn || simParams->tCoupleOn)
05667   {
05668     msg->put(numAtoms, langevinParams);
05669   }
05670   
05671   //  Send fixed atoms, if active
05672   if (simParams->fixedAtomsOn)
05673   {
05674     msg->put(numFixedAtoms);
05675     msg->put(numAtoms, fixedAtomFlags);
05676   msg->put(numFixedRigidBonds);
05677   }
05678   
05679   if (simParams->qmForcesOn)
05680   {
05681     msg->put(numAtoms, qmAtomGroup);
05682     msg->put(qmNumQMAtoms);
05683     msg->put(qmNumQMAtoms, qmAtmChrg);
05684     msg->put(qmNumQMAtoms, qmAtmIndx);
05685     msg->put(qmNoPC);
05686     msg->put(qmNumBonds);
05687     msg->put(qmMeNumBonds);
05688     msg->put(qmMeNumBonds, qmMeMMindx);
05689     msg->put(qmMeNumBonds, qmMeQMGrp);
05690     msg->put(qmPCFreq);
05691     msg->put(qmNumGrps);
05692     msg->put(qmNumGrps, qmGrpID);
05693     msg->put(qmNumGrps, qmCustPCSizes);
05694     msg->put(qmTotCustPCs);
05695     msg->put(qmTotCustPCs, qmCustomPCIdxs);
05696   }
05697   
05698   //fepb
05699   // send fep atom info
05700   if (simParams->alchOn || simParams->lesOn || simParams->pairInteractionOn) {
05701     msg->put(numFepInitial);
05702     msg->put(numFepFinal);
05703     msg->put(numAtoms*sizeof(char), (char*)fepAtomFlags);
05704   }
05705   //fepe
05706 
05707   if (simParams->soluteScalingOn) {
05708     msg->put(numAtoms*sizeof(char), (char*)ssAtomFlags);
05709     msg->put(ss_num_vdw_params);
05710     msg->put(params->get_num_vdw_params()*sizeof(int), (char*)ss_vdw_type);
05711     msg->put(numAtoms*sizeof(int), (char*)ss_index);
05712   }
05713 
05714   #ifdef OPENATOM_VERSION
05715   // needs to be refactored into its own openatom version
05716   if (simParams->openatomOn ) {
05717     msg->put(numFepInitial);
05718     msg->put(numAtoms*sizeof(char), (char*)fepAtomFlags);
05719   }
05720   #endif //OPENATOM_VERSION
05721   
05722   // DRUDE: send data read from PSF
05723   msg->put(is_lonepairs_psf);
05724   if (is_lonepairs_psf) {
05725     msg->put(numLphosts);
05726     msg->put(numLphosts*sizeof(Lphost), (char*)lphosts);
05727   }
05728   msg->put(is_drude_psf);
05729   if (is_drude_psf) {
05730     msg->put(numAtoms*sizeof(DrudeConst), (char*)drudeConsts);
05731     msg->put(numTholes);
05732     msg->put(numTholes*sizeof(Thole), (char*)tholes);
05733     msg->put(numAnisos);
05734     msg->put(numAnisos*sizeof(Aniso), (char*)anisos);
05735   }
05736   msg->put(numZeroMassAtoms);
05737   // DRUDE
05738 
05739   //LCPO
05740   if (simParams->LCPOOn) {
05741     msg->put(numAtoms, (int*)lcpoParamType);
05742   }
05743   
05744   //Send GromacsPairStuff -- JLai
05745   if (simParams->goGroPair) {
05746     msg->put(numLJPair);
05747     msg->put(numLJPair,indxLJA);
05748     msg->put(numLJPair,indxLJB);
05749     msg->put(numLJPair,pairC6);
05750     msg->put(numLJPair,pairC12);
05751     msg->put(numLJPair,gromacsPair_type);
05752     msg->put((numAtoms),pointerToLJBeg);
05753     msg->put((numAtoms),pointerToLJEnd);
05754     msg->put(numGaussPair);
05755     msg->put(numGaussPair,indxGaussA);
05756     msg->put(numGaussPair,indxGaussB);
05757     msg->put(numGaussPair,gA);
05758     msg->put(numGaussPair,gMu1);
05759     msg->put(numGaussPair,giSigma1);
05760     msg->put(numGaussPair,gMu2);
05761     msg->put(numGaussPair,giSigma2);
05762     msg->put(numGaussPair,gRepulsive);
05763     msg->put((numAtoms),pointerToGaussBeg);
05764     msg->put((numAtoms),pointerToGaussEnd);
05765   }
05766 #endif
05767 
05768   // Broadcast the message to the other nodes
05769   msg->end();
05770   delete msg;
05771 
05772 #ifdef MEM_OPT_VERSION
05773 
05774   build_excl_check_signatures();
05775 
05776   //set num{Calc}Tuples(Bonds,...,Impropers) to 0
05777   numBonds = numCalcBonds = 0;
05778   numAngles = numCalcAngles = 0;
05779   numDihedrals = numCalcDihedrals = 0;
05780   numImpropers = numCalcImpropers = 0;
05781   numCrossterms = numCalcCrossterms = 0;
05782   numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;  
05783   // JLai
05784   numLJPair = numCalcLJPair = 0;
05785   // End of JLai
05786 
05787 #else
05788 
05789   //  Now build arrays of indexes into these arrays by atom      
05790   build_lists_by_atom();
05791 
05792 #endif
05793 }

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

Definition at line 1285 of file Molecule.h.

References numGridforceGrids.

Referenced by Node::reloadGridforceGrid().

01286   {
01287       if (grid && gridnum >= 0 && gridnum < numGridforceGrids) {
01288           gridfrcGrid[gridnum] = grid;
01289           return 0;
01290       } else {
01291           return -1;
01292       }
01293   }

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

Definition at line 817 of file Molecule.h.

Referenced by NamdState::loadStructure().

00817 { qmReplaceAll = newReplaceAll; } ;

void Molecule::setBFactorData ( molfile_atom_t *  atomarray  ) 

Definition at line 3243 of file Molecule.C.

Referenced by Molecule().

03243                                                       {
03244     bfactor = new float[numAtoms];
03245     for(int i=0; i<numAtoms; i++) {
03246         bfactor[i] = atomarray[i].bfactor;
03247     }
03248 }

void Molecule::setOccupancyData ( molfile_atom_t *  atomarray  ) 

Definition at line 3236 of file Molecule.C.

Referenced by Molecule().

03236                                                         {
03237     occupancy = new float[numAtoms];
03238     for(int i=0; i<numAtoms; i++) {
03239         occupancy[i] = atomarray[i].occupancy;
03240     }
03241 }


Friends And Related Function Documentation

friend class AngleElem [friend]

Definition at line 215 of file Molecule.h.

friend class AnisoElem [friend]

Definition at line 219 of file Molecule.h.

friend class BondElem [friend]

Definition at line 214 of file Molecule.h.

friend class CrosstermElem [friend]

Definition at line 220 of file Molecule.h.

friend class DihedralElem [friend]

Definition at line 216 of file Molecule.h.

friend class ExclElem [friend]

Definition at line 213 of file Molecule.h.

friend class GromacsPairElem [friend]

Definition at line 222 of file Molecule.h.

friend class ImproperElem [friend]

Definition at line 217 of file Molecule.h.

friend class TholeElem [friend]

Definition at line 218 of file Molecule.h.

friend class WorkDistrib [friend]

Definition at line 224 of file Molecule.h.


Member Data Documentation

Definition at line 563 of file Molecule.h.

Referenced by NamdState::loadStructure().

Definition at line 564 of file Molecule.h.

Referenced by NamdState::loadStructure().

Definition at line 565 of file Molecule.h.

Referenced by NamdState::loadStructure().

ConsTorqueParams* Molecule::consTorqueParams

Definition at line 616 of file Molecule.h.

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

Definition at line 684 of file Molecule.h.

Referenced by build_go_arrays(), and goInit().

Definition at line 685 of file Molecule.h.

Referenced by build_go_arrays(), and goInit().

Definition at line 676 of file Molecule.h.

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

Definition at line 678 of file Molecule.h.

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

Definition at line 680 of file Molecule.h.

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

Definition at line 677 of file Molecule.h.

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

Definition at line 679 of file Molecule.h.

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

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

Definition at line 1585 of file Molecule.h.

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

Definition at line 653 of file Molecule.h.

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

Definition at line 654 of file Molecule.h.

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

Definition at line 652 of file Molecule.h.

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

Definition at line 650 of file Molecule.h.

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

Definition at line 645 of file Molecule.h.

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

Definition at line 655 of file Molecule.h.

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

Definition at line 656 of file Molecule.h.

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

Definition at line 681 of file Molecule.h.

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

Definition at line 669 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 639 of file Molecule.h.

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

Definition at line 674 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 675 of file Molecule.h.

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

Definition at line 665 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 666 of file Molecule.h.

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

Definition at line 461 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 462 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 1465 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 1465 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 569 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Number of anisotropic terms.

Definition at line 583 of file Molecule.h.

Referenced by AnisoElem::getMoleculePointers().

Definition at line 631 of file Molecule.h.

Referenced by Controller::compareChecksums().