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 1046 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().

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

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

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

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

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

References atoms.

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

01057   {      
01058       return(atoms[anum].vdw_type);
01059   }

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 6310 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().

06316 {
06317     PDB *kPDB;
06318     register int i;             //  Loop counters
06319     register int j;
06320     register int k;
06321 
06322     DebugM(3,  "Entered build_gridforce_params multi...\n");
06323 //     DebugM(3, "\tgridfrcfile = " << gridfrcfile->data << endi);
06324 //     DebugM(3, "\tgridfrccol = " << gridfrccol->data << endi);
06325     
06326     MGridforceParams* mgridParams = simParams->mgridforcelist.get_first();
06327     numGridforceGrids = 0;
06328     while (mgridParams != NULL) {
06329         numGridforceGrids++;
06330         mgridParams = mgridParams->next;
06331     }
06332     
06333     DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n");
06334     gridfrcIndexes = new int32*[numGridforceGrids];
06335     gridfrcParams = new GridforceParams*[numGridforceGrids];
06336     gridfrcGrid = new GridforceGrid*[numGridforceGrids];
06337     numGridforces = new int[numGridforceGrids];
06338     
06339     int grandTotalGrids = 0;    // including all subgrids
06340     
06341     mgridParams = simParams->mgridforcelist.get_first();
06342     for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
06343         int current_index=0;    //  Index into values used
06344         int kcol = 5;           //  Column to look for force constant in
06345         int qcol = 0;           //  Column for charge (default 0: use electric charge)
06346         Real kval = 0;          //  Force constant value retreived
06347         char filename[129];     //  PDB filename
06348         char potfilename[129];  //  Potential file name
06349         
06350         if (mgridParams == NULL) {
06351             NAMD_die("Problem with mgridParams!");
06352         }
06353         
06354         // Now load values from mgridforcelist object
06355         if (mgridParams->gridforceFile == NULL)
06356         {
06357             if ( ! initial_pdb ) NAMD_die("Initial PDB file unavailable, gridforceFile required.");
06358             kPDB = initial_pdb;
06359         }
06360         else
06361         {
06362             DebugM(4, "mgridParams->gridforceFile = " << mgridParams->gridforceFile << "\n" << endi);
06363             
06364             if ( (cwd == NULL) || (mgridParams->gridforceFile[0] == '/') )
06365             {
06366                 strcpy(filename, mgridParams->gridforceFile);
06367             }
06368             else
06369             {
06370                 strcpy(filename, cwd);
06371                 strcat(filename, mgridParams->gridforceFile);
06372             }
06373         
06374             kPDB = new PDB(filename);
06375             if ( kPDB == NULL )
06376             {
06377                 NAMD_die("Memory allocation failed in Molecule::build_gridforce_params");
06378             }
06379            
06380             if (kPDB->num_atoms() != numAtoms)
06381             {
06382                 NAMD_die("Number of atoms in grid force PDB doesn't match coordinate PDB");
06383             }
06384         }
06385 
06386         //  Get the column that the force constant is going to be in.  It
06387         //  can be in any of the 5 floating point fields in the PDB, according
06388         //  to what the user wants.  The allowable fields are X, Y, Z, O, or
06389         //  B which correspond to the 1st, 2nd, ... 5th floating point fields.
06390         //  The default is the 5th field, which is beta (temperature factor)
06391         if (mgridParams->gridforceCol == NULL)
06392         {
06393             kcol = 5;
06394         }
06395         else
06396         {
06397             if (strcasecmp(mgridParams->gridforceCol, "X") == 0)
06398             {
06399                 kcol=1;
06400             }
06401             else if (strcasecmp(mgridParams->gridforceCol, "Y") == 0)
06402             {
06403                 kcol=2;
06404             }
06405             else if (strcasecmp(mgridParams->gridforceCol, "Z") == 0)
06406             {
06407                 kcol=3;
06408             }
06409             else if (strcasecmp(mgridParams->gridforceCol, "O") == 0)
06410             {
06411                 kcol=4;
06412             }
06413             else if (strcasecmp(mgridParams->gridforceCol, "B") == 0)
06414             {
06415                 kcol=5;
06416             }
06417             else
06418             {
06419                 NAMD_die("gridforcecol must have value of X, Y, Z, O, or B");
06420             }
06421         }
06422     
06423         //  Get the column that the charge is going to be in.
06424         if (mgridParams->gridforceQcol == NULL)
06425         {
06426             qcol = 0;   // Default: don't read charge from file, use electric charge
06427         }
06428         else
06429         {
06430             if (strcasecmp(mgridParams->gridforceQcol, "X") == 0)
06431             {
06432                 qcol=1;
06433             }
06434             else if (strcasecmp(mgridParams->gridforceQcol, "Y") == 0)
06435             {
06436                 qcol=2;
06437             }
06438             else if (strcasecmp(mgridParams->gridforceQcol, "Z") == 0)
06439             {
06440                 qcol=3;
06441             }
06442             else if (strcasecmp(mgridParams->gridforceQcol, "O") == 0)
06443             {
06444                 qcol=4;
06445             }
06446             else if (strcasecmp(mgridParams->gridforceQcol, "B") == 0)
06447             {
06448                 qcol=5;
06449             }
06450             else
06451             {
06452                 NAMD_die("gridforcechargecol must have value of X, Y, Z, O, or B");
06453             }
06454         }
06455     
06456         if (kcol == qcol) {
06457             NAMD_die("gridforcecol and gridforcechargecol cannot have same value");
06458         }
06459 
06460     
06461         //  Allocate an array that will store an index into the constraint
06462         //  parameters for each atom.  If the atom is not constrained, its
06463         //  value will be set to -1 in this array.
06464         gridfrcIndexes[gridnum] = new int32[numAtoms];
06465        
06466         if (gridfrcIndexes[gridnum] == NULL)
06467         {
06468             NAMD_die("memory allocation failed in Molecule::build_gridforce_params()");
06469         }
06470         
06471         //  Loop through all the atoms and find out which ones are constrained
06472         for (i=0; i<numAtoms; i++)
06473         {
06474             //  Get the k value based on where we were told to find it
06475             switch (kcol)
06476             {
06477             case 1:
06478                 kval = (kPDB->atom(i))->xcoor();
06479                 break;
06480             case 2:
06481                 kval = (kPDB->atom(i))->ycoor();
06482                 break;
06483             case 3:
06484                 kval = (kPDB->atom(i))->zcoor();
06485                 break;
06486             case 4:
06487                 kval = (kPDB->atom(i))->occupancy();
06488                 break;
06489             case 5:
06490                 kval = (kPDB->atom(i))->temperaturefactor();
06491                 break;
06492             }
06493            
06494             if (kval > 0.0)
06495             {
06496                 //  This atom is constrained
06497                 gridfrcIndexes[gridnum][i] = current_index;
06498                 current_index++;
06499             }
06500             else
06501             {
06502                 //  This atom is not constrained
06503                 gridfrcIndexes[gridnum][i] = -1;
06504             }
06505         }
06506     
06507         if (current_index == 0)
06508         {
06509             //  Constraints were turned on, but there weren't really any constrained
06510             iout << iWARN << "NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" << endi;
06511         }
06512         else
06513         {
06514             //  Allocate an array to hold the constraint parameters
06515             gridfrcParams[gridnum] = new GridforceParams[current_index];
06516             if (gridfrcParams[gridnum] == NULL)
06517             {
06518                 NAMD_die("memory allocation failed in Molecule::build_gridforce_params");
06519             }
06520         }
06521     
06522         numGridforces[gridnum] = current_index;
06523 
06524         //  Loop through all the atoms and assign the parameters for those
06525         //  that are constrained
06526         for (i=0; i<numAtoms; i++)
06527         {
06528             if (gridfrcIndexes[gridnum][i] != -1)
06529             {
06530                 //  This atom has grid force, so get the k value again
06531                 switch (kcol)
06532                 {
06533                 case 1:
06534                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->xcoor();
06535                     break;
06536                 case 2:
06537                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->ycoor();
06538                     break;
06539                 case 3:
06540                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->zcoor();
06541                     break;
06542                 case 4:
06543                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->occupancy();
06544                     break;
06545                 case 5:
06546                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->temperaturefactor();
06547                     break;
06548                 }
06549             
06550                 //  Also get charge column
06551                 switch (qcol)
06552                 {
06553                 case 0:
06554 #ifdef MEM_OPT_VERSION
06555                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]];
06556 #else
06557                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge;
06558 #endif
06559                     break;
06560                 case 1:
06561                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->xcoor();
06562                     break;
06563                 case 2:
06564                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->ycoor();
06565                     break;
06566                 case 3:
06567                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->zcoor();
06568                     break;
06569                 case 4:
06570                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->occupancy();
06571                     break;
06572                 case 5:
06573                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->temperaturefactor();
06574                     break;
06575                 }
06576             }
06577         }
06578        
06579         //  If we had to create new PDB objects, delete them now
06580         if (mgridParams->gridforceFile != NULL)
06581         {
06582             delete kPDB;
06583         }
06584     
06585         //  Now we fill in our grid information
06586     
06587         // Open potential file
06588         if ( (cwd == NULL) || (mgridParams->gridforceVfile[0] == '/') )
06589         {
06590             strcpy(potfilename, mgridParams->gridforceVfile);
06591         }
06592         else
06593         {
06594             strcpy(potfilename, cwd);
06595             strcat(potfilename, mgridParams->gridforceVfile);
06596         }
06597     
06598 //        iout << iINFO << "Allocating grid " << gridnum
06599 //             << "\n" << endi;
06600         
06601         DebugM(3, "allocating GridforceGrid(" << gridnum << ")\n");
06602         gridfrcGrid[gridnum] = GridforceGrid::new_grid(gridnum, potfilename, simParams, mgridParams);
06603         
06604         grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
06605         DebugM(4, "grandTotalGrids = " << grandTotalGrids << "\n" << endi);
06606         
06607         // Finally, get next mgridParams pointer
06608         mgridParams = mgridParams->next;
06609     }
06610 }

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

Referenced by NamdState::loadStructure().

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

void Molecule::freeOccupancyData (  )  [inline]

Definition at line 1028 of file Molecule.h.

Referenced by NamdState::loadStructure().

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

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

Definition at line 1103 of file Molecule.h.

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

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

Definition at line 1066 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

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

Definition at line 1145 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01146       { 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 1114 of file Molecule.h.

References atomTypePool, and NAMD_die().

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

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

Definition at line 1063 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

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

Definition at line 1143 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01144       { return bondsByAtom[anum]; } 

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

Definition at line 1021 of file Molecule.h.

Referenced by wrap_coor_int().

01021 { return cluster[anum]; }

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

Definition at line 1022 of file Molecule.h.

Referenced by wrap_coor_int().

01022 { return clusterSize[anum]; }

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

Definition at line 1260 of file Molecule.h.

Referenced by ComputeRestraints::doForce().

01261   {
01262     k = consParams[consIndexes[atomnum]].k;
01263     refPos = consParams[consIndexes[atomnum]].refPos;
01264   }

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

Definition at line 1334 of file Molecule.h.

References consTorqueIndexes, and consTorqueParams.

Referenced by ComputeConsTorque::doForce().

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

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

Definition at line 1075 of file Molecule.h.

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

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

Definition at line 1151 of file Molecule.h.

01152       { return crosstermsByAtom[anum]; }  

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

Definition at line 872 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00872 {return cSMDcoffs;} ;

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

Definition at line 868 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00868 {return cSMDindex;} ;

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

Definition at line 867 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00867 { return cSMDindxLen;} ;

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

Definition at line 870 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00870 {return cSMDKs;} ;

int Molecule::get_cSMDnumInst (  )  [inline]

Definition at line 866 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00866 { return cSMDnumInst;} ;

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

Definition at line 869 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00869 {return cSMDpairs;} ;

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

Definition at line 871 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00871 {return cSMDVels;} ;

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

Definition at line 1072 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

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

Definition at line 1147 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

01148       { return dihedralsByAtom[anum]; }

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

Definition at line 1100 of file Molecule.h.

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

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

Definition at line 1171 of file Molecule.h.

Referenced by dumpbench().

01171                                                                {      
01172       return &all_exclusions[anum];             
01173   }

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

Definition at line 1110 of file Molecule.h.

References exclusions.

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

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

Definition at line 1153 of file Molecule.h.

References exclusionsByAtom.

01154       { return exclusionsByAtom[anum]; }

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

Definition at line 1383 of file Molecule.h.

References NAMD_die(), and simParams.

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

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

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

Referenced by ComputeNonbondedCUDA::build_exclusions().

01156       { return fullExclusionsByAtom[anum]; }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

01617 { 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 1605 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().

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

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

Definition at line 1274 of file Molecule.h.

References numGridforceGrids.

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

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

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

Definition at line 1268 of file Molecule.h.

Referenced by ComputeGridForce::do_calc().

01269   {
01270       k = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].k;
01271       q = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].q;
01272   }

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

Referenced by dumpbench().

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

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

Definition at line 1149 of file Molecule.h.

Referenced by dumpbench().

01150       { return impropersByAtom[anum]; }  

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

Definition at line 1079 of file Molecule.h.

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

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

Definition at line 1157 of file Molecule.h.

01158       { return modExclusionsByAtom[anum]; }

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

Definition at line 1319 of file Molecule.h.

Referenced by Sequencer::addMovDragToPosition().

01320   {
01321     v = movDragParams[movDragIndexes[atomnum]].v;
01322   }

Bool Molecule::get_noPC (  )  [inline]

Definition at line 854 of file Molecule.h.

Referenced by ComputeQM::initialize().

00854 { 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 818 of file Molecule.h.

00818 {return qmAtomGroup[indx]; } ;

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

Definition at line 865 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00865 { return qmcSMD;} ;

int* Molecule::get_qmCustomPCIdxs (  )  [inline]

Definition at line 863 of file Molecule.h.

Referenced by ComputeQM::initialize().

00863 { return qmCustomPCIdxs; } ;

int* Molecule::get_qmCustPCSizes (  )  [inline]

Definition at line 862 of file Molecule.h.

Referenced by ComputeQM::initialize().

00862 { return qmCustPCSizes; } ;

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

Definition at line 832 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00832 { return qmDummyBondVal; } ;

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

Definition at line 837 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00837 { return qmDummyElement; } ;

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

Definition at line 824 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00824 {return qmElementArray;} ;

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

Definition at line 835 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00835 { return qmGrpBonds; } ;

Real* Molecule::get_qmGrpChrg (  )  [inline]

Definition at line 828 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00828 {return qmGrpChrg; } ;

Real* Molecule::get_qmGrpID (  )  [inline]

Definition at line 827 of file Molecule.h.

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

00827 { return qmGrpID; } ;

Real* Molecule::get_qmGrpMult (  )  [inline]

Definition at line 829 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00829 {return qmGrpMult; } ;

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

Definition at line 833 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00833 { return qmGrpNumBonds; } ;

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

Definition at line 826 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00826 {return qmGrpSizes; } ;

int Molecule::get_qmLSSFreq (  )  [inline]

Definition at line 845 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00845 { return qmLSSFreq; } ;

int* Molecule::get_qmLSSIdxs (  )  [inline]

Definition at line 843 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00843 { return qmLSSIdxs;} ;

Mass* Molecule::get_qmLSSMass (  )  [inline]

Definition at line 844 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00844 { return qmLSSMass; } ;

int* Molecule::get_qmLSSRefIDs (  )  [inline]

Definition at line 847 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00847 { return qmLSSRefIDs ; } ;

Mass* Molecule::get_qmLSSRefMass (  )  [inline]

Definition at line 849 of file Molecule.h.

00849 {return qmLSSRefMass; } ;

int* Molecule::get_qmLSSRefSize (  )  [inline]

Definition at line 848 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00848 { return qmLSSRefSize ; } ;

int Molecule::get_qmLSSResSize (  )  [inline]

Definition at line 846 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00846 { return qmLSSResidueSize; } ;

int* Molecule::get_qmLSSSize (  )  [inline]

Definition at line 842 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00842 { return qmLSSSize;} ;

int* Molecule::get_qmMeMMindx (  )  [inline]

Definition at line 856 of file Molecule.h.

Referenced by ComputeQM::initialize().

00856 { return qmMeMMindx; } ;

int Molecule::get_qmMeNumBonds (  )  [inline]

Definition at line 855 of file Molecule.h.

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

00855 { return qmMeNumBonds; } ;

Real* Molecule::get_qmMeQMGrp (  )  [inline]

Definition at line 857 of file Molecule.h.

Referenced by ComputeQM::initialize().

00857 { return qmMeQMGrp; } ;

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

Definition at line 834 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00834 { return qmMMBond; } ;

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

Definition at line 836 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00836 { return qmMMBondedIndx; } ;

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

Definition at line 839 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00839 { return qmMMChargeTarget; } ;

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

Definition at line 840 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

00840 { return qmMMNumTargs; } ;

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

Definition at line 850 of file Molecule.h.

00850 { return qmClassicSolv;} ;

int Molecule::get_qmNumBonds (  )  [inline]

Definition at line 831 of file Molecule.h.

00831 { return qmNumBonds; } ;

int Molecule::get_qmNumGrps (  )  const [inline]

Definition at line 825 of file Molecule.h.

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

00825 {return qmNumGrps; } ;

int Molecule::get_qmPCFreq (  )  [inline]

Definition at line 859 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00859 { return qmPCFreq; } ;

Bool Molecule::get_qmReplaceAll (  )  [inline]

Definition at line 852 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

00852 {return qmReplaceAll; } ;

int Molecule::get_qmTotCustPCs (  )  [inline]

Definition at line 861 of file Molecule.h.

Referenced by ComputeQM::initialize().

00861 { 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 1325 of file Molecule.h.

Referenced by Sequencer::addRotDragToPosition().

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

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

Definition at line 1349 of file Molecule.h.

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

01350         {
01351                 return(ssAtomFlags[anum]);
01352         }

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

Definition at line 1300 of file Molecule.h.

01301   {
01302     refPos = stirParams[stirIndexes[atomnum]].refPos;
01303   }

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

Definition at line 1312 of file Molecule.h.

Referenced by ComputeStir::doForce().

01313   {
01314     return stirParams[stirIndexes[atomnum]].startTheta;
01315   }

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

Definition at line 1106 of file Molecule.h.

01106 {return acceptors;}

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

Definition at line 1089 of file Molecule.h.

Referenced by buildAngleData().

01089 {return angles;}

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

Definition at line 1088 of file Molecule.h.

Referenced by buildBondData().

01088 {return bonds;}

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

Definition at line 1092 of file Molecule.h.

Referenced by buildCrosstermData().

01092 {return crossterms;}

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

Definition at line 1091 of file Molecule.h.

Referenced by buildDihedralData().

01091 {return dihedrals;}

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

Definition at line 1105 of file Molecule.h.

01105 {return donors;}

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

Definition at line 1090 of file Molecule.h.

Referenced by buildImproperData().

01090 {return impropers;}

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

Definition at line 1096 of file Molecule.h.

01096 { return lphosts; }

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

Definition at line 480 of file Molecule.h.

References drude_constants::alpha.

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

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

Definition at line 489 of file Molecule.h.

Referenced by buildAtomData().

00489 { return atomNames; }

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

Definition at line 488 of file Molecule.h.

References atoms.

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

00488 { return atoms; }

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

Definition at line 492 of file Molecule.h.

Referenced by buildAtomData().

00492 { return atomSegResids; }

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

Definition at line 1030 of file Molecule.h.

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

01030 { return (const float *)bfactor; }

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

Definition at line 476 of file Molecule.h.

Referenced by HomePatch::setLcpoType().

00476                                  {
00477     return lcpoParamType;
00478   }

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

Definition at line 1026 of file Molecule.h.

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

01026 { 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 1244 of file Molecule.h.

References consTorqueIndexes, FALSE, and numConsTorque.

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

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

Definition at line 1195 of file Molecule.h.

References FALSE, and numConstraints.

Referenced by ComputeRestraints::doForce().

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

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

Definition at line 1445 of file Molecule.h.

References numExPressureAtoms.

Referenced by Sequencer::langevinPiston().

01446   {
01447     return (numExPressureAtoms && exPressureAtomFlags[atomnum]);
01448   }

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

Definition at line 1405 of file Molecule.h.

References numFixedAtoms.

Referenced by WorkDistrib::createAtomLists().

01406   {
01407     return (numFixedAtoms && fixedAtomFlags[atomnum]);
01408   }

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

Definition at line 1179 of file Molecule.h.

References FALSE, and numGridforceGrids.

Referenced by ComputeGridForce::do_calc().

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

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

Definition at line 1212 of file Molecule.h.

References FALSE, and numMovDrag.

Referenced by Sequencer::addMovDragToPosition().

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

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

Definition at line 1228 of file Molecule.h.

References FALSE, and numRotDrag.

Referenced by Sequencer::addRotDragToPosition().

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

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

Definition at line 1426 of file Molecule.h.

References FALSE, and numStirredAtoms.

Referenced by ComputeStir::doForce().

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

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

Definition at line 1441 of file Molecule.h.

References numFixedAtoms.

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

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

Referenced by WorkDistrib::createAtomLists().

01295   {
01296     return(langevinParams ? langevinParams[atomnum] : 0.);
01297   }

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

Definition at line 521 of file Molecule.h.

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

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

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

int Molecule::num_fixed_atoms (  )  const [inline]

Definition at line 495 of file Molecule.h.

References numFixedAtoms, and simParams.

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

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

int Molecule::num_fixed_groups (  )  const [inline]

Definition at line 501 of file Molecule.h.

References num_fixed_atoms(), and numFixedGroups.

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

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

int64_t Molecule::num_group_deg_freedom (  )  const [inline]

Definition at line 508 of file Molecule.h.

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

Referenced by Controller::receivePressure().

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

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

Prepares Live Solvent Selection

Data gathering from PDB to find QM atom and bond info

Multiplicity of each QM group

Charge of each QM group

Populate arrays that are used throughout the the calculations.

Overides Link Atom element with user selection.

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

Live Solvent Selection

Custom Point Charge selection

Topology preparation

Definition at line 109 of file MoleculeQM.C.

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

Referenced by NamdState::loadStructure().

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

void Molecule::print_atoms ( Parameters params  ) 

Definition at line 5282 of file Molecule.C.

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

Referenced by NamdState::loadStructure().

05283 {
05284 #ifdef MEM_OPT_VERSION
05285     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
05286 #else   
05287   register int i;
05288   Real sigma;
05289   Real epsilon;
05290   Real sigma14;
05291   Real epsilon14;
05292 
05293   DebugM(2,"ATOM LIST\n" \
05294       << "******************************************\n" \
05295                   << "NUM  NAME TYPE RES  MASS    CHARGE CHARGE   FEP-CHARGE"  \
05296       << "SIGMA   EPSILON SIGMA14 EPSILON14\n" \
05297         << endi);
05298 
05299   for (i=0; i<numAtoms; i++)
05300   {
05301     params->get_vdw_params(&sigma, &epsilon, &sigma14, &epsilon14, 
05302         atoms[i].vdw_type);
05303 
05304     DebugM(2,i+1 << " " << atomNames[i].atomname  \
05305               << " " << atomNames[i].atomtype << " " \
05306               << atomNames[i].resname  << " " << atoms[i].mass  \
05307         << " " << atoms[i].charge << " " << sigma \
05308         << " " << epsilon << " " << sigma14 \
05309         << " " << epsilon14 << "\n" \
05310         << endi);
05311   }
05312 #endif  
05313 }

void Molecule::print_bonds ( Parameters params  ) 

Definition at line 5325 of file Molecule.C.

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

Referenced by NamdState::loadStructure().

05326 {
05327 #ifdef MEM_OPT_VERSION
05328     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
05329 #else   
05330   register int i;
05331   Real k;
05332   Real x0;
05333 
05334   DebugM(2,"BOND LIST\n" << "********************************\n" \
05335       << "ATOM1 ATOM2 TYPE1 TYPE2      k        x0" \
05336       << endi);
05337 
05338   for (i=0; i<numBonds; i++)
05339   {
05340     params->get_bond_params(&k, &x0, bonds[i].bond_type);
05341 
05342     DebugM(2,bonds[i].atom1+1 << " " \
05343        << bonds[i].atom2+1 << " "   \
05344        << atomNames[bonds[i].atom1].atomtype << " "  \
05345        << atomNames[bonds[i].atom2].atomtype << " " << k \
05346        << " " << x0 << endi);
05347   }
05348   
05349 #endif  
05350 }

void Molecule::print_exclusions (  ) 

Definition at line 5362 of file Molecule.C.

References DebugM, and endi().

Referenced by NamdState::loadStructure().

05363 {
05364 #ifdef MEM_OPT_VERSION
05365     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
05366 #else
05367   register int i;
05368 
05369   DebugM(2,"EXPLICIT EXCLUSION LIST\n" \
05370       << "********************************\n" \
05371             << "ATOM1 ATOM2 " \
05372       << endi);
05373 
05374   for (i=0; i<numExclusions; i++)
05375   {
05376     DebugM(2,exclusions[i].atom1+1 << "  " \
05377        << exclusions[i].atom2+1 << endi);
05378   }
05379 #endif
05380 }

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

01307   {
01308     stirParams[stirIndexes[atomnum]].startTheta = theta;
01309   }

void Molecule::read_alch_unpert_angles ( FILE *  fd  ) 

Definition at line 1728 of file Molecule.C.

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

01728                                                {
01729   int atom_nums[3];  //  Atom numbers for the three atoms
01730   register int j;    //  Loop counter
01731   int num_read=0;    //  Number of angles read so far
01732 
01733   alch_unpert_angles=new Angle[num_alch_unpert_Angles];
01734 
01735   if (alch_unpert_angles == NULL) {
01736     NAMD_die("memory allocation failed in Molecule::read_alch_unpert_angles");
01737   }
01738 
01739   while (num_read < num_alch_unpert_Angles) {
01740     for (j=0; j<3; j++) {
01741       atom_nums[j]=NAMD_read_int(fd, "ANGLES")-1;
01742 
01743       if (atom_nums[j] >= numAtoms) {
01744         char err_msg[128];
01745         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);
01746         NAMD_die(err_msg);
01747       }
01748     }
01749 
01750     alch_unpert_angles[num_read].atom1=atom_nums[0];
01751     alch_unpert_angles[num_read].atom2=atom_nums[1];
01752     alch_unpert_angles[num_read].atom3=atom_nums[2];
01753 
01754     ++num_read;
01755   }
01756   return;
01757 }

void Molecule::read_alch_unpert_bonds ( FILE *  fd  ) 

Definition at line 1605 of file Molecule.C.

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

01605                                               {
01606   int atom_nums[2];  // Atom indexes for the bonded atoms
01607   register int j;      // Loop counter
01608   int num_read=0;    // Number of bonds read so far
01609 
01610   alch_unpert_bonds=new Bond[num_alch_unpert_Bonds];
01611 
01612   if (alch_unpert_bonds == NULL) {
01613     NAMD_die("memory allocations failed in Molecule::read_alch_unpert_bonds");
01614   }
01615 
01616   while (num_read < num_alch_unpert_Bonds) {
01617     for (j=0; j<2; j++) {
01618       atom_nums[j]=NAMD_read_int(fd, "BONDS")-1;
01619 
01620       if (atom_nums[j] >= numAtoms) {
01621         char err_msg[128];
01622 
01623         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);
01624         NAMD_die(err_msg);
01625       }
01626     }
01627 
01628     Bond *b = &(alch_unpert_bonds[num_read]);
01629     b->atom1=atom_nums[0];
01630     b->atom2=atom_nums[1];
01631 
01632     ++num_read;
01633   }
01634   return;
01635 }

void Molecule::read_alch_unpert_dihedrals ( FILE *  fd  ) 

Definition at line 1876 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.

01876                                                   {
01877   int atom_nums[4];  // The 4 atom indexes
01878   int num_read=0;    // number of dihedrals read so far
01879   register int j;    // loop counter
01880 
01881   alch_unpert_dihedrals = new Dihedral[num_alch_unpert_Dihedrals];
01882 
01883   if (alch_unpert_dihedrals == NULL) {
01884     NAMD_die("memory allocation failed in Molecule::read_alch_unpert_dihedrals");
01885   }
01886 
01887   while (num_read < num_alch_unpert_Dihedrals) {
01888     for (j=0; j<4; j++) {
01889       atom_nums[j]=NAMD_read_int(fd, "DIHEDRALS")-1;
01890 
01891       if (atom_nums[j] >= numAtoms) {
01892         char err_msg[128];
01893 
01894         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);
01895         NAMD_die(err_msg);
01896       }
01897     }
01898 
01899     alch_unpert_dihedrals[num_read].atom1=atom_nums[0];
01900     alch_unpert_dihedrals[num_read].atom2=atom_nums[1];
01901     alch_unpert_dihedrals[num_read].atom3=atom_nums[2];
01902     alch_unpert_dihedrals[num_read].atom4=atom_nums[3];
01903 
01904     num_read++;
01905   }
01906   return;
01907 }

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 5748 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, 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().

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

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

Referenced by Node::reloadCharges().

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

Definition at line 1451 of file Molecule.h.

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

01452   {
01453     return(rigidBondLengths[atomnum]);
01454   }

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 5393 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, 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().

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

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

Definition at line 1283 of file Molecule.h.

References numGridforceGrids.

Referenced by Node::reloadGridforceGrid().

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

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

Definition at line 815 of file Molecule.h.

Referenced by NamdState::loadStructure().

00815 { qmReplaceAll = newReplaceAll; } ;

void Molecule::setBFactorData ( molfile_atom_t *  atomarray  ) 

Definition at line 3188 of file Molecule.C.

Referenced by Molecule().

03188                                                       {
03189     bfactor = new float[numAtoms];
03190     for(int i=0; i<numAtoms; i++) {
03191         bfactor[i] = atomarray[i].bfactor;
03192     }
03193 }

void Molecule::setOccupancyData ( molfile_atom_t *  atomarray  ) 

Definition at line 3181 of file Molecule.C.

Referenced by Molecule().

03181                                                         {
03182     occupancy = new float[numAtoms];
03183     for(int i=0; i<numAtoms; i++) {
03184         occupancy[i] = atomarray[i].occupancy;
03185     }
03186 }


Friends And Related Function Documentation

friend class AngleElem [friend]

Definition at line 215 of file Molecule.h.

friend class AnisoElem [friend]

Definition at line 219 of file Molecule.h.

friend class BondElem [friend]

Definition at line 214 of file Molecule.h.

friend class CrosstermElem [friend]

Definition at line 220 of file Molecule.h.

friend class DihedralElem [friend]

Definition at line 216 of file Molecule.h.

friend class ExclElem [friend]

Definition at line 213 of file Molecule.h.

friend class GromacsPairElem [friend]

Definition at line 222 of file Molecule.h.

friend class ImproperElem [friend]

Definition at line 217 of file Molecule.h.

friend class TholeElem [friend]

Definition at line 218 of file Molecule.h.

friend class WorkDistrib [friend]

Definition at line 224 of file Molecule.h.


Member Data Documentation

Definition at line 561 of file Molecule.h.

Referenced by NamdState::loadStructure().

Definition at line 562 of file Molecule.h.

Referenced by NamdState::loadStructure().

Definition at line 563 of file Molecule.h.

Referenced by NamdState::loadStructure().

ConsTorqueParams* Molecule::consTorqueParams

Definition at line 614 of file Molecule.h.

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

Definition at line 682 of file Molecule.h.

Referenced by build_go_arrays(), and goInit().

Definition at line 683 of file Molecule.h.

Referenced by build_go_arrays(), and goInit().

Definition at line 674 of file Molecule.h.

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

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 675 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().

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

Definition at line 1583 of file Molecule.h.

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

Definition at line 651 of file Molecule.h.

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

Definition at line 652 of file Molecule.h.

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

Definition at line 650 of file Molecule.h.

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

Definition at line 648 of file Molecule.h.

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

Definition at line 643 of file Molecule.h.

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

Definition at line 653 of file Molecule.h.

Referenced by build_go_sigmas2(), get_go_force2(), 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 679 of file Molecule.h.

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

Definition at line 667 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 637 of file Molecule.h.

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

Definition at line 672 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 673 of file Molecule.h.

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

Definition at line 663 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 664 of file Molecule.h.

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

Definition at line 459 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 460 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 1463 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 1463 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Definition at line 567 of file Molecule.h.

Referenced by receive_Molecule(), and send_Molecule().

Number of anisotropic terms.

Definition at line 581 of file Molecule.h.

Referenced by AnisoElem::getMoleculePointers().

Definition at line 629 of file Molecule.h.

Referenced by Controller::compareChecksums().