Molecule Class Reference

#include <Molecule.h>

Public Member Functions
void	compute_LJcorrection ()
BigReal	getEnergyTailCorr (const BigReal, const int)
BigReal	getVirialTailCorr (const BigReal)
int const *	getLcpoParamType ()
BigReal	GetAtomAlpha (int i) const
Atom *	getAtoms () const
AtomNameInfo *	getAtomNames () const
AtomSegResInfo *	getAtomSegResInfo () 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 Real *	get_qmAtomGroup () const
Real	get_qmAtomGroup (int indx) const
Real *	get_qmAtmChrg ()
const int *	get_qmAtmIndx ()
int	get_numQMAtoms ()
const char *const *	get_qmElements ()
int	get_qmNumGrps () const
const int *	get_qmGrpSizes ()
Real *	get_qmGrpID ()
Real *	get_qmGrpChrg ()
Real *	get_qmGrpMult ()
int	get_qmNumBonds ()
const BigReal *	get_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 ()
Mass *	get_qmLSSMass ()
int	get_qmLSSFreq ()
int	get_qmLSSResSize ()
int *	get_qmLSSRefIDs ()
int *	get_qmLSSRefSize ()
Mass *	get_qmLSSRefMass ()
std::map< int, int > &	get_qmMMSolv ()
Bool	get_qmReplaceAll ()
Bool	get_noPC ()
int	get_qmMeNumBonds ()
int *	get_qmMeMMindx ()
Real *	get_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 Real *	get_cSMDKs ()
const Real *	get_cSMDVels ()
const Real *	get_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 *, AmberParm7Reader::Ambertoppar *)
void	read_parm (AmberParm7Reader::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
Bond *	get_bond (int bnum) const
Angle *	get_angle (int anum) const
Improper *	get_improper (int inum) const
Dihedral *	get_dihedral (int dnum) const
Crossterm *	get_crossterm (int inum) const
Lphost *	get_lphost (int atomid) const
Bond *	getAllBonds () const
Angle *	getAllAngles () const
Improper *	getAllImpropers () const
Dihedral *	getAllDihedrals () const
Crossterm *	getAllCrossterms () const
Lphost *	getAllLphosts () const
Bond *	get_donor (int dnum) const
Bond *	get_acceptor (int dnum) const
Bond *	getAllDonors () const
Bond *	getAllAcceptors () const
Exclusion *	get_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
int32 *	get_bonds_for_atom (int anum)
int32 *	get_angles_for_atom (int anum)
int32 *	get_dihedrals_for_atom (int anum)
int32 *	get_impropers_for_atom (int anum)
int32 *	get_crossterms_for_atom (int anum)
int32 *	get_exclusions_for_atom (int anum)
const int32 *	get_full_exclusions_for_atom (int anum) const
const int32 *	get_mod_exclusions_for_atom (int anum) const
int	checkexcl (int atom1, int atom2) const
const ExclusionCheck *	get_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
GridforceGrid *	get_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
Bond *	alch_unpert_bonds
Angle *	alch_unpert_angles
Dihedral *	alch_unpert_dihedrals
int	numLonepairs
	Number of lone pairs. More...
int	numDrudeAtoms
	Number of Drude particles. More...
int	numTholes
	Number of Thole terms. More...
int	numAnisos
	Number of anisotropic terms. More...
int	numLphosts
	Number of lone pair host records in PSF. More...
int	numZeroMassAtoms
	Number of atoms with zero mass. More...
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
int32 *	consForceIndexes
Vector *	consForce
int32 *	consTorqueIndexes
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
int32 *	atomChainTypes
int32 *	goSigmaIndices
int32 *	goResidIndices
Real *	goSigmas
bool *	goWithinCutoff
Real *	goCoordinates
int *	goResids
PDB *	goPDB
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
Real *	pairC6
Real *	pairC12
int *	gromacsPair_type
int *	pointerToGaussBeg
int *	pointerToGaussEnd
int	numGaussPair
int *	indxGaussA
int *	indxGaussB
Real *	gA
Real *	gMu1
Real *	giSigma1
Real *	gMu2
Real *	giSigma2
Real *	gRepulsive
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 151 of file Molecule.h.

Constructor & Destructor Documentation

Molecule::Molecule	(	SimParameters *	simParams,
		Parameters *	param
	)

Definition at line 429 of file Molecule.C.

{
  initialize(simParams,param);
}

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

Definition at line 442 of file Molecule.C.

References SimParameters::LCPOOn, and SimParameters::useCompressedPsf.

{
  initialize(simParams,param);

#ifdef MEM_OPT_VERSION
  if(simParams->useCompressedPsf)
      read_mol_signatures(filename, param, cfgList);
#else
        read_psf_file(filename, param); 
 //LCPO
  if (simParams->LCPOOn)
    assignLCPOTypes( 0 );
#endif      
 }

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

Definition at line 464 of file Molecule.C.

References SimParameters::LCPOOn, and NAMD_die().

{
#ifdef MEM_OPT_VERSION
  NAMD_die("Sorry, plugin IO is not supported in the memory optimized version.");
#else
    initialize(simParams, param);
    numAtoms = natoms;
    int optflags = MOLFILE_BADOPTIONS;
    molfile_atom_t *atomarray = (molfile_atom_t *) malloc(natoms*sizeof(molfile_atom_t));
    memset(atomarray, 0, natoms*sizeof(molfile_atom_t));

    //1a. read basic atoms information
    int rc = pIOHdl->read_structure(pIOFileHdl, &optflags, atomarray);
    if (rc != MOLFILE_SUCCESS && rc != MOLFILE_NOSTRUCTUREDATA) {
        free(atomarray);
        NAMD_die("ERROR: plugin failed reading structure data");
    }
    if(optflags == MOLFILE_BADOPTIONS) {
        free(atomarray);
        NAMD_die("ERROR: plugin didn't initialize optional data flags");
    }
    if(optflags & MOLFILE_OCCUPANCY) {
        setOccupancyData(atomarray);
    }
    if(optflags & MOLFILE_BFACTOR) {
        setBFactorData(atomarray);
    }
    //1b. load basic atoms information to the molecule object
    plgLoadAtomBasics(atomarray);    
    free(atomarray);

    //2a. read bonds
    //indices are one-based in read_bonds
    int *from, *to;
    float *bondorder;
    int *bondtype, nbondtypes;
    char **bondtypename;
    if(pIOHdl->read_bonds!=NULL) {
        if(pIOHdl->read_bonds(pIOFileHdl, &numBonds, &from, &to, &bondorder,
                                 &bondtype, &nbondtypes, &bondtypename)){
            NAMD_die("ERROR: failed reading bond information.");
        }
    }    
    //2b. load bonds information to the molecule object
    if(numBonds!=0) {
        plgLoadBonds(from,to);
    }

    //3a. read other bonded structures
    int *plgAngles, *plgDihedrals, *plgImpropers, *plgCterms;
    int ctermcols, ctermrows;
    int *angletypes, numangletypes, *dihedraltypes, numdihedraltypes;
    int *impropertypes, numimpropertypes; 
    char **angletypenames, **dihedraltypenames, **impropertypenames;

    plgAngles=plgDihedrals=plgImpropers=plgCterms=NULL;
    if(pIOHdl->read_angles!=NULL) {
        if(pIOHdl->read_angles(pIOFileHdl,
                  &numAngles, &plgAngles,
                  &angletypes, &numangletypes, &angletypenames,
                  &numDihedrals, &plgDihedrals,
                  &dihedraltypes, &numdihedraltypes, &dihedraltypenames,
                  &numImpropers, &plgImpropers,
                  &impropertypes, &numimpropertypes, &impropertypenames,
                  &numCrossterms, &plgCterms, &ctermcols, &ctermrows)) {
            NAMD_die("ERROR: failed reading angle information.");
        }
    }
    //3b. load other bonded structures to the molecule object
    if(numAngles!=0) plgLoadAngles(plgAngles);
    if(numDihedrals!=0) plgLoadDihedrals(plgDihedrals);
    if(numImpropers!=0) plgLoadImpropers(plgImpropers);
    if(numCrossterms!=0) plgLoadCrossterms(plgCterms);

  numRealBonds = numBonds;
  build_atom_status();
  //LCPO
  if (simParams->LCPOOn)
    assignLCPOTypes( 2 );
#endif
}

Molecule::Molecule	(	SimParameters *	,
		Parameters *	,
		Ambertoppar *
	)

Molecule::Molecule	(	SimParameters *	,
		Parameters *	,
		AmberParm7Reader::Ambertoppar *
	)

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

Molecule::~Molecule

(

)

Definition at line 558 of file Molecule.C.

References atoms, atomSigPool, exclusions, and exclusionsByAtom.

{
  /*  Check to see if each array was ever allocated.  If it was   */
  /*  then free it            */
  if (atoms != NULL)
    delete [] atoms;

  if (atomNames != NULL)
  {
    // subarrarys allocated from arena - automatically deleted
    delete [] atomNames;
  }
  delete nameArena;

  if (resLookup != NULL)
    delete resLookup;

  // DRUDE: free arrays read from PSF
  if (drudeConsts != NULL) delete [] drudeConsts;
  if (lphosts != NULL) delete [] lphosts;
  if (anisos != NULL) delete [] anisos;
  if (tholes != NULL) delete [] tholes;
  if (lphostIndexes != NULL) delete [] lphostIndexes;
  // DRUDE

  //LCPO
  if (lcpoParamType != NULL) delete [] lcpoParamType;

  #ifdef MEM_OPT_VERSION
  if(eachAtomSig) delete [] eachAtomSig;
  if(atomSigPool) delete [] atomSigPool;
  #else
  if (bonds != NULL)
    delete [] bonds;

  if (angles != NULL)
    delete [] angles;

  if (dihedrals != NULL)
    delete [] dihedrals;

  if (impropers != NULL)
    delete [] impropers;

  if (crossterms != NULL)
    delete [] crossterms;

  if (exclusions != NULL)
    delete [] exclusions;
  #endif

  if (donors != NULL)
    delete [] donors;

  if (acceptors != NULL)
    delete [] acceptors;  

  #ifdef MEM_OPT_VERSION
  if(exclSigPool) delete [] exclSigPool;
  if(exclChkSigPool) delete [] exclChkSigPool;
  if(eachAtomExclSig) delete [] eachAtomExclSig;
  if(fixedAtomsSet) delete fixedAtomsSet;
  if(constrainedAtomsSet) delete constrainedAtomsSet;
  #else
  if (bondsByAtom != NULL)
       delete [] bondsByAtom;
  
  if (anglesByAtom != NULL)
       delete [] anglesByAtom;
  
  if (dihedralsByAtom != NULL)
       delete [] dihedralsByAtom;
  
  if (impropersByAtom != NULL)
       delete [] impropersByAtom;
  
  if (crosstermsByAtom != NULL)
       delete [] crosstermsByAtom;  

  if (exclusionsByAtom != NULL)
       delete [] exclusionsByAtom;
  
  if (fullExclusionsByAtom != NULL)
       delete [] fullExclusionsByAtom;
  
  if (modExclusionsByAtom != NULL)
       delete [] modExclusionsByAtom;
  
  if (all_exclusions != NULL)
       delete [] all_exclusions;

  // JLai
  if (gromacsPairByAtom != NULL)
      delete [] gromacsPairByAtom;
  // End of JLai

  // DRUDE
  if (tholesByAtom != NULL)
       delete [] tholesByAtom;
  if (anisosByAtom != NULL)
       delete [] anisosByAtom;
  // DRUDE
  #endif

  //LCPO
  if (lcpoParamType != NULL)
    delete [] lcpoParamType;

  if (fixedAtomFlags != NULL)
       delete [] fixedAtomFlags;

  if (stirIndexes != NULL)
    delete [] stirIndexes;


  #ifdef MEM_OPT_VERSION
  if(clusterSigs != NULL){      
      delete [] clusterSigs;
  }  
  #else
  if (cluster != NULL)
       delete [] cluster;  
  #endif
  if (clusterSize != NULL)
       delete [] clusterSize;

  if (exPressureAtomFlags != NULL)
       delete [] exPressureAtomFlags;

  if (rigidBondLengths != NULL)
       delete [] rigidBondLengths;

//fepb
  if (fepAtomFlags != NULL)
       delete [] fepAtomFlags;
  if (alch_unpert_bonds != NULL)
       delete [] alch_unpert_bonds;
  if (alch_unpert_angles != NULL)
       delete [] alch_unpert_angles;
  if (alch_unpert_dihedrals != NULL)
       delete [] alch_unpert_dihedrals;
//fepe

//soluteScaling
  if (ssAtomFlags != NULL)
       delete [] ssAtomFlags;
  if (ss_vdw_type != NULL)
       delete [] ss_vdw_type;
  if (ss_index != NULL)
       delete [] ss_index;
//soluteScaling

  if (qmAtomGroup != NULL)
       delete [] qmAtomGroup;
  
  if (qmAtmIndx != NULL)
       delete [] qmAtmIndx;
  
  if (qmAtmChrg != NULL)
       delete [] qmAtmChrg;
  
  
  if (qmGrpNumBonds != NULL)
       delete [] qmGrpNumBonds;
  
  if (qmGrpSizes != NULL)
       delete [] qmGrpSizes;
  
  if (qmDummyBondVal != NULL)
       delete [] qmDummyBondVal;
  
  if (qmMMNumTargs != NULL)
       delete [] qmMMNumTargs;
  
  if (qmGrpID != NULL)
       delete [] qmGrpID;
  
  if (qmGrpChrg != NULL)
       delete [] qmGrpChrg;
  
  if (qmGrpMult != NULL)
       delete [] qmGrpMult;
  
  if (qmMeMMindx != NULL)
       delete [] qmMeMMindx;
  
  if (qmMeQMGrp != NULL)
       delete [] qmMeQMGrp;
  
  if (qmLSSSize != NULL)
       delete [] qmLSSSize;
  
  if (qmLSSIdxs != NULL)
       delete [] qmLSSIdxs;
  
  if (qmLSSMass != NULL)
       delete [] qmLSSMass;
  
  if (qmLSSRefSize != NULL)
       delete [] qmLSSRefSize;
  
  if (qmLSSRefIDs != NULL)
       delete [] qmLSSRefIDs;
  
  if (qmLSSRefMass != NULL)
       delete [] qmLSSRefMass;
  
  if (qmMMBond != NULL) {
      for(int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
          if (qmMMBond[grpIndx] != NULL)
              delete [] qmMMBond[grpIndx];
      }
      delete [] qmMMBond;
  }
  
  if (qmGrpBonds != NULL) {
      for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
          if (qmGrpBonds[grpIndx] != NULL)
              delete [] qmGrpBonds[grpIndx];
      }
      delete [] qmGrpBonds;
  }
  
  if (qmMMBondedIndx != NULL) {
      for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
          if (qmMMBondedIndx[grpIndx] != NULL)
              delete [] qmMMBondedIndx[grpIndx];
      }
      delete [] qmMMBondedIndx;
  }
  
  if (qmMMChargeTarget != NULL) {
      for(int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
          if (qmMMChargeTarget[grpIndx] != NULL)
              delete [] qmMMChargeTarget[grpIndx];
      }
      delete [] qmMMChargeTarget;
  }
  
    if (qmElementArray != NULL){
        for(int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
            if (qmElementArray[atmInd] != NULL)
                delete [] qmElementArray[atmInd];
        }
        delete [] qmElementArray;
    }
    
    if (qmDummyElement != NULL){
        for(int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
            if (qmDummyElement[atmInd] != NULL)
                delete [] qmDummyElement[atmInd];
        }
        delete [] qmDummyElement;
    }

    if (qmCustPCSizes != NULL){
        delete [] qmCustPCSizes;
    }
    
    if (qmCustomPCIdxs != NULL){
        delete [] qmCustomPCIdxs;
    }
    
    if (cSMDindex != NULL) {
      for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
          if (cSMDindex[grpIndx] != NULL)
              delete [] cSMDindex[grpIndx];
      }
      delete [] cSMDindex;
    }
    
    if (cSMDpairs != NULL) {
      for(int instIndx = 0 ; instIndx < cSMDnumInst; instIndx++) {
          if (cSMDpairs[instIndx] != NULL)
              delete [] cSMDpairs[instIndx];
      }
      delete [] cSMDpairs;
    }
    
    if (cSMDindxLen != NULL)
        delete [] cSMDindxLen;
    if (cSMDKs != NULL)
        delete [] cSMDKs;
    if (cSMDVels != NULL)
        delete [] cSMDVels;
    if (cSMDcoffs != NULL)
        delete [] cSMDcoffs;
    
  #ifndef MEM_OPT_VERSION
  delete arena;
  delete exclArena;
  #endif
}

Member Function Documentation

Real Molecule::atomcharge ( int  anum ) const

inline

Definition at line 1052 of file Molecule.h.

References atoms, and charge.

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

  {
    #ifdef MEM_OPT_VERSION
    return atomChargePool[eachAtomCharge[anum]];
    #else
    return(atoms[anum].charge);
    #endif
  }

Real Molecule::atommass ( int  anum ) const

inline

Definition at line 1042 of file Molecule.h.

References atoms.

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

  {
    #ifdef MEM_OPT_VERSION
    return atomMassPool[eachAtomMass[anum]];
    #else
    return(atoms[anum].mass);
    #endif
  }

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

Definition at line 1537 of file GoMolecule.C.

References bond::atom1, angle::atom1, dihedral::atom1, bond::atom2, angle::atom2, dihedral::atom2, angle::atom3, dihedral::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().

{
  int bondNum;   //  Bonds in bonded list
  int angleNum;  //  Angles in angle list
  int dihedralNum;   //  Dihedrals in dihedral list
  int *bonds;
  int *angles;
  int *dihedrals;
  Bond *bond;     //  Temporary bond structure
  Angle *angle;   //  Temporary angle structure
  Dihedral *dihedral; //  Temporary dihedral structure

  DebugM(2,"atoms_1to4(" << atom1 << "," << atom2 << ")" << std::endl);

  bonds = get_bonds_for_atom(atom1);
  bondNum = *bonds;
  while(bondNum != -1) {
    bond = get_bond(bondNum);
    DebugM(2,"bond  atom1:" << bond->atom1 << ", atom2:" << bond->atom2 << std::endl);
    if (atom2 == bond->atom1 || atom2 == bond->atom2) {
      return TRUE;
    }
    bondNum = *(++bonds);
  }

  bonds = get_bonds_for_atom(atom2);
  bondNum = *bonds;
  while(bondNum != -1) {
    bond = get_bond(bondNum);
    DebugM(2,"bond  atom1:" << bond->atom1 << ", atom2:" << bond->atom2 << std::endl);
    if (atom1 == bond->atom1 || atom1 == bond->atom2) {
      return TRUE;
    }
    bondNum = *(++bonds);
  }

  angles = get_angles_for_atom(atom1);
  angleNum = *angles;
  while(angleNum != -1) {
    angle = get_angle(angleNum);
    DebugM(2,"angle  atom1:" << angle->atom1 << ", atom2:" << angle->atom2 << ", atom3:" << angle->atom3 << std::endl);
    if (atom2 == angle->atom1 || atom2 == angle->atom2 || atom2 == angle->atom3) {
      return TRUE;
    }
    angleNum = *(++angles);
  }

  angles = get_angles_for_atom(atom2);
  angleNum = *angles;
  while(angleNum != -1) {
    angle = get_angle(angleNum);
    DebugM(2,"angle  atom1:" << angle->atom1 << ", atom2:" << angle->atom2 << ", atom3:" << angle->atom3 << std::endl);
    if (atom1 == angle->atom1 || atom1 == angle->atom2 || atom1 == angle->atom3) {
      return TRUE;
    }
    angleNum = *(++angles);
  }

  dihedrals = get_dihedrals_for_atom(atom1);
  dihedralNum = *dihedrals;
  while(dihedralNum != -1) {
    dihedral = get_dihedral(dihedralNum);
    DebugM(2,"dihedral  atom1:" << dihedral->atom1 << ", atom2:" << dihedral->atom2 << ", atom3:" << dihedral->atom3 << ", atom4:" << dihedral->atom4 << std::endl);
    if (atom2 == dihedral->atom1 || atom2 == dihedral->atom2 \
        || atom2 == dihedral->atom3 || atom2 == dihedral->atom4) {
      return TRUE;
    }
    dihedralNum = *(++dihedrals);
  }

  dihedrals = get_dihedrals_for_atom(atom2);
  dihedralNum = *dihedrals;
  while(dihedralNum != -1) {
    dihedral = get_dihedral(dihedralNum);
    DebugM(2,"dihedral  atom1:" << dihedral->atom1 << ", atom2:" << dihedral->atom2 << ", atom3:" << dihedral->atom3 << ", atom4:" << dihedral->atom4 << std::endl);
    if (atom1 == dihedral->atom1 || atom1 == dihedral->atom2 \
        || atom1 == dihedral->atom3 || atom1 == dihedral->atom4) {
      return TRUE;
    }
    dihedralNum = *(++dihedrals);
  }
  
  return FALSE;
}

Index Molecule::atomvdwtype ( int  anum ) const

inline

Definition at line 1062 of file Molecule.h.

References atoms.

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

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

void Molecule::build_alch_unpert_bond_lists

(

char *

)

Referenced by NamdState::loadStructure().

void Molecule::build_constant_forces

(

char *

)

Referenced by NamdState::loadStructure().

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

Referenced by NamdState::loadStructure().

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

Referenced by NamdState::loadStructure().

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

Referenced by NamdState::loadStructure().

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

Referenced by NamdState::loadStructure().

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

Referenced by NamdState::loadStructure().

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

Referenced by NamdState::loadStructure().

void Molecule::build_go_arrays	(	StringList *	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, NAMD_die(), NAMD_FILENAME_BUFFER_SIZE, StringList::next, PDB::num_atoms(), numAtoms, numGoAtoms, PDBAtom::residueseq(), PDBAtom::xcoor(), PDBAtom::ycoor(), and PDBAtom::zcoor().

Referenced by NamdState::loadStructure().

{
  DebugM(3,"->build_go_arrays" << std::endl);
  //PDB *goPDB;      //  Pointer to PDB object to use
  int bcol = 4;      //  Column that data is in
  int32 chainType = 0;      //  b value from PDB file
  int i;         //  Loop counter
  int j;         //  Loop counter
  BigReal atomAtomDist;     //  Distance between two atoms -- JLai put back 
  int resid1;    //  Residue ID for first atom
  int resid2;    //  Residue ID for second atom
  int residDiff;     //  Difference between resid1 and resid2
  int goIndex;       //  Index into the goCoordinates array
  int goIndx;        //  Index into the goCoordinates array
  char filename[NAMD_FILENAME_BUFFER_SIZE];    //  Filename
  
  //JLai
  BigReal nativeEnergy, *native;
  BigReal nonnativeEnergy, *nonnative;
  nativeEnergy = 0;
  nonnativeEnergy = 0;
  native = &nativeEnergy;
  nonnative = &nonnativeEnergy;

  long nativeContacts = 0;
  long nonnativeContacts = 0;

  //  Get the PDB object that contains the Go coordinates.  If
  //  the user gave another file name, use it.  Otherwise, just use
  //  the PDB file that has the initial coordinates.  
  if (goCoordFile == NULL)
    {
      //goPDB = initial_pdb;
      NAMD_die("Error: goCoordFile is NULL - build_go_arrays");
    }
  else
  {
    if (goCoordFile->next != NULL)
      {
        NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
      }
    
    if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
      {
        strcpy(filename, goCoordFile->data);
      }
    else
      {
        strcpy(filename, cwd);
        strcat(filename, goCoordFile->data);
      }
    
    goPDB = new PDB(filename);
    if ( goPDB == NULL )
      {
        NAMD_die("goPDB memory allocation failed in Molecule::build_go_arrays");
      }
    
    if (goPDB->num_atoms() != numAtoms)
      {
        NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
      }
  }
  
  //  Allocate the array to hold Go atom indices into the sigma array
  goSigmaIndices = new int32[numAtoms];
  
  if (goSigmaIndices == NULL) {
    NAMD_die("goSigmaIndices memory allocation failed in Molecule::build_go_arrays");
  }
  
  numGoAtoms = 0;
  
  //  Loop through all the atoms and get the Go chain types
  for (i=0; i<numAtoms; i++) {
    chainType = (int32)(goPDB->atom(i))->occupancy();
    if ( chainType != 0 ) {
      DebugM(3,"build_go_arrays - atom:" << i << std::endl);
      goSigmaIndices[i] = numGoAtoms;
      numGoAtoms++;
    }
    else {
      goSigmaIndices[i] = -1;
    }
  }

  // Allocate the array to hold the sigma values for each Go atom pair
  //  Allocate the array to hold the chain types
  atomChainTypes = new int32[numGoAtoms];

  if (atomChainTypes == NULL) {
    NAMD_die("atomChainTypes memory allocation failed in Molecule::build_go_arrays");
  }

  // Allocate the array to hold (x,y,z) coordinates for all of the Go atoms
  goCoordinates = new Real[numGoAtoms*3];

  if (goCoordinates == NULL) {
    NAMD_die("goCoordinates memory allocation failed in Molecule::build_go_arrays");
  }

  goResids = new int[numGoAtoms];

  // Allocate the array to hold PDB residu IDs for all of the Go atoms
  if (goResids == NULL) {
    NAMD_die("goResids memory allocation failed in Molecule::build_go_arrays");
  }
  
  for (i=0; i<numAtoms; i++) {
    goIndex = goSigmaIndices[i];
    if (goIndex != -1) {
      //  Assign the chainType value!
      //  Get the chainType from the occupancy field
      atomChainTypes[goIndex] = (int32)(goPDB->atom(i))->occupancy();
      goCoordinates[goIndex*3] = goPDB->atom(i)->xcoor();
      goCoordinates[goIndex*3 + 1] = goPDB->atom(i)->ycoor();
      goCoordinates[goIndex*3 + 2] = goPDB->atom(i)->zcoor();
      goResids[goIndex] = goPDB->atom(i)->residueseq();
    }
  }
      // JLai
  energyNative = 0;
  energyNonnative = 0;
  //printf("INIT ENERGY: (N) %f (NN) %f\n", energyNative, energyNonnative);
  for (i=0; i<numAtoms-1; i++) {
    goIndex = goSigmaIndices[i];
    if (goIndex != -1) {
      for (j=i+1; j<numAtoms;j++) {
        goIndx = goSigmaIndices[j];
        if (goIndx != -1) {
          resid1 = (goPDB->atom(i))->residueseq();
          resid2 = (goPDB->atom(j))->residueseq();
          residDiff = resid2 - resid1;
          if (residDiff < 0) residDiff = -residDiff;
          if (atomChainTypes[goIndex] && atomChainTypes[goIndx] &&
              !(this->go_restricted(atomChainTypes[goIndex],atomChainTypes[goIndx],residDiff)) &&
              !atoms_1to4(i,j)) {
            atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
                                pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
                                pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
            if (atomAtomDist <= this->get_go_cutoff(atomChainTypes[goIndex],atomChainTypes[goIndx]) ) {
              nativeContacts++;
            } else {
              nonnativeContacts++;
            }    
          }
        }
      }
    }
  }
  iout << iINFO << "Number of UNIQUE    native contacts: " << nativeContacts     << "\n" << endi;
  iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts  << "\n" << endi;
  
  //  If we had to create a PDB object, delete it now
  if (goCoordFile != NULL) {
    delete goPDB;
  }
  
  return;
}

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

                                            {
  iout << iINFO << "Building Go Parameters" << "\n" << endi;
#ifdef MEM_OPT_VERSION
  NAMD_die("Go forces are not supported in memory-optimized builds.");
#else
  build_lists_by_atom();
#endif
  int iterator = 0;
    do
    {
      iout << iINFO << "Reading Go File: " << iterator << "\n" << endi;
      read_go_file(g->data);
      g = g->next;
      iterator++;
    } while ( g != NULL && iterator < 100);    
}

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, NAMD_die(), NAMD_FILENAME_BUFFER_SIZE, StringList::next, PDB::num_atoms(), numAtoms, and numGoAtoms.

Referenced by NamdState::loadStructure().

{
  DebugM(3,"->build_go_sigmas" << std::endl);
  PDB *goPDB;      //  Pointer to PDB object to use
  int bcol = 4;      //  Column that data is in
  int32 chainType = 0;      //  b value from PDB file
  int i;         //  Loop counter
  int j;         //  Loop counter
  int resid1;    //  Residue ID for first atom
  int resid2;    //  Residue ID for second atom
  int residDiff;     //  Difference between resid1 and resid2
  Real sigma;    //  Sigma calculated for a Go atom pair
  Real atomAtomDist;     //  Distance between two atoms
  Real exp_a;            //  First exponent in L-J like formula
  Real exp_b;            //  Second exponent in L-J like formula
  char filename[NAMD_FILENAME_BUFFER_SIZE];    //  Filename
  
  //JLai
  BigReal nativeEnergy, *native;
  BigReal nonnativeEnergy, *nonnative;
  nativeEnergy = 0;
  nonnativeEnergy = 0;
  native = &nativeEnergy;
  nonnative = &nonnativeEnergy;

  long nativeContacts = 0;
  long nonnativeContacts = 0;

  //  Get the PDB object that contains the Go coordinates.  If
  //  the user gave another file name, use it.  Otherwise, just use
  //  the PDB file that has the initial coordinates.  
  if (goCoordFile == NULL)
    {
      //goPDB = initial_pdb;
      NAMD_die("Error: goCoordFile is NULL - build_go_sigmas");
    }
  else
  {
    if (goCoordFile->next != NULL)
      {
        NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
      }
    
    if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
      {
        strcpy(filename, goCoordFile->data);
      }
    else
      {
        strcpy(filename, cwd);
        strcat(filename, goCoordFile->data);
      }
    
    goPDB = new PDB(filename);
    if ( goPDB == NULL )
      {
        NAMD_die("Memory allocation failed in Molecule::build_go_sigmas");
      }
    
    if (goPDB->num_atoms() != numAtoms)
      {
        NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
      }
  }
  //  Allocate the array to hold the chain types
  atomChainTypes = new int32[numAtoms];
  //  Allocate the array to hold Go atom indices into the sigma array
  goSigmaIndices = new int32[numAtoms];
  
  if (atomChainTypes == NULL) {
    NAMD_die("memory allocation failed in Molecule::build_go_sigmas");
  }
  
  numGoAtoms = 0;
  
  //  Loop through all the atoms and get the Go chain types
  for (i=0; i<numAtoms; i++) {
    //  Get the chainType from the occupancy field
    chainType = (int32)(goPDB->atom(i))->occupancy();
    //  Assign the chainType value
    if ( chainType != 0 ) {
      //DebugM(3,"build_go_sigmas - atom:" << i << ", chainType:" << chainType << std::endl);
      atomChainTypes[i] = chainType;
      goSigmaIndices[i] = numGoAtoms;
      numGoAtoms++;
    }
    else {
      atomChainTypes[i] = 0;
      goSigmaIndices[i] = -1;
    }
    //printf("CT: %d %d %d %d\n",i,numGoAtoms,atomChainTypes[i],goSigmaIndices[i]);
  }

  // Allocate the array to hold the sigma values for each Go atom pair
  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] = -1.0;
      goWithinCutoff[i*numGoAtoms + j] = false;
    }
  }
  //  Loop through all atom pairs and calculate sigmas
  DebugM(3,"    numAtoms=" << numAtoms << std::endl);
  for (i=0; i<numAtoms; i++) {
    //DebugM(3,"    i=" << i << std::endl);
    resid1 = (goPDB->atom(i))->residueseq();
    //DebugM(3,"    resid1=" << resid1 << std::endl);
    //if ( goSigmaIndices[i] != -1) {
    //  goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[i]] = 0.0;
    //}
     for (j=i+1; j<numAtoms; j++) {
      //DebugM(3,"    j=" << j << std::endl);
      resid2 = (goPDB->atom(j))->residueseq();
      //printf("GSIi %d %d %d\n",i,numAtoms,goSigmaIndices[i]);
      //printf("SAN CHECK: %d\n",goSigmaIndices[37]);
      //printf("GSIj %d %d %d\n",j,numAtoms,goSigmaIndices[j]);
      //printf("ATOMS_1to4 %d\n",!atoms_1to4(i,j));
      //DebugM(3,"    resid2=" << resid2 << std::endl);
      //  if goSigmaIndices aren't defined, don't set anything in goSigmas
      if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 && !atoms_1to4(i,j) ) {
        //printf("TAKING DIFFERENCE\n");
        residDiff = resid2 - resid1;
        //printf("RESIDDIFF %d\n",residDiff);
        if (residDiff < 0) residDiff = -residDiff;
        //printf("RESIDDIFF2 %d\n",residDiff);
        //  if this is a Go atom pair that is not restricted
        //    calculate sigma
        //  sigmas are initially set to -1.0 if the atom pair fails go_restricted
        //printf("CHECKING LOOPING\n");
        if ( atomChainTypes[i] && atomChainTypes[j] &&
             !(this->go_restricted(atomChainTypes[i],atomChainTypes[j],residDiff)) ) {
          atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
                              pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
                              pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
          if ( atomAtomDist <= this->get_go_cutoff(atomChainTypes[i],atomChainTypes[j]) ) {
            exp_a = this->get_go_exp_a(atomChainTypes[i],atomChainTypes[j]);
            exp_b = this->get_go_exp_b(atomChainTypes[i],atomChainTypes[j]);
            sigma = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
            goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = sigma;
            goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = sigma;
            goWithinCutoff[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = true;
            goWithinCutoff[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = true;
            nativeContacts++;
          } else {
            goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = 0.0;
            goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = 0.0;
            nonnativeContacts++;
          }
        } else {
          goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = -1.0;
          goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = -1.0;
        }
      } 
    }
  }

  iout << iINFO << "Number of UNIQUE    native contacts: " << nativeContacts << "\n" << endi;
  iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts << "\n" << endi;
  
  //  If we had to create a PDB object, delete it now
  if (goCoordFile != NULL) {
    delete goPDB;
  }
  
  return;
}

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

Definition at line 747 of file GoMolecule.C.

References go_pair::A, ResizeArray< T >::add(), PDB::atom(), atomChainTypes, atoms_1to4(), go_pair::B, ResizeArray< T >::begin(), StringList::data, DebugM, ResizeArray< T >::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, NAMD_die(), NAMD_FILENAME_BUFFER_SIZE, StringList::next, PDB::num_atoms(), numAtoms, numGoAtoms, numLJPair, pointerToGoBeg, pointerToGoEnd, and sort.

Referenced by NamdState::loadStructure().

{
  DebugM(3,"->build_go_sigmas2" << std::endl);
  PDB *goPDB;      //  Pointer to PDB object to use
  int bcol = 4;      //  Column that data is in
  int32 chainType = 0;      //  b value from PDB file
  int32 residType = 0;      //  resid value from PDB file
  int i;         //  Loop counter
  int j;         //  Loop counter
  int resid1;    //  Residue ID for first atom
  int resid2;    //  Residue ID for second atom
  int residDiff;     //  Difference between resid1 and resid2
  Real sigma;    //  Sigma calculated for a Go atom pair
  Real atomAtomDist;     //  Distance between two atoms
  Real exp_a;            //  First exponent in L-J like formula
  Real exp_b;            //  Second exponent in L-J like formula
  char filename[NAMD_FILENAME_BUFFER_SIZE];    //  Filename
  
  //JLai
  int numLJPair = 0;
  long nativeContacts = 0;
  long nonnativeContacts = 0;

  //  Get the PDB object that contains the Go coordinates.  If
  //  the user gave another file name, use it.  Otherwise, just use
  //  the PDB file that has the initial coordinates.  
  if (goCoordFile == NULL)
    {
      //goPDB = initial_pdb;
      NAMD_die("Error: goCoordFile is NULL - build_go_sigmas2");
    }
  else
  {
    if (goCoordFile->next != NULL)
      {
        NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
      }
    
    if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
      {
        strcpy(filename, goCoordFile->data);
      }
    else
      {
        strcpy(filename, cwd);
        strcat(filename, goCoordFile->data);
      }
    
    goPDB = new PDB(filename);
    if ( goPDB == NULL )
      {
        NAMD_die("Memory allocation failed in Molecule::build_go_sigmas2");
      }
    
    if (goPDB->num_atoms() != numAtoms)
      {
        NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
      }
  }
  //  Allocate the array to hold the chain types
  atomChainTypes = new int32[numAtoms];
  //  Allocate the array to hold Go atom indices into the sigma array
  goSigmaIndices = new int32[numAtoms];
  //  Allocate the array to hold resid 
  goResidIndices = new int32[numAtoms];

  if (atomChainTypes == NULL) {
    NAMD_die("memory allocation failed in Molecule::build_go_sigmas2");
  }
  
  numGoAtoms = 0;
  
  //  Loop through all the atoms and get the Go chain types
  for (i=0; i<numAtoms; i++) {
    //  Get the chainType from the occupancy field
    chainType = (int32)(goPDB->atom(i))->occupancy();
    //  Get the resid from the resid field
    residType = (int32)(goPDB->atom(i))->residueseq();
    //  Assign the chainType value
    if ( chainType != 0 ) {
      //DebugM(3,"build_go_sigmas2 - atom:" << i << ", chainType:" << chainType << std::endl);
      atomChainTypes[i] = chainType;
      goSigmaIndices[i] = numGoAtoms;
      goResidIndices[i] = residType;
      numGoAtoms++;
    }
    else {
      atomChainTypes[i] = 0;
      goSigmaIndices[i] = -1;
      goResidIndices[i] = -1;
    }
  }

  //Define:
  ResizeArray<GoPair> tmpGoPair;
  
  //  Loop through all atom pairs and calculate sigmas
  // Store sigmas into sorted array
  DebugM(3,"    numAtoms=" << numAtoms << std::endl);
  for (i=0; i<numAtoms; i++) {
    resid1 = (goPDB->atom(i))->residueseq();
     for (j=i+1; j<numAtoms; j++) {
      resid2 = (goPDB->atom(j))->residueseq();
      if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 && !atoms_1to4(i,j) ) {
        residDiff = resid2 - resid1;
        if (residDiff < 0) residDiff = -residDiff;
        if ( atomChainTypes[i] && atomChainTypes[j] &&
             !(this->go_restricted(atomChainTypes[i],atomChainTypes[j],residDiff)) ) {
          atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
                              pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
                              pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
          if ( atomAtomDist <= this->get_go_cutoff(atomChainTypes[i],atomChainTypes[j]) ) {
            exp_a = this->get_go_exp_a(atomChainTypes[i],atomChainTypes[j]);
            exp_b = this->get_go_exp_b(atomChainTypes[i],atomChainTypes[j]);
            sigma = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
            double tmpA = pow(sigma,exp_a);
            double tmpB = pow(sigma,exp_b);
            GoPair gp;
            GoPair gp2;
            gp.goIndxA = i;
            gp.goIndxB = j;
            gp.A = tmpA;
            gp.B = tmpB;
            tmpGoPair.add(gp);
            gp2.goIndxA = j;
            gp2.goIndxB = i;
            gp2.A = tmpA;
            gp2.B = tmpB;
            tmpGoPair.add(gp2);
            nativeContacts++;
          } else {
            nonnativeContacts++;
          }
        } 
      } 
    }
  }

  iout << iINFO << "Number of UNIQUE    native contacts: " << nativeContacts << "\n" << endi;
  iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts << "\n" << endi;

  // Copies the resizeArray into a block of continuous memory
  std::sort(tmpGoPair.begin(),tmpGoPair.end(),goPairCompare);
  goNumLJPair = 2*nativeContacts;
  goIndxLJA = new int[goNumLJPair];
  goIndxLJB = new int[goNumLJPair];
  goSigmaPairA = new double[goNumLJPair];
  goSigmaPairB = new double[goNumLJPair];
  for(i=0; i< goNumLJPair; i++) {
    goIndxLJA[i] = tmpGoPair[i].goIndxA;
    goIndxLJB[i] = tmpGoPair[i].goIndxB;
    goSigmaPairA[i] = tmpGoPair[i].A;
    goSigmaPairB[i] = tmpGoPair[i].B;
  }

  pointerToGoBeg = new int[numAtoms];
  pointerToGoEnd = new int[numAtoms];
  int oldIndex = -1;
  for(i=0; i<numAtoms; i++) {
    pointerToGoBeg[i] = -1;
    pointerToGoEnd[i] = -2;
  }
  for(i=0; i < goNumLJPair; i++) {
    if(pointerToGoBeg[goIndxLJA[i]] == -1) {
      pointerToGoBeg[goIndxLJA[i]] = i;
      oldIndex = goIndxLJA[i];
    }
    pointerToGoEnd[oldIndex] = i;
  }

  //  If we had to create a PDB object, delete it now
  if (goCoordFile != NULL) {
    delete goPDB;
  }
    
  return;
}

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

Definition at line 6373 of file Molecule.C.

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

Referenced by NamdState::loadStructure().

{
    PDB *kPDB;
    register int i;             //  Loop counters
    register int j;
    register int k;

    DebugM(3,  "Entered build_gridforce_params multi...\n");
//     DebugM(3, "\tgridfrcfile = " << gridfrcfile->data << endi);
//     DebugM(3, "\tgridfrccol = " << gridfrccol->data << endi);
    
    MGridforceParams* mgridParams = simParams->mgridforcelist.get_first();
    numGridforceGrids = 0;
    while (mgridParams != NULL) {
        numGridforceGrids++;
        mgridParams = mgridParams->next;
    }
    
    DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n");
    gridfrcIndexes = new int32*[numGridforceGrids];
    gridfrcParams = new GridforceParams*[numGridforceGrids];
    gridfrcGrid = new GridforceGrid*[numGridforceGrids];
    numGridforces = new int[numGridforceGrids];
    
    int grandTotalGrids = 0;    // including all subgrids
    
    mgridParams = simParams->mgridforcelist.get_first();
    for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
        int current_index=0;    //  Index into values used
        int kcol = 5;           //  Column to look for force constant in
        int qcol = 0;           //  Column for charge (default 0: use electric charge)
        Real kval = 0;          //  Force constant value retreived
        char filename[NAMD_FILENAME_BUFFER_SIZE];       //  PDB filename
        char potfilename[NAMD_FILENAME_BUFFER_SIZE];    //  Potential file name
        
        if (mgridParams == NULL) {
            NAMD_die("Problem with mgridParams!");
        }
        
        // Now load values from mgridforcelist object
        if (mgridParams->gridforceFile == NULL)
        {
            if ( ! initial_pdb ) NAMD_die("Initial PDB file unavailable, gridforceFile required.");
            kPDB = initial_pdb;
        }
        else
        {
            DebugM(4, "mgridParams->gridforceFile = " << mgridParams->gridforceFile << "\n" << endi);
            
            if ( (cwd == NULL) || (mgridParams->gridforceFile[0] == '/') )
            {
                strcpy(filename, mgridParams->gridforceFile);
            }
            else
            {
                strcpy(filename, cwd);
                strcat(filename, mgridParams->gridforceFile);
            }
        
            kPDB = new PDB(filename);
            if ( kPDB == NULL )
            {
                NAMD_die("Memory allocation failed in Molecule::build_gridforce_params");
            }
           
            if (kPDB->num_atoms() != numAtoms)
            {
                NAMD_die("Number of atoms in grid force PDB doesn't match coordinate PDB");
            }
        }

        //  Get the column that the force constant is going to be in.  It
        //  can be in any of the 5 floating point fields in the PDB, according
        //  to what the user wants.  The allowable fields are X, Y, Z, O, or
        //  B which correspond to the 1st, 2nd, ... 5th floating point fields.
        //  The default is the 5th field, which is beta (temperature factor)
        if (mgridParams->gridforceCol == NULL)
        {
            kcol = 5;
        }
        else
        {
            if (strcasecmp(mgridParams->gridforceCol, "X") == 0)
            {
                kcol=1;
            }
            else if (strcasecmp(mgridParams->gridforceCol, "Y") == 0)
            {
                kcol=2;
            }
            else if (strcasecmp(mgridParams->gridforceCol, "Z") == 0)
            {
                kcol=3;
            }
            else if (strcasecmp(mgridParams->gridforceCol, "O") == 0)
            {
                kcol=4;
            }
            else if (strcasecmp(mgridParams->gridforceCol, "B") == 0)
            {
                kcol=5;
            }
            else
            {
                NAMD_die("gridforcecol must have value of X, Y, Z, O, or B");
            }
        }
    
        //  Get the column that the charge is going to be in.
        if (mgridParams->gridforceQcol == NULL)
        {
            qcol = 0;   // Default: don't read charge from file, use electric charge
        }
        else
        {
            if (strcasecmp(mgridParams->gridforceQcol, "X") == 0)
            {
                qcol=1;
            }
            else if (strcasecmp(mgridParams->gridforceQcol, "Y") == 0)
            {
                qcol=2;
            }
            else if (strcasecmp(mgridParams->gridforceQcol, "Z") == 0)
            {
                qcol=3;
            }
            else if (strcasecmp(mgridParams->gridforceQcol, "O") == 0)
            {
                qcol=4;
            }
            else if (strcasecmp(mgridParams->gridforceQcol, "B") == 0)
            {
                qcol=5;
            }
            else
            {
                NAMD_die("gridforcechargecol must have value of X, Y, Z, O, or B");
            }
        }
    
        if (kcol == qcol) {
            NAMD_die("gridforcecol and gridforcechargecol cannot have same value");
        }

    
        //  Allocate an array that will store an index into the constraint
        //  parameters for each atom.  If the atom is not constrained, its
        //  value will be set to -1 in this array.
        gridfrcIndexes[gridnum] = new int32[numAtoms];
       
        if (gridfrcIndexes[gridnum] == NULL)
        {
            NAMD_die("memory allocation failed in Molecule::build_gridforce_params()");
        }
        
        //  Loop through all the atoms and find out which ones are constrained
        for (i=0; i<numAtoms; i++)
        {
            //  Get the k value based on where we were told to find it
            switch (kcol)
            {
            case 1:
                kval = (kPDB->atom(i))->xcoor();
                break;
            case 2:
                kval = (kPDB->atom(i))->ycoor();
                break;
            case 3:
                kval = (kPDB->atom(i))->zcoor();
                break;
            case 4:
                kval = (kPDB->atom(i))->occupancy();
                break;
            case 5:
                kval = (kPDB->atom(i))->temperaturefactor();
                break;
            }
           
            if (kval > 0.0)
            {
                //  This atom is constrained
                gridfrcIndexes[gridnum][i] = current_index;
                current_index++;
            }
            else
            {
                //  This atom is not constrained
                gridfrcIndexes[gridnum][i] = -1;
            }
        }
    
        if (current_index == 0)
        {
            //  Constraints were turned on, but there weren't really any constrained
            iout << iWARN << "NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" << endi;
        }
        else
        {
            //  Allocate an array to hold the constraint parameters
            gridfrcParams[gridnum] = new GridforceParams[current_index];
            if (gridfrcParams[gridnum] == NULL)
            {
                NAMD_die("memory allocation failed in Molecule::build_gridforce_params");
            }
        }
    
        numGridforces[gridnum] = current_index;

        //  Loop through all the atoms and assign the parameters for those
        //  that are constrained
        for (i=0; i<numAtoms; i++)
        {
            if (gridfrcIndexes[gridnum][i] != -1)
            {
                //  This atom has grid force, so get the k value again
                switch (kcol)
                {
                case 1:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->xcoor();
                    break;
                case 2:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->ycoor();
                    break;
                case 3:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->zcoor();
                    break;
                case 4:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->occupancy();
                    break;
                case 5:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->temperaturefactor();
                    break;
                }
            
                //  Also get charge column
                switch (qcol)
                {
                case 0:
#ifdef MEM_OPT_VERSION
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]];
#else
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge;
#endif
                    break;
                case 1:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->xcoor();
                    break;
                case 2:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->ycoor();
                    break;
                case 3:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->zcoor();
                    break;
                case 4:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->occupancy();
                    break;
                case 5:
                    gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->temperaturefactor();
                    break;
                }
            }
        }
       
        //  If we had to create new PDB objects, delete them now
        if (mgridParams->gridforceFile != NULL)
        {
            delete kPDB;
        }
    
        //  Now we fill in our grid information
    
        // Open potential file
        if ( (cwd == NULL) || (mgridParams->gridforceVfile[0] == '/') )
        {
            strcpy(potfilename, mgridParams->gridforceVfile);
        }
        else
        {
            strcpy(potfilename, cwd);
            strcat(potfilename, mgridParams->gridforceVfile);
        }
    
//        iout << iINFO << "Allocating grid " << gridnum
//             << "\n" << endi;
        
        DebugM(3, "allocating GridforceGrid(" << gridnum << ")\n");
        gridfrcGrid[gridnum] = GridforceGrid::new_grid(gridnum, potfilename, simParams, mgridParams);
        
        grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
        DebugM(4, "grandTotalGrids = " << grandTotalGrids << "\n" << endi);
        
        // Finally, get next mgridParams pointer
        mgridParams = mgridParams->next;
    }
}

void Molecule::build_gro_pair

(

)

void Molecule::build_langevin_params	(	BigReal	coupling,
		BigReal	drudeCoupling,
		Bool	doHydrogen
	)

Referenced by NamdState::loadStructure().

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

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

Referenced by NamdState::loadStructure().

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

Referenced by NamdState::loadStructure().

void Molecule::build_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 *
	)

Referenced by NamdState::loadStructure().

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

void Molecule::compute_LJcorrection

(

)

Referenced by NamdState::loadStructure().

void Molecule::delete_alch_bonded

(

void

)

Referenced by NamdState::loadStructure().

void Molecule::delete_qm_bonded

(

void

)

Definition at line 1547 of file MoleculeQM.C.

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

                                     {

#ifdef MEM_OPT_VERSION
    NAMD_die("QMMM interface is not supported in memory optimized builds.");
#else
    
    DebugM(3,"Cleaning QM bonds, angles, dihedrals, impropers and crossterms.\n")
    
    // Bonds
    Bond *nonQMBonds;
    nonQMBonds = new Bond[numBonds] ;
    int nonQMBondCount = 0;
    qmDroppedBonds = 0;
    for (int i = 0; i < numBonds; i++) {
        
        int part1 = qmAtomGroup[bonds[i].atom1];
        int part2 = qmAtomGroup[bonds[i].atom2];
        if (part1 > 0 && part2 > 0 ) {
            
            // If the user defined extra bonds, we check if the QM-QM bond is an "extra"
            // bond, and do not delete it.
            if (simParams->extraBondsOn) {
//                 std::cout << "Checking Bond: " << bonds[i].atom1 << "," << bonds[i].atom2 << std::endl ;
                
                for (int ebi=0; ebi < qmExtraBonds.size() ; ebi++) {
                    
                    if( (qmExtraBonds[ebi].atom1 == bonds[i].atom1 || 
                        qmExtraBonds[ebi].atom1 == bonds[i].atom2) && 
                    (qmExtraBonds[ebi].atom2 == bonds[i].atom1 || 
                        qmExtraBonds[ebi].atom2 == bonds[i].atom2) ) {
//                         std::cout << "This is an extra bond! We will keep it." << std::endl;
                        nonQMBonds[nonQMBondCount++] = bonds[i];
                        break;
                    }
                }
            }
            
            qmDroppedBonds++;
        } else {
            // Just a simple sanity check.
            // If there are no QM bonds defined by the user but we find a
            // bond between a QM and an MM atom, error out.
            if ((part1 > 0 || part2 > 0) && qmNumBonds == 0) {
                iout << iERROR << " Atoms " << bonds[i].atom1
                << " and " << bonds[i].atom2 << " form a QM-MM bond!\n" << endi ;
                NAMD_die("No QM-MM bonds were defined by the user, but one was found.");
            }
            nonQMBonds[nonQMBondCount++] = bonds[i];
        }
    }
    numBonds = nonQMBondCount;
    delete [] bonds;
    bonds = new Bond[numBonds];
    for (int i = 0; i < nonQMBondCount; i++) {
        bonds[i]=nonQMBonds[i];
    }
    delete [] nonQMBonds;
  numRealBonds = numBonds;
    
  // Angles
  Angle *nonQMAngles;
  nonQMAngles = new Angle[numAngles];
  int nonQMAngleCount = 0;
  qmDroppedAngles = 0;
  for (int i = 0; i < numAngles; i++) {
    int part1 = qmAtomGroup[angles[i].atom1];
    int part2 = qmAtomGroup[angles[i].atom2];
    int part3 = qmAtomGroup[angles[i].atom3];
    if (part1 > 0 && part2 > 0 && part3 > 0) {
      //CkPrintf("-----ANGLE ATOMS %i %i %i partitions %i %i %i\n",angles[i].atom1, angles[i].atom2, angles[i].atom3, part1, part2, part3);
      qmDroppedAngles++;
    }
    else {
      nonQMAngles[nonQMAngleCount++] = angles[i];
    }
  }
  numAngles = nonQMAngleCount;
  delete [] angles;
  angles = new Angle[numAngles];
  for (int i = 0; i < nonQMAngleCount; i++) {
    angles[i]=nonQMAngles[i];
  }
  delete [] nonQMAngles;


  // Dihedrals
  Dihedral *nonQMDihedrals;
  nonQMDihedrals = new Dihedral[numDihedrals];
  int nonQMDihedralCount = 0;
  qmDroppedDihedrals = 0;
  for (int i = 0; i < numDihedrals; i++) {
    int part1 = qmAtomGroup[dihedrals[i].atom1];
    int part2 = qmAtomGroup[dihedrals[i].atom2];
    int part3 = qmAtomGroup[dihedrals[i].atom3];
    int part4 = qmAtomGroup[dihedrals[i].atom4];
    if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0) {
      //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);
      qmDroppedDihedrals++;
    }
    else {
      nonQMDihedrals[nonQMDihedralCount++] = dihedrals[i];
    }
  }
  numDihedrals = nonQMDihedralCount;
  delete [] dihedrals;
  dihedrals = new Dihedral[numDihedrals];
  for (int i = 0; i < numDihedrals; i++) {
    dihedrals[i]=nonQMDihedrals[i];
  }
  delete [] nonQMDihedrals;

  // Impropers
  Improper *nonQMImpropers;
  nonQMImpropers = new Improper[numImpropers];
  int nonQMImproperCount = 0;
  qmDroppedImpropers = 0;
  for (int i = 0; i < numImpropers; i++) {
    int part1 = qmAtomGroup[impropers[i].atom1];
    int part2 = qmAtomGroup[impropers[i].atom2];
    int part3 = qmAtomGroup[impropers[i].atom3];
    int part4 = qmAtomGroup[impropers[i].atom4];
    if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0) {
      //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);
      qmDroppedImpropers++;
    }
    else {
      nonQMImpropers[nonQMImproperCount++] = impropers[i];
    }
  }
  numImpropers = nonQMImproperCount;
  delete [] impropers;
  impropers = new Improper[numImpropers];
  for (int i = 0; i < numImpropers; i++) {
    impropers[i]=nonQMImpropers[i];
  }
  delete [] nonQMImpropers;
  
  // Crossterms 
  Crossterm *nonQMCrossterms;
  nonQMCrossterms = new Crossterm[numCrossterms];
  int nonQMCrosstermCount = 0;
  qmDroppedCrossterms = 0 ;
  for (int i = 0; i < numCrossterms; i++) {
    int part1 = qmAtomGroup[crossterms[i].atom1];
    int part2 = qmAtomGroup[crossterms[i].atom2];
    int part3 = qmAtomGroup[crossterms[i].atom3];
    int part4 = qmAtomGroup[crossterms[i].atom4];
    int part5 = qmAtomGroup[crossterms[i].atom5];
    int part6 = qmAtomGroup[crossterms[i].atom6];
    int part7 = qmAtomGroup[crossterms[i].atom7];
    int part8 = qmAtomGroup[crossterms[i].atom8];
    if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0 &&
        part5 > 0 && part6 > 0 && part7 > 0 && part8 > 0) {
      //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);
      qmDroppedCrossterms++;
    }
    else {
      nonQMCrossterms[nonQMCrosstermCount++] = crossterms[i];
    }
  }
  numCrossterms = nonQMCrosstermCount;
  delete [] crossterms;
  crossterms = new Crossterm[numCrossterms] ;
  for (int i=0; i < numCrossterms; i++){
      crossterms[i] = nonQMCrossterms[i] ;
  }
  delete [] nonQMCrossterms ;
  
  if (!CkMyPe()) {
      iout << iINFO << "The QM region will remove " << qmDroppedBonds << " bonds, " << 
          qmDroppedAngles << " angles, " << qmDroppedDihedrals << " dihedrals, "
          << qmDroppedImpropers << " impropers and " << qmDroppedCrossterms 
          << " crossterms.\n" << endi ;
  }
  
#endif
}

void Molecule::freeBFactorData ( )

inline

Definition at line 1038 of file Molecule.h.

Referenced by NamdState::loadStructure().

{ delete [] bfactor; bfactor=NULL; }

void Molecule::freeOccupancyData ( )

inline

Definition at line 1034 of file Molecule.h.

Referenced by NamdState::loadStructure().

{ delete [] occupancy; occupancy=NULL; }

Bond* Molecule::get_acceptor ( int  dnum ) const

inline

Definition at line 1109 of file Molecule.h.

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

Angle* Molecule::get_angle ( int  anum ) const

inline

Definition at line 1072 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

int32* Molecule::get_angles_for_atom ( int  anum )

inline

Definition at line 1151 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

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

                                                       {

  if (atomNames == NULL || resLookup == NULL)
  {
    NAMD_die("Tried to find atom from name on node other than node 0");
  }
  int i = 0;
  int end = 0;
  if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
  if ( index >= 0 && index < ( end - i ) ) return ( index + i );
  return -1;
}

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

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

                                                               {

  if (atomNames == NULL || resLookup == NULL)
  {
    NAMD_die("Tried to find atom from name on node other than node 0");
  }

  int i = 0;
  int end = 0;
  if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
  for ( ; i < end; ++i ) {
    #ifdef MEM_OPT_VERSION    
    Index idx = atomNames[i].atomnameIdx;
    if(!strcasecmp(aname, atomNamePool[idx])) return i;
    #else
    if ( ! strcasecmp(aname,atomNames[i].atomname) ) return i;
    #endif
  }
  return -1;
}

const char* Molecule::get_atomtype ( int  anum ) const

inline

Definition at line 1120 of file Molecule.h.

References atomTypePool, and NAMD_die().

  {
    if (atomNames == NULL)
    {
      NAMD_die("Tried to find atom type on node other than node 0");
    }

    #ifdef MEM_OPT_VERSION    
    return atomTypePool[atomNames[anum].atomtypeIdx];
    #else
    return(atomNames[anum].atomtype);
    #endif
  }

Bond* Molecule::get_bond ( int  bnum ) const

inline

Definition at line 1069 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

int32* Molecule::get_bonds_for_atom ( int  anum )

inline

Definition at line 1149 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

      { return bondsByAtom[anum]; } 

int Molecule::get_cluster ( int  anum ) const

inline

Definition at line 1027 of file Molecule.h.

Referenced by wrap_coor_int().

{ return cluster[anum]; }

int Molecule::get_clusterSize ( int  anum ) const

inline

Definition at line 1028 of file Molecule.h.

Referenced by wrap_coor_int().

{ return clusterSize[anum]; }

void Molecule::get_cons_params ( Real &  k, Vector &  refPos, int  atomnum  ) const

inline

Definition at line 1266 of file Molecule.h.

Referenced by ComputeRestraints::doForce().

  {
    k = consParams[consIndexes[atomnum]].k;
    refPos = consParams[consIndexes[atomnum]].refPos;
  }

void Molecule::get_constorque_params ( BigReal &  v, Vector &  a, Vector &  p, int  atomnum  ) const

inline

Definition at line 1340 of file Molecule.h.

References consTorqueIndexes, and consTorqueParams.

Referenced by ComputeConsTorque::doForce().

  {
    v = consTorqueParams[consTorqueIndexes[atomnum]].v;
    a = consTorqueParams[consTorqueIndexes[atomnum]].a;
    p = consTorqueParams[consTorqueIndexes[atomnum]].p;
  }

Crossterm* Molecule::get_crossterm ( int  inum ) const

inline

Definition at line 1081 of file Molecule.h.

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

int32* Molecule::get_crossterms_for_atom ( int  anum )

inline

Definition at line 1157 of file Molecule.h.

      { return crosstermsByAtom[anum]; }  

const Real* Molecule::get_cSMDcoffs ( )

inline

Definition at line 875 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return cSMDcoffs;} ;

int const* const* Molecule::get_cSMDindex ( )

inline

Definition at line 871 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return cSMDindex;} ;

int const* Molecule::get_cSMDindxLen ( )

inline

Definition at line 870 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return cSMDindxLen;} ;

const Real* Molecule::get_cSMDKs ( )

inline

Definition at line 873 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return cSMDKs;} ;

int Molecule::get_cSMDnumInst ( )

inline

Definition at line 869 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return cSMDnumInst;} ;

int const* const* Molecule::get_cSMDpairs ( )

inline

Definition at line 872 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return cSMDpairs;} ;

const Real* Molecule::get_cSMDVels ( )

inline

Definition at line 874 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return cSMDVels;} ;

Dihedral* Molecule::get_dihedral ( int  dnum ) const

inline

Definition at line 1078 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

int32* Molecule::get_dihedrals_for_atom ( int  anum )

inline

Definition at line 1153 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

      { return dihedralsByAtom[anum]; }

Bond* Molecule::get_donor ( int  dnum ) const

inline

Definition at line 1106 of file Molecule.h.

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

const ExclusionCheck* Molecule::get_excl_check_for_atom ( int  anum ) const

inline

Definition at line 1177 of file Molecule.h.

Referenced by dumpbench().

                                                               {      
      return &all_exclusions[anum];             
  }

Exclusion* Molecule::get_exclusion ( int  ex ) const

inline

Definition at line 1116 of file Molecule.h.

References exclusions.

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

int32* Molecule::get_exclusions_for_atom ( int  anum )

inline

Definition at line 1159 of file Molecule.h.

References exclusionsByAtom.

      { return exclusionsByAtom[anum]; }

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

inline

Definition at line 1389 of file Molecule.h.

References NAMD_die(), order, and simParams.

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

  {
    int typeSum = 0;
    for ( int i=0; i < order; ++i ) {
      typeSum += (fepAtomFlags[atomID[i]] == 2 ? -1 : fepAtomFlags[atomID[i]]);
    }
    // Increase the cutoff if scaling purely alchemical bonds. 
    // This really only applies when order = 2.
    if ( simParams->alchBondDecouple ) order++;

    // The majority of interactions are type 0, so catch those first.
    if ( typeSum == 0 || abs(typeSum) == order ) return 0;
    else if ( 0 < typeSum && typeSum < order ) return 1;
    else if ( -order < typeSum && typeSum < 0 ) return 2;

    // Alchemify should always keep this from bombing, but just in case...
    NAMD_die("Unexpected alchemical bonded interaction!");
    return 0;
  }

unsigned char Molecule::get_fep_type ( int  anum ) const

inline

Definition at line 1349 of file Molecule.h.

Referenced by CrosstermElem::computeForce(), TholeElem::computeForce(), WorkDistrib::createAtomLists(), and ComputeHomeTuples< TholeElem, Thole, TholeValue >::loadTuples().

  {
    return(fepAtomFlags[anum]);
  }

const int32* Molecule::get_full_exclusions_for_atom ( int  anum ) const

inline

Definition at line 1161 of file Molecule.h.

Referenced by ComputeNonbondedCUDA::build_exclusions().

      { return fullExclusionsByAtom[anum]; }

Real Molecule::get_go_cutoff ( int  chain1, int  chain2  )

inline

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

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

Real Molecule::get_go_epsilon ( int  chain1, int  chain2  )

inline

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

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

Real Molecule::get_go_epsilonRep ( int  chain1, int  chain2  )

inline

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

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

int Molecule::get_go_exp_a ( int  chain1, int  chain2  )

inline

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

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

int Molecule::get_go_exp_b ( int  chain1, int  chain2  )

inline

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

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

int Molecule::get_go_exp_rep ( int  chain1, int  chain2  )

inline

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

{ 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(), goForce, goSigmaIndices, goSigmas, goWithinCutoff, and numGoAtoms.

{
  BigReal goForce = 0.0;
  Real pow1;
  Real pow2;
  //  determine which Go chain pair we are working with
  //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ")" << std::endl);
  int32 chain1 = atomChainTypes[atom1];
  int32 chain2 = atomChainTypes[atom2];

  //DebugM(3,"  chain1:" << chain1 << ", chain2:" << chain2 << std::endl);
  if (chain1 == 0 || chain2 == 0)  return 0.0;

  //  retrieve Go cutoff for this chain pair
  //TMP// JLai -- I'm going to replace this with a constant accessor.  This is just a temporary thing
  Real goCutoff = const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2);
  //DebugM(3,"  goCutoff:" << goCutoff << std::endl);
  if (goCutoff == 0)  return 0.0;
  //  if repulsive then calculate repulsive
  //  sigmas are initially set to -1.0 if the atom pair fails go_restricted
  if (goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]] != -1.0) {
    if (!goWithinCutoff[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]]) {
      Real epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1,chain2);
      Real sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1,chain2);
      int exp_rep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
      pow1 = pow(sigmaRep/r,exp_rep);
      goForce = 4*((exp_rep/(r*r)) * epsilonRep * pow1);
      *goNative = 0.0;
      *goNonnative = (4 * epsilonRep * pow1 );
      //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", epsilonRep:" << epsilonRep << ", sigmaRep:" << sigmaRep << ", r:" << r << ", goForce:" << goForce << std::endl);
    }
    //  if attractive then calculate attractive
    else {
      int goSigmaIndex1 = goSigmaIndices[atom1];
      int goSigmaIndex2 = goSigmaIndices[atom2];
      if (goSigmaIndex1 != -1 && goSigmaIndex2 != -1) {
        Real epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
        int exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
        int exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
        Real sigma_ij = goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]];
        // Positive gradient of potential, not negative gradient of potential
        pow1 = pow(sigma_ij/r,exp_a);
        pow2 = pow(sigma_ij/r,exp_b);
        goForce = ((4/(r*r)) * epsilon * (exp_a * pow1 - exp_b * pow2));
        //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", sigma_ij:" << sigma_ij << ", r:" << r << ", goForce:" << goForce << std::endl);
        *goNative = (4 * epsilon * ( pow1 -  pow2 ) );
        *goNonnative = 0.0;
      }
    }
  }
  //DebugM(3,"goForce:" << goForce << std::endl);
  return goForce;
}

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, pointerToGoEnd, and z.

{
 
  // Check to see if restricted.  If so, escape function early
  int32 chain1 = atomChainTypes[atom1];
  int32 chain2 = atomChainTypes[atom2];

  if(chain1 == 0 || chain2 == 0) return 0.0;
  Molecule *mol = const_cast<Molecule*>(this);
  Real goCutoff = mol->get_go_cutoff(chain1,chain2);
  if(goCutoff == 0) return 0.0;

  int resid1 = goResidIndices[atom1];
  int resid2 = goResidIndices[atom2];
  int residDiff = abs(resid1 - resid2);
  if((mol->go_restricted(chain1,chain2,residDiff))) {
    return 0.0;
  }

  int LJIndex = -1;
  int LJbegin = pointerToGoBeg[atom1];
  int LJend = pointerToGoEnd[atom1];
  for(int i = LJbegin; i <= LJend; i++) {
    if(goIndxLJB[i] == atom2) {
      LJIndex = i;
    }
  }
  
  BigReal r2 = x*x + y*y + z*z;
  BigReal r = sqrt(r2);

  if (LJIndex == -1) {
    int exp_rep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
    BigReal epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1, chain2);
    BigReal sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1, chain2);
    double sigmaRepPow = pow(sigmaRep,exp_rep);
    BigReal LJ = (4*epsilonRep*exp_rep*sigmaRepPow*pow(r,-(exp_rep+1)));
    *goNative = 0;
    *goNonnative = (4*epsilonRep*sigmaRepPow*pow(r,-exp_rep));
    //*goNonnative = (4*epsilonRep * pow(sigmaRep/r,exp_rep));
    return (LJ/r);
  } else {
    // Code to calculate distances because the pair was found in one of the lists
    int exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
    int exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
    // We want the force, so we have to take the n+1 derivative
    BigReal powA = pow(r,-(exp_a + 1));
    BigReal powB = pow(r,-(exp_b + 1));
    BigReal powaa = pow(r,-exp_a);
    BigReal powbb = pow(r,-exp_b);
    BigReal epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
    BigReal LJ = 4 * epsilon * (exp_a*goSigmaPairA[LJIndex]*powA - exp_b*goSigmaPairB[LJIndex]*powB);
    *goNative =  4 * epsilon * (goSigmaPairA[LJIndex]*powaa - goSigmaPairB[LJIndex]*powbb);
    *goNonnative = 0;
    return (LJ/r);
  }
  return 0;
}

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, goForce, goResids, and goSigmaIndices.

{
  int resid1;
  int resid2;
  int residDiff;
  Real xcoorDiff;
  Real ycoorDiff;
  Real zcoorDiff;
  Real atomAtomDist;
  Real exp_a;
  Real exp_b;
  Real sigma_ij;
  Real epsilon;
  Real epsilonRep;
  Real sigmaRep;
  Real expRep;
  Real pow1;
  Real pow2;
  
  BigReal goForce = 0.0;
  *goNative = 0;
  *goNonnative = 0;

  //  determine which Go chain pair we are working with
  DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ")" << std::endl);
  int goIndex1 = goSigmaIndices[atom1];
  int goIndex2 = goSigmaIndices[atom2];

  int32 chain1 = atomChainTypes[goIndex1];
  int32 chain2 = atomChainTypes[goIndex2];

  DebugM(3,"  chain1:" << chain1 << ", chain2:" << chain2 << std::endl);
  if (chain1 == 0 || chain2 == 0)  return 0.0;

  //  retrieve Go cutoff for this chain pair
  Real goCutoff = const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2);
  DebugM(3,"  goCutoff:" << goCutoff << std::endl);
  if (goCutoff == 0)  return 0.0;

  //  sigmas are initially set to -1.0 if the atom pair fails go_restricted
  //  no goSigmas array anymore
  //Real sigma_ij = goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]];

  // XXX - used to be a condition for the following if
  //if the atoms are within 4 of each other
  //if ( !atoms_1to4(atom1,atom2) ) {

  //  if goSigmaIndices aren't defined, don't calculate forces
  if ( goIndex1 != -1 && goIndex2 != -1 ) {
    resid1 = goResids[goIndex1];
    resid2 = goResids[goIndex2];
    residDiff = resid2 - resid1;
    if (residDiff < 0) residDiff = -residDiff;
    //  if this is a Go atom pair that is not restricted
    if ( !(const_cast<Molecule*>(this)->go_restricted(chain1,chain2,residDiff)) ) {
      xcoorDiff = goCoordinates[goIndex1*3] - goCoordinates[goIndex2*3];
      ycoorDiff = goCoordinates[goIndex1*3 + 1] - goCoordinates[goIndex2*3 + 1];
      zcoorDiff = goCoordinates[goIndex1*3 + 2] - goCoordinates[goIndex2*3 + 2];
      atomAtomDist = sqrt(xcoorDiff*xcoorDiff + ycoorDiff*ycoorDiff + zcoorDiff*zcoorDiff);
      
      //  if attractive then calculate attractive
      if ( atomAtomDist <= const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2) ) {
        exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
        exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
        sigma_ij = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
        
        // [JLai] print out atoms involved in native contacts
        // printf("ATOM1: %d C1: %d ATOM2: %d C2: %d\n", atom1,chain1,atom2,chain2);

        epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
        pow1 = pow(sigma_ij/r,static_cast<double>(exp_a));
        pow2 = pow(sigma_ij/r,static_cast<double>(exp_b));
        //goForce = ((4/r) * epsilon * (exp_a * pow(sigma_ij/r,exp_a) - exp_b * pow(sigma_ij/r,exp_b)));
        goForce = ((4/(r*r)) * epsilon * (exp_a * pow1 - exp_b * pow2));
        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);
        //goEnergy = (4 * epsilon * ( pow(sigma_ij/r,exp_a) -  pow(sigma_ij/r,exp_b) ) ); // JLai I changed some of the expressions
        *goNative = (4 * epsilon * ( pow1 -  pow2 ) ); 
        *goNonnative = 0;
      }
      
      //  if repulsive then calculate repulsive
      else {
        epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1,chain2);
        sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1,chain2);
        expRep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
        pow1 = pow(sigmaRep/r,(BigReal)expRep);
        //goForce = ((12.0/r) * epsilonRep * pow(sigmaRep/r,12.0));
        goForce = ((4/(r*r)) * expRep * epsilonRep * pow1);
        DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", epsilonRep:" << epsilonRep << ", sigmaRep:" << sigmaRep << ", r:" << r << ", goForce:" << goForce << std::endl);
        //goEnergy = (4 * epsilonRep * pow(sigmaRep/r,12.0)); // JLai I changed some of the expressions
        *goNonnative = (4 * epsilonRep * pow1); 
        *goNative = 0;
      }
    }
  }
  
  //DebugM(3,"goForce:" << goForce << std::endl);
  return goForce;
}

Real Molecule::get_go_sigmaRep ( int  chain1, int  chain2  )

inline

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

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

GridforceGrid* Molecule::get_gridfrc_grid ( int  gridnum ) const

inline

Definition at line 1280 of file Molecule.h.

References numGridforceGrids.

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

  {
      GridforceGrid *result = NULL;
      if (gridnum >= 0 && gridnum < numGridforceGrids) {
          result = gridfrcGrid[gridnum];
      }
      return result;
  }

void Molecule::get_gridfrc_params ( Real &  k, Charge &  q, int  atomnum, int  gridnum  ) const

inline

Definition at line 1274 of file Molecule.h.

Referenced by ComputeGridForce::do_calc().

  {
      k = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].k;
      q = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].q;
  }

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, pointerToLJEnd, and z.

{
  //Initialize return energies to zero
  *pairLJEnergy = 0.0;
  *pairGaussEnergy = 0.0;

  // Linear search for LJ data
  int LJIndex = -1;
  int LJbegin = pointerToLJBeg[atom1];
  int LJend = pointerToLJEnd[atom1];
  for(int i = LJbegin; i <= LJend; i++) {
    if(indxLJB[i] == atom2) {
      LJIndex = i;
      break;
    }
  }

  // Linear search for Gaussian data
  int GaussIndex = -1;
  int Gaussbegin = pointerToGaussBeg[atom1];
  int Gaussend = pointerToGaussEnd[atom1];
  for(int i = Gaussbegin; i <= Gaussend; i++) {
    if(indxGaussB[i] == atom2) {
      GaussIndex = i;
      break;
    }
  }

  if( LJIndex == -1 && GaussIndex == -1) {
    return 0;
  } else {
    // Code to calculate distances because the pair was found in one of the lists
    BigReal r2 = x*x + y*y + z*z;
    BigReal r = sqrt(r2);
    BigReal ri = 1/r;
    BigReal ri6 = (ri*ri*ri*ri*ri*ri);
    BigReal ri12 = ri6*ri6;
    BigReal ri13 = ri12*ri;
    BigReal LJ = 0;
    BigReal Gauss = 0;
    // Code to calculate LJ 
    if (LJIndex != -1) {
      BigReal ri7 = ri6*ri;
      LJ = (12*(pairC12[LJIndex]*ri13) - 6*(pairC6[LJIndex]*ri7));
      *pairLJEnergy = (pairC12[LJIndex]*ri12 - pairC6[LJIndex]*ri6);
      //std::cout << pairC12[LJIndex] << " " << pairC6[LJIndex] << " " << ri13 << " " << ri7 << " " << LJ << " " << r << "\n";
    }
    // Code to calculate Gaussian
    if (GaussIndex != -1) {
      BigReal gr = 12*gRepulsive[GaussIndex]*ri13;
      BigReal r1prime = r - gMu1[GaussIndex];
      BigReal tmp1 = r1prime * r1prime;
      BigReal r2prime = r - gMu2[GaussIndex];
      BigReal tmp2 = r2prime * r2prime;
      BigReal tmp_gauss1 = 0;
      BigReal one_gauss1 = 1;
      BigReal tmp_gauss2 = 0;
      BigReal one_gauss2 = 1;
      if (giSigma1[GaussIndex] != 0) {
        tmp_gauss1 = exp(-tmp1*giSigma1[GaussIndex]);
        one_gauss1 = 1 - tmp_gauss1;
      }
      if (giSigma2[GaussIndex] != 0) {
        tmp_gauss2 = exp(-tmp2*giSigma2[GaussIndex]);
        one_gauss2 = 1 - tmp_gauss2;
      } 
      BigReal A = gA[GaussIndex];
      Gauss = gr*one_gauss1*one_gauss2 - 2*A*tmp_gauss1*one_gauss2*r1prime*giSigma1[GaussIndex] \
        - 2*tmp_gauss1*one_gauss2*r1prime*giSigma1[GaussIndex]*gRepulsive[GaussIndex]*ri12 - 2*A*tmp_gauss2*one_gauss1*r2prime*giSigma2[GaussIndex] \
        - 2*tmp_gauss2*one_gauss1*r2prime*giSigma2[GaussIndex]*gRepulsive[GaussIndex]*ri12;
      *pairGaussEnergy = A*(-1+(one_gauss1)*(one_gauss2)*(1+gRepulsive[GaussIndex]*ri12/A));
    }
    //std::cout << "Net force: " << (LJ + Gauss) << " with ri " << (LJ + Gauss)*ri << "\n";
    return (LJ + Gauss)*ri;
  }
  return 0;
}

int Molecule::get_groupSize

(

int

)

Improper* Molecule::get_improper ( int  inum ) const

inline

Definition at line 1075 of file Molecule.h.

Referenced by dumpbench().

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

int32* Molecule::get_impropers_for_atom ( int  anum )

inline

Definition at line 1155 of file Molecule.h.

Referenced by dumpbench().

      { return impropersByAtom[anum]; }  

Lphost* Molecule::get_lphost ( int  atomid ) const

inline

Definition at line 1085 of file Molecule.h.

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

const int32* Molecule::get_mod_exclusions_for_atom ( int  anum ) const

inline

Definition at line 1163 of file Molecule.h.

      { return modExclusionsByAtom[anum]; }

int Molecule::get_mother_atom

(

int

)

void Molecule::get_movdrag_params ( Vector &  v, int  atomnum  ) const

inline

Definition at line 1325 of file Molecule.h.

Referenced by Sequencer::addMovDragToPosition().

  {
    v = movDragParams[movDragIndexes[atomnum]].v;
  }

Bool Molecule::get_noPC ( )

inline

Definition at line 857 of file Molecule.h.

Referenced by ComputeQM::initialize().

{ return qmNoPC; } ;

int Molecule::get_numQMAtoms ( )

inline

Definition at line 826 of file Molecule.h.

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

{return qmNumQMAtoms; } ;

Real* Molecule::get_qmAtmChrg ( )

inline

Definition at line 823 of file Molecule.h.

Referenced by ComputeQMMgr::calcMOPAC(), ComputeQMMgr::calcORCA(), ComputeQMMgr::calcUSR(), ComputePme::doQMWork(), ComputeQM::initialize(), HomePatch::qmSwapAtoms(), and ComputeQM::saveResults().

{return qmAtmChrg; } ;

const int* Molecule::get_qmAtmIndx ( )

inline

Definition at line 824 of file Molecule.h.

Referenced by ComputeQMMgr::calcMOPAC(), ComputeQMMgr::calcORCA(), ComputeQMMgr::calcUSR(), ComputePme::doQMWork(), ComputeQM::initialize(), HomePatch::qmSwapAtoms(), and ComputeQM::saveResults().

{return qmAtmIndx; } ;

const Real* Molecule::get_qmAtomGroup ( ) const

inline

Definition at line 820 of file Molecule.h.

Referenced by ComputeQMMgr::calcMOPAC(), ComputeQMMgr::calcORCA(), ComputeQMMgr::calcUSR(), ComputePme::doQMWork(), ComputeQM::initialize(), HomePatch::qmSwapAtoms(), and ComputeQM::saveResults().

{return qmAtomGroup; } ;

Real Molecule::get_qmAtomGroup ( int  indx ) const

inline

Definition at line 821 of file Molecule.h.

{return qmAtomGroup[indx]; } ;

Bool Molecule::get_qmcSMD ( )

inline

Definition at line 868 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmcSMD;} ;

int* Molecule::get_qmCustomPCIdxs ( )

inline

Definition at line 866 of file Molecule.h.

Referenced by ComputeQM::initialize().

{ return qmCustomPCIdxs; } ;

int* Molecule::get_qmCustPCSizes ( )

inline

Definition at line 865 of file Molecule.h.

Referenced by ComputeQM::initialize().

{ return qmCustPCSizes; } ;

const BigReal* Molecule::get_qmDummyBondVal ( )

inline

Definition at line 835 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmDummyBondVal; } ;

const char* const* Molecule::get_qmDummyElement ( )

inline

Definition at line 840 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmDummyElement; } ;

const char* const* Molecule::get_qmElements ( )

inline

Definition at line 827 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{return qmElementArray;} ;

const int* const* Molecule::get_qmGrpBonds ( )

inline

Definition at line 838 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmGrpBonds; } ;

Real* Molecule::get_qmGrpChrg ( )

inline

Definition at line 831 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return qmGrpChrg; } ;

Real* Molecule::get_qmGrpID ( )

inline

Definition at line 830 of file Molecule.h.

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

{ return qmGrpID; } ;

Real* Molecule::get_qmGrpMult ( )

inline

Definition at line 832 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return qmGrpMult; } ;

const int* Molecule::get_qmGrpNumBonds ( )

inline

Definition at line 836 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmGrpNumBonds; } ;

const int* Molecule::get_qmGrpSizes ( )

inline

Definition at line 829 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{return qmGrpSizes; } ;

int Molecule::get_qmLSSFreq ( )

inline

Definition at line 848 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSFreq; } ;

int* Molecule::get_qmLSSIdxs ( )

inline

Definition at line 846 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSIdxs;} ;

Mass* Molecule::get_qmLSSMass ( )

inline

Definition at line 847 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSMass; } ;

int* Molecule::get_qmLSSRefIDs ( )

inline

Definition at line 850 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSRefIDs ; } ;

Mass* Molecule::get_qmLSSRefMass ( )

inline

Definition at line 852 of file Molecule.h.

{return qmLSSRefMass; } ;

int* Molecule::get_qmLSSRefSize ( )

inline

Definition at line 851 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSRefSize ; } ;

int Molecule::get_qmLSSResSize ( )

inline

Definition at line 849 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSResidueSize; } ;

int* Molecule::get_qmLSSSize ( )

inline

Definition at line 845 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSSize;} ;

int* Molecule::get_qmMeMMindx ( )

inline

Definition at line 859 of file Molecule.h.

Referenced by ComputeQM::initialize().

{ return qmMeMMindx; } ;

int Molecule::get_qmMeNumBonds ( )

inline

Definition at line 858 of file Molecule.h.

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

{ return qmMeNumBonds; } ;

Real* Molecule::get_qmMeQMGrp ( )

inline

Definition at line 860 of file Molecule.h.

Referenced by ComputeQM::initialize().

{ return qmMeQMGrp; } ;

const int* const* Molecule::get_qmMMBond ( )

inline

Definition at line 837 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmMMBond; } ;

const int* const* Molecule::get_qmMMBondedIndx ( )

inline

Definition at line 839 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmMMBondedIndx; } ;

const int* const* Molecule::get_qmMMChargeTarget ( )

inline

Definition at line 842 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmMMChargeTarget; } ;

const int* Molecule::get_qmMMNumTargs ( )

inline

Definition at line 843 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmMMNumTargs; } ;

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

inline

Definition at line 853 of file Molecule.h.

{ return qmClassicSolv;} ;

int Molecule::get_qmNumBonds ( )

inline

Definition at line 834 of file Molecule.h.

{ return qmNumBonds; } ;

int Molecule::get_qmNumGrps ( ) const

inline

Definition at line 828 of file Molecule.h.

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

{return qmNumGrps; } ;

int Molecule::get_qmPCFreq ( )

inline

Definition at line 862 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmPCFreq; } ;

Bool Molecule::get_qmReplaceAll ( )

inline

Definition at line 855 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return qmReplaceAll; } ;

int Molecule::get_qmTotCustPCs ( )

inline

Definition at line 864 of file Molecule.h.

Referenced by ComputeQM::initialize().

{ return qmTotCustPCs; } ;

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

Definition at line 144 of file Molecule.C.

References NAMD_die().

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

                                            {

  if (atomNames == NULL || resLookup == NULL)
  {
    NAMD_die("Tried to find atom from name on node other than node 0");
  }
  int i = 0;
  int end = 0;
  if ( resLookup->lookup(segid,resid,&i,&end) ) return 0;
  return ( end - i );
}

void Molecule::get_rotdrag_params ( BigReal &  v, Vector &  a, Vector &  p, int  atomnum  ) const

inline

Definition at line 1331 of file Molecule.h.

Referenced by Sequencer::addRotDragToPosition().

  {
    v = rotDragParams[rotDragIndexes[atomnum]].v;
    a = rotDragParams[rotDragIndexes[atomnum]].a;
    p = rotDragParams[rotDragIndexes[atomnum]].p;
  }

unsigned char Molecule::get_ss_type ( int  anum ) const

inline

Definition at line 1355 of file Molecule.h.

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

        {
                return(ssAtomFlags[anum]);
        }

void Molecule::get_stir_refPos ( Vector &  refPos, int  atomnum  ) const

inline

Definition at line 1306 of file Molecule.h.

  {
    refPos = stirParams[stirIndexes[atomnum]].refPos;
  }

Real Molecule::get_stir_startTheta ( int  atomnum ) const

inline

Definition at line 1318 of file Molecule.h.

Referenced by ComputeStir::doForce().

  {
    return stirParams[stirIndexes[atomnum]].startTheta;
  }

Bond* Molecule::getAllAcceptors ( ) const

inline

Definition at line 1112 of file Molecule.h.

{return acceptors;}

Angle* Molecule::getAllAngles ( ) const

inline

Definition at line 1095 of file Molecule.h.

Referenced by buildAngleData().

{return angles;}

Bond* Molecule::getAllBonds ( ) const

inline

Definition at line 1094 of file Molecule.h.

Referenced by buildBondData().

{return bonds;}

Crossterm* Molecule::getAllCrossterms ( ) const

inline

Definition at line 1098 of file Molecule.h.

Referenced by buildCrosstermData().

{return crossterms;}

Dihedral* Molecule::getAllDihedrals ( ) const

inline

Definition at line 1097 of file Molecule.h.

Referenced by buildDihedralData().

{return dihedrals;}

Bond* Molecule::getAllDonors ( ) const

inline

Definition at line 1111 of file Molecule.h.

{return donors;}

Improper* Molecule::getAllImpropers ( ) const

inline

Definition at line 1096 of file Molecule.h.

Referenced by buildImproperData().

{return impropers;}

Lphost* Molecule::getAllLphosts ( ) const

inline

Definition at line 1102 of file Molecule.h.

{ return lphosts; }

BigReal Molecule::GetAtomAlpha ( int  i ) const

inline

Definition at line 483 of file Molecule.h.

References drude_constants::alpha.

                                    {
    return drudeConsts[i].alpha;
  }

AtomNameInfo* Molecule::getAtomNames ( ) const

inline

Definition at line 492 of file Molecule.h.

Referenced by buildAtomData().

{ return atomNames; }

Atom* Molecule::getAtoms ( ) const

inline

Definition at line 491 of file Molecule.h.

References atoms.

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

{ return atoms; }

AtomSegResInfo* Molecule::getAtomSegResInfo ( ) const

inline

Definition at line 495 of file Molecule.h.

Referenced by buildAtomData().

{ return atomSegResids; }

const float* Molecule::getBFactorData ( )

inline

Definition at line 1036 of file Molecule.h.

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

{ return (const float *)bfactor; }

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

Referenced by Controller::printEnergies().

int const* Molecule::getLcpoParamType ( )

inline

Definition at line 479 of file Molecule.h.

Referenced by HomePatch::setLcpoType().

                                 {
    return lcpoParamType;
  }

const float* Molecule::getOccupancyData ( )

inline

Definition at line 1032 of file Molecule.h.

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

{ return (const float *)occupancy; }

BigReal Molecule::getVirialTailCorr

(

const BigReal

)

Referenced by Controller::calcPressure().

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

{
  int i;      //  Loop counter
  for(i=0; i<MAX_RESTRICTIONS; i++) {
    if (go_array[(MAX_GO_CHAINS*chain1) + chain2].restrictions[i]  == rDiff) {
      return TRUE;
    } else if (go_array[(MAX_GO_CHAINS*chain1) + chain2].restrictions[i] == -1) {
      return FALSE;
    }
  }
  return FALSE;
}

void Molecule::goInit

(

)

Definition at line 54 of file GoMolecule.C.

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

                      {
  numGoAtoms=0;
  energyNative=0;
  energyNonnative=0;
  atomChainTypes=NULL;
  goSigmaIndices=NULL;
  goSigmas=NULL;
  goWithinCutoff=NULL;
  goCoordinates=NULL;
  goResids=NULL;
  goPDB=NULL;
}

void Molecule::initialize

(

)

Bool Molecule::is_atom_constorqued ( int  atomnum ) const

inline

Definition at line 1250 of file Molecule.h.

References consTorqueIndexes, FALSE, and numConsTorque.

  {
    if (numConsTorque)
    {
      //  Check the index to see if it is constrained
      return(consTorqueIndexes[atomnum] != -1);
    }
    else
    {
      //  No constraints at all, so just return FALSE
      return(FALSE);
    }
  }

Bool Molecule::is_atom_constrained ( int  atomnum ) const

inline

Definition at line 1201 of file Molecule.h.

References FALSE, and numConstraints.

Referenced by ComputeRestraints::doForce().

  {
    if (numConstraints)
    {
      //  Check the index to see if it is constrained
      return(consIndexes[atomnum] != -1);
    }
    else
    {
      //  No constraints at all, so just return FALSE
      return(FALSE);
    }
  }

Bool Molecule::is_atom_exPressure ( int  atomnum ) const

inline

Definition at line 1451 of file Molecule.h.

References numExPressureAtoms.

Referenced by Sequencer::langevinPiston().

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

Bool Molecule::is_atom_fixed ( int  atomnum ) const

inline

Definition at line 1411 of file Molecule.h.

References numFixedAtoms.

Referenced by WorkDistrib::createAtomLists().

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

Bool Molecule::is_atom_gridforced ( int  atomnum, int  gridnum  ) const

inline

Definition at line 1185 of file Molecule.h.

References FALSE, and numGridforceGrids.

Referenced by ComputeGridForce::do_calc().

  {
      if (numGridforceGrids)
      {
          return(gridfrcIndexes[gridnum][atomnum] != -1);
      }
      else
      {
          return(FALSE);
      }
  }

Bool Molecule::is_atom_movdragged ( int  atomnum ) const

inline

Definition at line 1218 of file Molecule.h.

References FALSE, and numMovDrag.

Referenced by Sequencer::addMovDragToPosition().

  {
    if (numMovDrag)
    {
      //  Check the index to see if it is constrained
      return(movDragIndexes[atomnum] != -1);
    }
    else
    {
      //  No constraints at all, so just return FALSE
      return(FALSE);
    }
  }

Bool Molecule::is_atom_rotdragged ( int  atomnum ) const

inline

Definition at line 1234 of file Molecule.h.

References FALSE, and numRotDrag.

Referenced by Sequencer::addRotDragToPosition().

  {
    if (numRotDrag)
    {
      //  Check the index to see if it is constrained
      return(rotDragIndexes[atomnum] != -1);
    }
    else
    {
      //  No constraints at all, so just return FALSE
      return(FALSE);
    }
  }

Bool Molecule::is_atom_stirred ( int  atomnum ) const

inline

Definition at line 1432 of file Molecule.h.

References FALSE, and numStirredAtoms.

Referenced by ComputeStir::doForce().

  {
    if (numStirredAtoms)
    {
      //  Check the index to see if it is constrained
      return(stirIndexes[atomnum] != -1);
    }
    else
    {
      //  No constraints at all, so just return FALSE
      return(FALSE);
    }
  }

Bool Molecule::is_drude

(

int

)

Bool Molecule::is_group_fixed ( int  atomnum ) const

inline

Definition at line 1447 of file Molecule.h.

References numFixedAtoms.

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

Bool Molecule::is_hydrogen

(

int

)

Bool Molecule::is_hydrogenGroupParent

(

int

)

Bool Molecule::is_lp

(

int

)

Bool Molecule::is_oxygen

(

int

)

Bool Molecule::is_water

(

int

)

Referenced by outputCompressedFile(), and wrap_coor_int().

Real Molecule::langevin_param ( int  atomnum ) const

inline

Definition at line 1300 of file Molecule.h.

Referenced by WorkDistrib::createAtomLists().

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

int64_t Molecule::num_deg_freedom ( int  isInitialReport = 0 ) const

inline

Definition at line 524 of file Molecule.h.

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

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

                                                         {
    // local variables prefixed by s_
    int64_t s_NumDegFreedom = 3 * (int64_t) numAtoms;
    int64_t s_NumFixedAtoms = num_fixed_atoms();
    if (s_NumFixedAtoms) s_NumDegFreedom -= 3 * s_NumFixedAtoms;
    if (numLonepairs) s_NumDegFreedom -= 3 * numLonepairs;
    if ( ! (s_NumFixedAtoms || numConstraints
          || simParams->comMove || simParams->langevinOn) ) {
      s_NumDegFreedom -= 3;
    }
    if ( ! isInitialReport && simParams->pairInteractionOn) {
      //
      // DJH: a kludge?  We want to initially report system degrees of freedom
      //
      // this doesn't attempt to deal with fixed atoms or constraints
      s_NumDegFreedom = 3 * numFepInitial;
    }
    int s_NumFixedRigidBonds = 
      (simParams->fixedAtomsOn ? numFixedRigidBonds : 0);
    if (simParams->watmodel == WAT_TIP4) {
      // numLonepairs is subtracted here because all lonepairs have a rigid bond
      // to oxygen, but all of the LP degrees of freedom are dealt with above
      s_NumDegFreedom -= (numRigidBonds - s_NumFixedRigidBonds - numLonepairs);
    }
    else {
      // Non-polarized systems don't have LPs.
      // For Drude model, bonds that attach LPs are not counted as rigid;
      // LPs have already been subtracted from degrees of freedom.
      s_NumDegFreedom -= (numRigidBonds - s_NumFixedRigidBonds);
    }
    return s_NumDegFreedom;
  }

int Molecule::num_fixed_atoms ( ) const

inline

Definition at line 498 of file Molecule.h.

References numFixedAtoms, and simParams.

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

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

int Molecule::num_fixed_groups ( ) const

inline

Definition at line 504 of file Molecule.h.

References num_fixed_atoms(), and numFixedGroups.

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

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

int64_t Molecule::num_group_deg_freedom ( ) const

inline

Definition at line 511 of file Molecule.h.

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

Referenced by Controller::receivePressure().

                                        {
    // local variables prefixed by s_
    int64_t s_NumGroupDegFreedom = 3 * (int64_t) numHydrogenGroups;
    int64_t s_NumFixedAtoms = num_fixed_atoms();
    int s_NumFixedGroups = num_fixed_groups();
    if (s_NumFixedGroups) s_NumGroupDegFreedom -= 3 * s_NumFixedGroups;
    if ( ! (s_NumFixedAtoms || numConstraints
          || simParams->comMove || simParams->langevinOn) ) {
      s_NumGroupDegFreedom -= 3;
    }
    return s_NumGroupDegFreedom;
  }

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 Parameters::assign_vdw_index(), qmSolvData::atmIDs, PDB::atom(), bond::atom1, bond::atom2, ResizeArray< T >::begin(), bond::bond_type, atom_constants::charge, charge, StringList::data, DebugM, PDBAtom::element(), ResizeArray< T >::end(), endi(), SimParameters::extraBondsOn, SortedArray< Type >::find(), ConfigList::find(), SimParameters::fixedAtomsOn, get_residue_size(), ObjectArena< Type >::getNewArray(), iERROR(), iINFO(), SortedArray< Type >::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< T >::size(), split(), SimParameters::stepsPerCycle, and PDBAtom::temperaturefactor().

Referenced by NamdState::loadStructure().

                                                                        {
#ifdef MEM_OPT_VERSION
    NAMD_die("QMMM interface is not supported in memory optimized builds.");
#else
    
    PDB *pdbP;      //  Pointer to PDB object to use
    
    qmNumQMAtoms = 0;
    
    qmNumGrps = 0 ;
    
    iout << iINFO << "Using the following PDB file for QM parameters: " << 
    pdbFileName << "\n" << endi;
    if (pdbFileName)
        pdbP = new PDB(pdbFileName);
    else
        iout << "pdbFile name not defined!\n\n" << endi ;
    
    if (pdbP->num_atoms() != numAtoms)
    {
       NAMD_die("Number of atoms in QM parameter PDB file doesn't match coordinate PDB");
    }
    
    qmElementArray = new char*[numAtoms] ;
    
    qmAtomGroup = new Real[numAtoms] ;
    
    BigReal bondVal;
    Real atmGrp;
    
    std::set<Real> atmGrpSet;
    
    std::vector<Real> grpChrgVec;
    
    // Value in the bond column fro QM-MM bonds.
    std::vector<BigReal> qmBondValVec;
    // Identifiers of different QM regions
    std::vector<Real> qmGroupIDsVec;
    // This maps the qm group ID with the serail index of this group in 
    // various arrays.
    std::map<Real,int> qmGrpIDMap ;
    
    std::vector<int> qmAtmIndxVec, qmGrpSizeVec;
    
    // Used to parse user selection in different places.
    std::vector<std::string> strVec;
    
    qmNoPC = simParams->qmNoPC;
    
    // We set to zero by default, So there is no need for extra processing.
    qmPCFreq = 0;
    if (simParams->qmPCSelFreq > 1)
        qmPCFreq = simParams->qmPCSelFreq;
    
    
    
    
    qmLSSTotalNumAtms = 0;
    SortedArray< qmSolvData> lssGrpRes;
    std::map<Real,std::vector<int> > lssGrpRefIDs;
    refSelStrMap lssRefUsrSel;
    int totalNumRefAtms = 0;
    int lssClassicResIndx = 0 ;
    int lssCurrClassResID = -1 ;
    char lssCurrClassResSegID[5];
    if (simParams->qmLSSOn) {
        DebugM(4, "LSS is ON! Processing QM solvent.\n") ;
        
        if (resLookup == NULL)
            NAMD_die("We need residue data to conduct LSS.");
         
        if (simParams->qmLSSMode == QMLSSMODECOM) {
            
            StringList *current = cfgList->find("QMLSSRef");
            for ( ; current; current = current->next ) {
                
                strVec = split( std::string(current->data) , " ");
                
                if (strVec.size() != 3 ) {
                    iout << iERROR << "Format error in QM LSS size: " 
                    << current->data
                    << "\n" << endi;
                    NAMD_die("Error processing QM information.");
                }
                
                std::stringstream storConv ;
                
                storConv << strVec[0] ;
                Real grpID ;
                storConv >> grpID;
                if (storConv.fail()) {
                    iout << iERROR << "Error parsing QMLSSRef selection: " 
                    << current->data
                    << "\n" << endi;
                    NAMD_die("Error processing QM information.");
                }
                
                std::string segName = strVec[1].substr(0,4);
                
                storConv.clear() ;
                storConv << strVec[2];
                int resID ;
                storConv >> resID;
                if (storConv.fail()) {
                    iout << iERROR << "Error parsing QMLSSRef selection: " 
                    << current->data
                    << "\n" << endi;
                    NAMD_die("Error processing QM information.");
                }
                
                auto it = lssRefUsrSel.find(grpID) ;
                if (it == lssRefUsrSel.end())
                    lssRefUsrSel.insert(refSelStrPair(grpID,refSelStrVec()));
                
                lssRefUsrSel[grpID].push_back(refSelStr(segName, resID));
            }
            
            for (auto it=lssRefUsrSel.begin(); it!=lssRefUsrSel.end(); it++) {
                iout << iINFO << "LSS user selection for COM composition of group "
                << it->first << ":\n" << endi ;
                for (int i=0; i<it->second.size();i++) {
                    iout << iINFO << "Segment " << it->second[i].segid 
                    << " ; residue " << it->second[i].resid << "\n" << endi ;
                }
            }
        }
    }
    
    
    
    
    
    for (int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
        
        // If we are looking for QM-MM bonds, then we need to store extra info.
        if (simParams->qmBondOn) {
            
            // We store both the qm group and the bond value
            if ( strcmp(simParams->qmColumn,"beta") == 0 ){
                atmGrp = pdbP->atom(atmInd)->temperaturefactor() ;
                
                bondVal = pdbP->atom(atmInd)->occupancy() ;
            }
            else {
                atmGrp = pdbP->atom(atmInd)->occupancy() ;
                
                bondVal = pdbP->atom(atmInd)->temperaturefactor() ;
            }
            
            qmBondValVec.push_back(bondVal);
        }
        else {
            
            if ( strcmp(simParams->qmColumn,"beta") == 0 ){
                atmGrp = pdbP->atom(atmInd)->temperaturefactor() ;
            }
            else {
                atmGrp = pdbP->atom(atmInd)->occupancy() ;
            }
        }
        
        // We store all atom QM Group IDs in an array for 
        // future transport to all PEs.
        qmAtomGroup[atmInd] = atmGrp;
        
        // We store all atom's elements for quick access in the QM code.
        // Elements are used to tell the QM code the atom type.
        qmElementArray[atmInd] = new char[3];
        strncpy(qmElementArray[atmInd],pdbP->atom(atmInd)->element(),3);
        
        // For QM atoms, we keep global atom indices and parse different QM Group
        // IDs, keeping track of each groups size
        if (atmGrp > 0){
            
            if (simParams->fixedAtomsOn) {
                if ( fixedAtomFlags[atmInd] == 1 ) {
                    iout << iERROR << "QM atom cannot be fixed in space!\n" 
                    << endi;
                    NAMD_die("Error processing QM information.");
                }
            }
            
            // Changes the VdW type of QM atoms.
            // This is may be used to counter balance the overpolarization 
            // that QM atoms sufer.
            if (simParams->qmVDW) {
                // Creating a new type string
                std::string qmType(qmElementArray[atmInd]) ;
                // Erases empty spaces
                qmType.erase(
                        std::remove_if(qmType.begin(), qmType.end(), isspace ),
                    qmType.end());
                // pre-pends a "q" to the element name.
                qmType = std::string("q") + qmType;
            
//                 iout << "QM VdW type: " << atoms[atmInd].vdw_type 
//                 << " atom type: " << atomNames[atmInd].atomtype 
//                 << " new type "<< qmType << "\n" << endi;
                
                /*  Look up the vdw constants for this atom    */
                // I am casting a non-const here because the function doesn't actually
                // change the string value, but it doesn't ask for a const char* as 
                // an argument.
                params->assign_vdw_index(const_cast<char*>( qmType.c_str()), 
                                         &(atoms[atmInd]));
                
//                 iout << "----> new VdW type: " << atoms[atmInd].vdw_type << "\n" << endi;
            }
            
            // Indexes all global indices of QM atoms.
            qmAtmIndxVec.push_back(atmInd);
            
            auto retVal = atmGrpSet.insert(atmGrp);
            
            // If a new QM group is found, store the reverse reference in a map
            // and increase the total count.
            if (retVal.second) {
                // This mak makes the reverse identification from the group ID
                // to the sequential order in which each group was first found.
                qmGrpIDMap.insert(std::pair<BigReal,int>(atmGrp, qmNumGrps));
                
                // This vector keeps track of the group ID for each group
                qmGroupIDsVec.push_back(atmGrp);
                
                // This counter is used to keep track of the sequential order in
                // which QM groups were first seen in the reference PDB file.
                qmNumGrps++ ;
                
                // If this is a new QM group, initialize a new variable in the 
                // vector to keep track of the number of atoms in this group.
                qmGrpSizeVec.push_back(1);
                
                // For each new QM group, a new entry in the total charge
                // vector is created
                grpChrgVec.push_back( atoms[atmInd].charge );
            }
            else {
                // If the QM group has already been seen, increase the count
                // of atoms in that group.
                qmGrpSizeVec[ qmGrpIDMap[atmGrp] ]++ ;
                
                // If the QM group exists, adds current atom charge to total charge.
                grpChrgVec[ qmGrpIDMap[atmGrp] ] += atoms[atmInd].charge;
            }
            
            // In case we are using live solvent selection
            if(simParams->qmLSSOn) {
                
                int resID = pdbP->atom(atmInd)->residueseq();
                char segName[5];
                strncpy(segName, pdbP->atom(atmInd)->segmentname(),5);
                
                int resSize = get_residue_size(segName,resID);
                
                int i =-1, end =-1;
                
                resLookup->lookup(segName,resID,&i, &end);
                
                if ( strncasecmp(pdbP->atom(atmInd)->residuename(),
                           simParams->qmLSSResname, 4) == 0) {
                    
                    // We try to insert every residue from every atom
                    qmSolvData *resP = lssGrpRes.find(qmSolvData(resID, i, resSize));
                    
                    if (resP != NULL) {
                        resP->atmIDs.push_back(atmInd) ;
//                         DebugM(3, "Existing residue " << resP->resID 
//                         << " from segID " << resP->segName
//                         << " received atom "
//                         << atmInd << "\n" );
                    }
                    else {
                        lssGrpRes.insert(qmSolvData(resID, i, resSize, segName, atmGrp));
//                         DebugM(3, lssGrpRes.size() << ") Inserted new residue: " 
//                         << resID << " from segID " << segName
//                         << " with atom "
//                         << i << "\n" ) ;
                    }
                }
                else {
                    // We store the QM atoms, per group, which are NOT part of
                    // solvent molecules.
                    
                    // Checks if we have a vector for this QM group.
                    auto itNewGrp = lssGrpRefIDs.find(atmGrp) ;
                    if (itNewGrp == lssGrpRefIDs.end()) {
                        lssGrpRefIDs.insert(std::pair<Real, std::vector<int> >(
                            atmGrp, std::vector<int>() ) );
                    }
                    
                    switch (simParams->qmLSSMode)
                    {
                    
                    case QMLSSMODECOM:
                    {
                        // If we are using COM selection, checks if the atom
                        // is part of the user selected residues
                        for(int i=0; i<lssRefUsrSel[atmGrp].size(); i++) {
                            if (lssRefUsrSel[atmGrp][i].resid == resID &&
                                strncmp(lssRefUsrSel[atmGrp][i].segid.c_str(),
                                        segName,5) == 0 ) {
                                
                                lssGrpRefIDs[atmGrp].push_back(atmInd);
                                totalNumRefAtms++;
                            }
                        }
                        
                    } break;
                    
                    case QMLSSMODEDIST:
                    {
                    
                        lssGrpRefIDs[atmGrp].push_back(atmInd);
                        totalNumRefAtms++;
                    
                    } break;
                    
                    }
                }
                
            }
            
        }
        else if (atmGrp == 0) {
            
            if(simParams->qmLSSOn) {
                
                if ( strncasecmp(pdbP->atom(atmInd)->residuename(),
                           simParams->qmLSSResname, 4) == 0) {
                    
                    int resID = pdbP->atom(atmInd)->residueseq();
                    char segName[5];
                    strncpy(segName, pdbP->atom(atmInd)->segmentname(),5);
                    
                    if (lssCurrClassResID < 0) {
                        lssCurrClassResID = resID ;
                        strncpy(lssCurrClassResSegID, segName,5);
                        lssClassicResIndx = 0;
                    }
                    else if (lssCurrClassResID != resID ||
                            strcmp(lssCurrClassResSegID, segName) != 0 ) {
                        lssCurrClassResID = resID ;
                        strncpy(lssCurrClassResSegID, segName,5);
                        lssClassicResIndx++;
                    }
                    
                    qmClassicSolv.insert(std::pair<int,int>(atmInd,lssClassicResIndx));
                    
                }
            }
        }
        else if(atmGrp < 0) {
            iout << iERROR << "QM group ID cannot be less than zero!\n" << endi;
            NAMD_die("Error processing QM information.");
        }
    }
    
    // Sanity check
    if (simParams->qmLSSOn) {
        if (lssGrpRes.size() == 0)
            NAMD_die("LSS was activated but there are no QM solvent molecules!\n");
        
        for (auto it=lssRefUsrSel.begin(); it!=lssRefUsrSel.end(); it++) {
            auto itTarget = qmGrpIDMap.find(it->first);
            if (itTarget == qmGrpIDMap.end()) {
                iout << iERROR << "QM group ID for LSS could not be found in input!"
                << " QM group ID: " << it->first << "\n" << endi;
                NAMD_die("Error processing QM information.");
            }
        }
        
        DebugM(3,"We found " << lssClassicResIndx << " classical solvent molecules totalling "
        << qmClassicSolv.size() << " atoms.\n" );
        
    }
    
    qmNumQMAtoms = qmAtmIndxVec.size();
    
    if (qmNumQMAtoms == 0)
        NAMD_die("No QM atoms were found in this QM simulation!") ;
    
    iout << iINFO << "Number of QM atoms (excluding Dummy atoms): " << 
    qmNumQMAtoms << "\n" << endi;
    
    qmAtmChrg = new Real[qmNumQMAtoms] ;
    qmAtmIndx = new int[qmNumQMAtoms] ;
    for (int i=0; i<qmNumQMAtoms; i++) {
        // qmAtmChrg gets initialized with the PSF charges at the end of this
        // function, but values may change as QM steps are taken.
        qmAtmChrg[i] = 0;  
        qmAtmIndx[i] = qmAtmIndxVec[i] ;
    }
    
    // This map relates the QM group index with a vector of pairs
    // of bonded MM-QM atoms (with the bonded atoms ins this order: 
    // MM first, QM second).
    std::map<int, std::vector<std::pair<int,int> > > qmGrpIDBonds ;
    int bondCounter = 0;
    if (simParams->qmBondOn) {
        
        // Checks all bonds for QM-MM pairs.
        // Makes sanity checks against QM-QM pairs and MM-MM pairs that
        // are flagged by the user to be bonds between QM and MM regions.
        // QM-QM bonds will be removed in another function.
        for (int bndIt = 0; bndIt < numBonds; bndIt++) {
            
            bond curr = bonds[bndIt] ;
            
            // In case either atom has a non-zero
            if ( qmBondValVec[curr.atom1] != 0 &&
                 qmBondValVec[curr.atom2] != 0 ) {
                
                // Sanity checks (1 of 2)
                if (qmAtomGroup[curr.atom1] != 0 &&
                    qmAtomGroup[curr.atom2] != 0) {
                    iout << iERROR << "Atoms " << curr.atom1 << " and " << 
                    curr.atom2 << " are assigned as QM atoms.\n" << endi;
                    NAMD_die("Error in QM-MM bond assignment.") ;
                }
                
                // Sanity checks (2 of 2)
                if (qmAtomGroup[curr.atom1] == 0 &&
                    qmAtomGroup[curr.atom2] == 0) {
                    iout << iERROR << "Atoms " << curr.atom1 << " and " << 
                    curr.atom2 << " are assigned as MM atoms.\n" << endi;
                    NAMD_die("Error in QM-MM bond assignment.") ;
                }
                
                int currGrpID = 0;
                std::pair<int,int> newPair(0,0);
                
                // We create a QM-MM pair with the MM atom first
                if (qmAtomGroup[curr.atom1] != 0) {
                    newPair = std::pair<int,int>(curr.atom2, curr.atom1) ;
                    currGrpID = qmAtomGroup[curr.atom1];
                } else {
                    newPair = std::pair<int,int>(curr.atom1, curr.atom2) ;
                    currGrpID = qmAtomGroup[curr.atom2];
                } 
                
                int grpIndx = qmGrpIDMap[currGrpID] ;
                
                // Stores bonds in vectors key-ed by the QM group ID of the QM atom.
                auto retIter = qmGrpIDBonds.find(grpIndx) ;
                // In case thi QM-MM bonds belong to a QM group we have not seen 
                // yet, stores a new vector in the map.
                if (retIter == qmGrpIDBonds.end()) {
                    qmGrpIDBonds.insert(std::pair<BigReal,std::vector<std::pair<int,int> > >(
                    grpIndx, std::vector<std::pair<int,int> >() ) );
                }
                
                qmGrpIDBonds[grpIndx].push_back( newPair );
                
                bondCounter++ ;
            }
            
            
        }
        
//         iout << "Finished processing QM-MM bonds.\n" << endi ;
        
        if (bondCounter == 0)
            iout << iWARN << "We found " << bondCounter << " QM-MM bonds.\n" << endi ;
        else
            iout << iINFO << "We found " << bondCounter << " QM-MM bonds.\n" << endi ;
        
    }
    
    // Initializes several arrays used to setup the QM simulation.
    qmNumBonds = bondCounter ;
    
    qmMMBond = new int*[qmNumBonds];
    
    qmDummyBondVal = new BigReal[qmNumBonds];
    
    qmMMChargeTarget = new int*[qmNumBonds] ;
    qmMMNumTargs = new int[qmNumBonds] ;
    
    qmDummyElement = new char*[qmNumBonds] ;
    
    
    qmNumGrps = atmGrpSet.size();
    
    qmGrpSizes = new int[qmNumGrps];
    
    qmGrpID = new Real[qmNumGrps];
    
    qmGrpChrg = new Real[qmNumGrps];
    
    qmGrpMult = new Real[qmNumGrps] ;
    
    qmGrpNumBonds = new int[qmNumGrps];
    
    qmGrpBonds = new int*[qmNumGrps];
    qmMMBondedIndx = new int*[qmNumGrps];
    
    
    
    
    // We first initialize the multiplicity vector with 
    // default values, then populate it with user defined values.
    for (int grpIter = 0; grpIter < qmNumGrps; grpIter++) {
        qmGrpMult[grpIter] = 0;
    }
    
    int multCount = 0;
    StringList *current = cfgList->find("QMMult");
    for ( ; current; current = current->next ) {
        
        auto strVec = split( std::string(current->data) , " ");
        
        if (strVec.size() != 2 ) {
            iout << iERROR << "Format error in QM Multiplicity string: " 
            << current->data
            << "\n" << endi;
            NAMD_die("Error processing QM information.");
        }
        
        std::stringstream storConv ;
        
        storConv << strVec[0] ;
        Real grpID ;
        storConv >> grpID;
        if (storConv.fail()) {
            iout << iERROR << "Error parsing QMMult selection: " 
            << current->data
            << "\n" << endi;
            NAMD_die("Error processing QM information.");
        }
        
        storConv.clear() ;
        storConv << strVec[1];
        Real multiplicity ;
        storConv >> multiplicity;
        if (storConv.fail()) {
            iout << iERROR << "Error parsing QMMult selection: " 
            << current->data
            << "\n" << endi;
            NAMD_die("Error processing QM information.");
        }
        
        auto it = qmGrpIDMap.find(grpID);
        
        if (it == qmGrpIDMap.end()) {
            iout << iERROR << "Error parsing QMMult selection. Could not find QM group ID: " 
            << grpID
            << "\n" << endi;
            NAMD_die("Error processing QM information.");
        }
        else {
            iout << iINFO << "Applying user defined multiplicity "
            << multiplicity << " to QM group ID " << grpID << "\n" << endi;
            qmGrpMult[it->second] = multiplicity;
        }
        
        multCount++;
    }
    
    if (multCount != qmNumGrps && simParams->qmFormat == QMFormatORCA) {
        NAMD_die("ORCA neds multiplicity values for all QM regions.");
    }
    
    
    
    
    for (size_t grpIndx = 0; grpIndx < qmGrpSizeVec.size(); grpIndx++) {
        
        bool nonInteger = true;
        if ((fabsf(roundf(grpChrgVec[grpIndx]) - grpChrgVec[grpIndx]) <= 0.001f) ) {
            grpChrgVec[grpIndx] = roundf(grpChrgVec[grpIndx]) ;
            nonInteger = false;
        }
        
        iout << iINFO << grpIndx + 1 << ") Group ID: " << qmGroupIDsVec[grpIndx]
        << " ; Group size: " << qmGrpSizeVec[grpIndx] << " atoms"
        << " ; Total PSF charge: " << grpChrgVec[grpIndx] << "\n" << endi ;
        
        if (nonInteger && simParams->PMEOn)
            NAMD_die("QM atoms do not add up to a whole charge, which is needed for PME.") ;
    }
    
    int chrgCount = 0;
    current = cfgList->find("QMCharge");
    for ( ; current; current = current->next ) {
        
        auto strVec = split( std::string(current->data) , " ");
        
        if (strVec.size() != 2 ) {
            iout << iERROR << "Format error in QM Charge string: " 
            << current->data
            << "\n" << endi;
            NAMD_die("Error processing QM information.");
        }
        
        std::stringstream storConv ;
        
        storConv << strVec[0] ;
        Real grpID ;
        storConv >> grpID;
        if (storConv.fail()) {
            iout << iERROR << "Error parsing QMCharge selection: " 
            << current->data
            << "\n" << endi;
            NAMD_die("Error processing QM information.");
        }
        
        storConv.clear() ;
        storConv << strVec[1];
        Real charge ;
        storConv >> charge;
        if (storConv.fail()) {
            iout << iERROR << "Error parsing QMCharge selection: " 
            << current->data
            << "\n" << endi;
            NAMD_die("Error processing QM information.");
        }
        
        auto it = qmGrpIDMap.find(grpID);
        
        if (it == qmGrpIDMap.end()) {
            iout << iERROR << "Error parsing QMCharge selection. Could not find QM group ID: " 
            << grpID
            << "\n" << endi;
            NAMD_die("Error processing QM information.");
        }
        else {
            iout << iINFO << "Found user defined charge "
            << charge << " for QM group ID " << grpID << ". Will ignore PSF charge.\n" << endi;
            grpChrgVec[it->second] = charge;
        }
        
        chrgCount++;
    }
    
    simParams->qmMOPACAddConfigChrg = false;
    // Checks if QM group charges can be defined for MOPAC.
    // Since charge can be defined in QMConfigLine, we need extra logic.
    if (simParams->qmFormat == QMFormatMOPAC) {
        
        // Checks if group charge was supplied in config line for MOPAC.
        std::string::size_type chrgLoc = std::string::npos ;
        
        // This will hold a sting with the first user supplied configuration line for MOPAC.
        std::string configLine(cfgList->find("QMConfigLine")->data) ;
        chrgLoc = configLine.find("CHARGE") ;
        
        if ( chrgLoc != std::string::npos ) {
            iout << iINFO << "Found user defined charge in command line. This \
will be used for all QM regions and will take precedence over all other charge \
definitions.\n" << endi;
        }
        else if ( chrgLoc == std::string::npos && (simParams->qmChrgFromPSF || chrgCount == qmNumGrps) ) {
        // If no charge was defined in the configuration line, gets from PSF and/or 
        // from user defined charges (through the QMCharge keyword).
            simParams->qmMOPACAddConfigChrg = true;
        }
        else
        {
        // If we could nont find a charge definition in the config line AND
        // no specific charge was selected for each QM region through QMCharge AND
        // the QMChargeFromPSF was not turned ON, the we stop NAMD and scream at the user.
//         if ( chrgLoc == std::string::npos && (chrgCount != qmNumGrps ) && !simParams->qmChrgFromPSF) 
            iout << iERROR << "Could not find charge for all QM groups. For MOPAC, \
charges can be defined through QMConfigLine, QMCharge or QMChargeFromPSF keywords.\n" << endi;
            NAMD_die("Error processing QM information.");
        }
        
    }
    
    if (simParams->qmFormat == QMFormatORCA ) {
        if ((chrgCount != qmNumGrps ) && !simParams->qmChrgFromPSF) {
            // If we are not supposed to get charges from the PSF, and not 
            // enough charges were set to cover all QM regions, 
            // we stop NAMD and scream at the user.
            iout << iERROR << "Could not find charge for all QM groups. For ORCA, \
charges can be defined through QMCharge or QMChargeFromPSF keywords.\n" << endi;
            NAMD_die("Error processing QM information.");
        }
    }
    
    // If mechanichal embedding was requested but we have QM-MM bonds, we need 
    // to send extra info to ComputeQM to preserve calculation speed.
    if (qmNumBonds > 0 && qmNoPC) {
        qmMeNumBonds = qmNumBonds;
        qmMeMMindx = new int[qmMeNumBonds] ;
        qmMeQMGrp = new Real[qmMeNumBonds] ;
    }
    else {
        qmMeNumBonds = 0 ;
        qmMeMMindx = NULL;
        qmMeQMGrp = NULL;
    }
    
    
    
    
    bondCounter = 0;
    for (int grpIter = 0; grpIter < qmNumGrps; grpIter++) {
        
        qmGrpID[grpIter] = qmGroupIDsVec[grpIter] ;
        
        qmGrpSizes[grpIter] = qmGrpSizeVec[grpIter] ;
        
        qmGrpChrg[grpIter] = grpChrgVec[grpIter];
        
//         iout << "Loaded " << qmGrpSizes[grpIter] << " atoms to this group.\n" << endi;
        
        int currNumbBonds = qmGrpIDBonds[grpIter].size() ;
        
        // Assigns the number of bonds that the current QM group has.
        qmGrpNumBonds[grpIter] = currNumbBonds;
        
        if (currNumbBonds > 0) {
            
            qmGrpBonds[grpIter] = new int[currNumbBonds];
            qmMMBondedIndx[grpIter] = new int[currNumbBonds];
            
            for (int bndIter = 0; bndIter<currNumbBonds; bndIter++) {
                
                // Adds the bonds to the overall sequential list.
                qmMMBond[bondCounter] = new int[2] ;
                qmMMBond[bondCounter][0] = qmGrpIDBonds[grpIter][bndIter].first ;
                qmMMBond[bondCounter][1] = qmGrpIDBonds[grpIter][bndIter].second ;
                
                // For the current QM group, and the current bond, gets the bond index.
                qmGrpBonds[grpIter][bndIter] =  bondCounter;
                
                // For the current QM group, and the current bond, gets the MM atom.
                qmMMBondedIndx[grpIter][bndIter] = qmGrpIDBonds[grpIter][bndIter].first ;
                
                // Assign the default value of dummy element
                qmDummyElement[bondCounter] = new char[3];
                strcpy(qmDummyElement[bondCounter],"H\0");
                
                // Determines the distance that will separate the new Dummy atom
                // and the Qm atom to which it will be bound.
                bondVal = 0;
                if (simParams->qmBondDist) {
                    if (qmBondValVec[qmMMBond[bondCounter][0]] != 
                        qmBondValVec[qmMMBond[bondCounter][1]]
                    ) {
                        iout << iERROR << "qmBondDist is ON but the values in the bond column are different for atom "
                        << qmMMBond[bondCounter][0] << " and " << qmMMBond[bondCounter][1]
                        << ".\n" << endi ;
                        NAMD_die("Error in QM-MM bond processing.");
                    }
                    
                    bondVal = qmBondValVec[qmMMBond[bondCounter][0]] ;
                } 
                else {
                    
                    if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"C") == 0 ) {
                        bondVal = 1.09 ;
                    }
                    else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"O") == 0 ) {
                        bondVal = 0.98 ;
                    }
                    else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"N") == 0 ) {
                        bondVal = 0.99 ;
                    }
                    else {
                        iout << iERROR << "qmBondDist is OFF but the QM atom has no default distance value. Atom "
                        << qmMMBond[bondCounter][1] << ", with element: " 
                        << qmElementArray[qmMMBond[bondCounter][1]] 
                        <<  ".\n" << endi ;
                        NAMD_die("Error in QM-MM bond processing.");
                    }
                    
                }
                
                iout << iINFO << "MM-QM pair: " << qmMMBond[bondCounter][0] << ":"
                << qmMMBond[bondCounter][1] 
                << " -> Value (distance or ratio): " << bondVal
                << " (QM Group " << grpIter << " ID " << qmGrpID[grpIter] << ")"
                << "\n" << endi ;
                
                qmDummyBondVal[bondCounter] = bondVal;
                
                // In case we are preparing for a mechanical embedding simulation
                // with no point charges, populate the following vectors
                if (qmMeNumBonds > 0) {
                    qmMeMMindx[bondCounter] = qmMMBond[bondCounter][0];
                    qmMeQMGrp[bondCounter] = qmGrpID[grpIter];
                }
                
                bondCounter++ ;
            }
        }
        
    }
    
    
    
    current = NULL;
    if (qmNumBonds > 0)
        current = cfgList->find("QMLinkElement");
    
    int numParsedLinkElem = 0;
    for ( ; current != NULL; current = current->next ) {
        
        DebugM(3,"Parsing link atom element: " << current->data << "\n" );
        
        strVec = split( std::string(current->data) , " ");
        
        // We need the two atoms that compose the QM-MM bonds and 
        // then the element.
        if (strVec.size() != 3 && qmNumBonds > 1) {
            iout << iERROR << "Format error in QM link atom element selection: " 
            << current->data
            << "\n" << endi;
            NAMD_die("Error processing QM information.");
        }
        
        // If there is only one QM-MM bond, we can accept the element only.
        if (strVec.size() != 1 && strVec.size() != 3 && qmNumBonds == 1) {
            iout << iERROR << "Format error in QM link atom element selection: " 
            << current->data
            << "\n" << endi;
            NAMD_die("Error processing QM information.");
        }
        
        std::stringstream storConv ;
        int bondAtom1, bondAtom2;
        std::string linkElement ;
        
        if (strVec.size() == 1) {
            linkElement = strVec[0].substr(0,2);
        }
        else if (strVec.size() == 3) {
            
            storConv << strVec[0] ;
            storConv >> bondAtom1;
            if (storConv.fail()) {
                iout << iERROR << "Error parsing link atom element selection: " 
                << current->data
                << "\n" << endi;
                NAMD_die("Error processing QM information.");
            }
            
            storConv.clear() ;
            storConv << strVec[1];
            storConv >> bondAtom2;
            if (storConv.fail()) {
                iout << iERROR << "Error parsing link atom element selection: " 
                << current->data
                << "\n" << endi;
                NAMD_die("Error processing QM information.");
            }
            
            linkElement = strVec[2].substr(0,2);
        }
        
        numParsedLinkElem++;
        
        if (numParsedLinkElem>qmNumBonds) {
            iout << iERROR << "More link atom elements were set than there"
            " are QM-MM bonds.\n" << endi;
            NAMD_die("Error processing QM information.");
        }
        
        int bondIter;
        
        if (strVec.size() == 1) {
            bondIter = 0;
        }
        else if (strVec.size() == 3) {
            
            Bool foundBond = false;
            
            for (bondIter=0; bondIter<qmNumBonds; bondIter++) {
                
                if ( (qmMMBond[bondIter][0] == bondAtom1 &&
                    qmMMBond[bondIter][1] == bondAtom2 ) ||
                    (qmMMBond[bondIter][0] == bondAtom2 &&
                    qmMMBond[bondIter][1] == bondAtom1 ) ) {
                    
                    foundBond = true;
                    break;
                }
            }
            
            if (! foundBond) {
                iout << iERROR << "Error parsing link atom element selection: " 
                << current->data
                << "\n" << endi;
                iout << iERROR << "Atom pair was not found to be a QM-MM bond.\n"
                << endi;
                NAMD_die("Error processing QM information.");
            }
        }
        
        strcpy(qmDummyElement[bondIter],linkElement.c_str());
        qmDummyElement[bondIter][2] = '\0';
        
        iout << iINFO << "Applying user defined link atom element "
        << qmDummyElement[bondIter] << " to QM-MM bond "
        << bondIter << ": " << qmMMBond[bondIter][0]
        << " - " << qmMMBond[bondIter][1]
        << "\n" << endi;
    }
    
    
    
    
    
    int32 **bondsWithAtomLocal = NULL ;
    int32 *byAtomSizeLocal = NULL;
    ObjectArena <int32 >* tmpArenaLocal = NULL;
    if (qmNumBonds > 0) {
        
        bondsWithAtomLocal = new int32 *[numAtoms];
        byAtomSizeLocal = new int32[numAtoms];
        tmpArenaLocal = new ObjectArena<int32>;
        
        //  Build the bond lists
       for (int i=0; i<numAtoms; i++)
       {
         byAtomSizeLocal[i] = 0;
       }
       for (int i=0; i<numRealBonds; i++)
       {
         byAtomSizeLocal[bonds[i].atom1]++;
         byAtomSizeLocal[bonds[i].atom2]++;
       }
       for (int i=0; i<numAtoms; i++)
       {
         bondsWithAtomLocal[i] = tmpArenaLocal->getNewArray(byAtomSizeLocal[i]+1);
         bondsWithAtomLocal[i][byAtomSizeLocal[i]] = -1;
         byAtomSizeLocal[i] = 0;
       }
       for (int i=0; i<numRealBonds; i++)
       {
         int a1 = bonds[i].atom1;
         int a2 = bonds[i].atom2;
         bondsWithAtomLocal[a1][byAtomSizeLocal[a1]++] = i;
         bondsWithAtomLocal[a2][byAtomSizeLocal[a2]++] = i;
       }
    }
    
    // In this loops we try to find other bonds in which the MM atoms from
    // QM-MM bonds may be involved. The other MM atoms (which we will call M2 and M3)
    // will be involved in charge manipulation. See ComputeQM.C for details.
    for (int qmBndIt = 0; qmBndIt < qmNumBonds; qmBndIt++) {
        
        // The charge targets are accumulated in a temporary vector and then 
        // transfered to an array that will be transmited to the ComputeQMMgr object.
        std::vector<int> chrgTrgt ;
        
        int MM1 = qmMMBond[qmBndIt][0], MM2, MM2BondIndx, MM3, MM3BondIndx;
        
        switch (simParams->qmBondScheme) {
            
            case QMSCHEMERCD:
            
            case QMSCHEMECS:
            {
                // Selects ALL MM2 atoms.
                for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
                    MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
                    
                    // Checks which of the atoms in the bond structure is the
                    // MM2 atom.
                    if (bonds[MM2BondIndx].atom1 == MM1)
                        MM2 = bonds[MM2BondIndx].atom2;
                    else
                        MM2 = bonds[MM2BondIndx].atom1;
                    
                    // In case the bonded atom is a QM atom,
                    // skips the index.
                    if (qmAtomGroup[MM2] > 0)
                        continue;
                    
                    chrgTrgt.push_back(MM2);
                }
                
            } break;
            
            case QMSCHEMEZ3:
            {
                // Selects all MM3 atoms
                for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
                    MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
                    
                    // Checks which of the atoms in the bond structure is the
                    // MM2 atom.
                    if (bonds[MM2BondIndx].atom1 == MM1)
                        MM2 = bonds[MM2BondIndx].atom2;
                    else
                        MM2 = bonds[MM2BondIndx].atom1;
                    
                    // In case the bonded atom is a QM atom,
                    // skips the index.
                    if (qmAtomGroup[MM2] > 0)
                        continue;
                    
                    for (int i=0; i<byAtomSizeLocal[MM2]; i++) {
                        MM3BondIndx = bondsWithAtomLocal[MM2][i] ;
                        
                        // Checks which of the atoms in the bond structure is the
                        // MM3 atom.
                        if (bonds[MM3BondIndx].atom1 == MM2)
                            MM3 = bonds[MM3BondIndx].atom2;
                        else
                            MM3 = bonds[MM3BondIndx].atom1;
                        
                        // In case the bonded atom is a QM atom,
                        // skips the index.
                        // We also keep the search from going back to MM1.
                        if (qmAtomGroup[MM3] > 0 || MM3 == MM1)
                            continue;
                        
                        chrgTrgt.push_back(MM3);
                    }
                    
                }
                
            };
            
            case QMSCHEMEZ2:
            {
                // Selects all MM2 atoms
                for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
                    MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
                    
                    // Checks which of the atoms in the bond structure is the
                    // MM2 atom.
                    if (bonds[MM2BondIndx].atom1 == MM1)
                        MM2 = bonds[MM2BondIndx].atom2;
                    else
                        MM2 = bonds[MM2BondIndx].atom1;
                    
                    // In case the bonded atom is a QM atom,
                    // skips the index.
                    if (qmAtomGroup[MM2] > 0)
                        continue;
                    
                    chrgTrgt.push_back(MM2);
                }
                
            };
            
            case QMSCHEMEZ1:
            {
                // Selects all MM1 atoms
                chrgTrgt.push_back(MM1);
            } break;
        }
        
        
        qmMMChargeTarget[qmBndIt] = new int[chrgTrgt.size()] ;
        qmMMNumTargs[qmBndIt] =  chrgTrgt.size();
        
        DebugM(3, "MM-QM bond " << qmBndIt << "; MM atom " 
        << qmMMBond[qmBndIt][0] << " conections: \n" );
        
        for (size_t i=0; i < chrgTrgt.size(); i++) {
            qmMMChargeTarget[qmBndIt][i] = chrgTrgt[i];
            DebugM(3,"MM Bonded to: " << chrgTrgt[i] << "\n" );
        }
        
        chrgTrgt.clear();
    }
    
    if (bondsWithAtomLocal != NULL)
        delete [] bondsWithAtomLocal;  bondsWithAtomLocal = 0;
    if (byAtomSizeLocal != NULL)
        delete [] byAtomSizeLocal;  byAtomSizeLocal = 0;
    if (tmpArenaLocal != NULL)
        delete tmpArenaLocal;  tmpArenaLocal = 0;
    
    
    
    
    if(simParams->qmLSSOn) {
        
        std::map<Real,int> grpLSSSize ;
        std::map<Real,int>::iterator itGrpSize;
        
        qmLSSTotalNumAtms = 0;
        qmLSSResidueSize = 0;
        
        if (simParams->qmLSSFreq == 0)
            qmLSSFreq = simParams->stepsPerCycle ;
        else
            qmLSSFreq = simParams->qmLSSFreq;
            
        #ifdef DEBUG_QM
        int resSize = -1;
        #endif
        
        std::map<Real, int> grpNumLSSRes;
        std::map<Real, int>::iterator itGrpNumRes;
        
        for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
            
            if (it->atmIDs.size() != it->size) {
                iout << iERROR << "The number of atoms loaded for residue "
                << it->resID << " does not match the expected for this residue type.\n" 
                << endi;
                NAMD_die("Error parsing data for LSS.");
            }
            
            qmLSSTotalNumAtms += it->size;
            
            #ifdef DEBUG_QM
            if (resSize < 0) resSize = it->size ;
            if (resSize > 0 and resSize != it->size) {
                iout << iERROR << "The number of atoms loaded for residue "
                << it->resID << " does not match previously loaded residues.\n" 
                << endi;
                NAMD_die("Error parsing data for LSS.");
            }
                
//             DebugM(3,"Residue " << it->resID << ": " << it->segName
//             << " - from " << it->begAtmID << " with size "
//             << it->size << " (QM ID: " << it->qmGrpID
//             << ") has " << it->atmIDs.size() << " atoms: \n" ) ;
//             for (int i=0; i<it->atmIDs.size(); i++) 
//                 DebugM(3, it->atmIDs[i] << "\n" );
            #endif
            
            // Calculating total number of atoms per group
            itGrpSize = grpLSSSize.find(it->qmGrpID) ;
            if (itGrpSize != grpLSSSize.end())
                itGrpSize->second += it->size;
            else
                grpLSSSize.insert(std::pair<Real,int>(it->qmGrpID, it->size));
            
            // Calculating total number of solvent residues per group
            itGrpNumRes = grpNumLSSRes.find(it->qmGrpID) ;
            if (itGrpNumRes != grpNumLSSRes.end())
                itGrpNumRes->second += 1;
            else
                grpNumLSSRes.insert(std::pair<Real,int>(it->qmGrpID, 1));
        }
        
        qmLSSResidueSize = lssGrpRes.begin()->size ;
        
        qmLSSSize = new int[qmNumGrps];
        
        qmLSSIdxs = new int[qmLSSTotalNumAtms];
        int lssAtmIndx = 0;
        
        switch (simParams->qmLSSMode) {
        
        case QMLSSMODECOM:
        {
            
            qmLSSRefSize = new int[qmNumGrps];
            
            qmLSSMass = new Mass[qmLSSTotalNumAtms];
            
            qmLSSRefIDs = new int[totalNumRefAtms];
            qmLSSRefMass = new Mass[totalNumRefAtms];
            int lssRefIndx = 0;
            
            for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
                
                itGrpSize = grpNumLSSRes.find(qmGrpID[grpIndx]) ;
                
                if (itGrpSize != grpNumLSSRes.end())
                    qmLSSSize[grpIndx] =  itGrpSize->second;
                else
                    qmLSSSize[grpIndx] =  0;
                
                for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
                    
                    if (it->qmGrpID == qmGrpID[grpIndx]) {
                        for (int i=0; i<it->atmIDs.size(); i++) {
                            qmLSSIdxs[lssAtmIndx] = it->atmIDs[i];
                            qmLSSMass[lssAtmIndx] = atoms[it->atmIDs[i]].mass;
                            lssAtmIndx++;
                        }
                    }
                }
                
                DebugM(4, "QM group " << qmGrpID[grpIndx]
                << " has " << qmLSSSize[grpIndx] << " solvent molecules, "
                << "totalling " << grpLSSSize[qmGrpID[grpIndx]]
                << " atoms (check " << lssAtmIndx << ").\n" );
                
                qmLSSRefSize[grpIndx] = lssGrpRefIDs[qmGrpID[grpIndx]].size();
                for(int i=0; i < qmLSSRefSize[grpIndx]; i++) {
                    qmLSSRefIDs[lssRefIndx] = lssGrpRefIDs[qmGrpID[grpIndx]][i];
                    qmLSSRefMass[lssRefIndx] =  atoms[qmLSSRefIDs[lssRefIndx]].mass;
                    lssRefIndx++;
                }
                
                DebugM(4, "QM group " << qmGrpID[grpIndx]
                << " has " << qmLSSRefSize[grpIndx] << " non-solvent atoms for LSS.\n" );
            }
            
        } break ;
        
        case QMLSSMODEDIST:
        {
            for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
                
                itGrpSize = grpNumLSSRes.find(qmGrpID[grpIndx]) ;
                
                if (itGrpSize != grpNumLSSRes.end())
                    qmLSSSize[grpIndx] =  itGrpSize->second;
                else
                    qmLSSSize[grpIndx] =  0;
                
                for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
                    
                    if (it->qmGrpID == qmGrpID[grpIndx]) {
                        for (int i=0; i<it->atmIDs.size(); i++) {
                            qmLSSIdxs[lssAtmIndx] = it->atmIDs[i];
                            lssAtmIndx++;
                        }
                    }
                }
                
                DebugM(4, "QM group " << qmGrpID[grpIndx]
                << " has " << qmLSSSize[grpIndx] << " solvent molecules, "
                << "totalling " << grpLSSSize[qmGrpID[grpIndx]]
                << " atoms (check " << lssAtmIndx << ").\n" );
            }
            
        } break ;
            
        }
    }
    
    
    
    
    PDB *customPCPDB;
    
    // In case we have a custom and fixed set of point charges for each QM group,
    // we process the files containing information.
    current = NULL;
    if (simParams->qmCustomPCSel) {
        current = cfgList->find("QMCustomPCFile");
    }
    
    std::map<Real,std::vector<int> > qmPCVecMap ;
    
    int numParsedPBDs = 0;
    for ( ; current != NULL; current = current->next ) {
        
        iout << iINFO << "Parsing QM Custom PC file " << current->data << "\n" << endi;
        customPCPDB = new PDB(current->data);
        
        if (customPCPDB->num_atoms() != numAtoms)
           NAMD_die("Number of atoms in QM Custom PC PDB file doesn't match coordinate PDB");
        
        std::vector< int > currPCSel ;
        Real currQMID = 0 ;
        int currGrpSize = 0 ;
        
        for (int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
            
            BigReal beta = customPCPDB->atom(atmInd)->temperaturefactor() ;
            BigReal occ = customPCPDB->atom(atmInd)->occupancy() ;
            
            if ( beta != 0 && occ != 0)
                NAMD_die("An atom cannot be marked as QM and poitn charge simultaneously!");
            
            // If this is not a QM atom and 
            if (occ != 0) {
                currPCSel.push_back(atmInd) ;
            }
            
            if (beta != 0) {
                if (pdbP->atom(atmInd)->temperaturefactor() != beta)
                    NAMD_die("QM Group selection is different in reference PDB and Custom-PC PDB!");
                
                if (currQMID == 0) {
                    // If this is the first QM atom we find, keep the QM Group ID.
                    currQMID = beta;
                }
                else {
                    // For all other atoms, check if it is the same group. It must be!!
                    if (currQMID != beta)
                        NAMD_die("Found two different QM group IDs in this file!");
                }
                
                currGrpSize++;
            }
            
        }
        
        if (currGrpSize != qmGrpSizeVec[ qmGrpIDMap[currQMID] ])
            NAMD_die("Reference PDB and Custom-PC PDB have different QM group sizes!") ;
        
        qmPCVecMap.insert(std::pair<Real,std::vector<int> >(
            currQMID, currPCSel ));
        
        numParsedPBDs++;
        delete customPCPDB;
    }
    
    delete pdbP;
    
    if (numParsedPBDs != qmNumGrps && simParams->qmCustomPCSel) {
        iout << iWARN << "The number of files provided for custom point "
        "charges is not equal to the number of QM groups!\n" << endi;
    }
    
    // Initializes an array with the number of Custom Point Charges per 
    // QM group.
    qmCustPCSizes = new int[qmNumGrps];
    for (int i=0; i<qmNumGrps; i++)
        qmCustPCSizes[i] = 0;
    
    qmTotCustPCs = 0;
    
    // Stores the size of each Custom PC vector in the array.
    // We may not have PCs for all QM groups.
    for (auto it = qmPCVecMap.begin(); it != qmPCVecMap.end(); it++) {
        qmTotCustPCs += (*it).second.size();
        int qmIndx = qmGrpIDMap[(*it).first];
        qmCustPCSizes[qmIndx] = (*it).second.size();
    }
    
    qmCustomPCIdxs = new int[qmTotCustPCs];
    
    if (simParams->qmCustomPCSel) {
        
        int auxIter = 0;
        for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
            
            Real currQMID = qmGrpID[grpIndx];
            
            iout << iINFO << "Loading " << qmPCVecMap[currQMID].size()
            << " custom point charges to QM Group " << grpIndx
            << " (ID: " << currQMID << ")\n" << endi;
            
            for (int pcIndx=0; pcIndx<qmPCVecMap[currQMID].size(); pcIndx++) {
                qmCustomPCIdxs[auxIter] = qmPCVecMap[currQMID][pcIndx] ;
                auxIter++;
            }
        }
    } 
    
    // Conditional SMD option
    qmcSMD = simParams->qmCSMD;
    if (qmcSMD) {
        read_qm_csdm_file(qmGrpIDMap);
    }
    
    
    if (simParams->qmElecEmbed) {
        // Modifies Atom charges for Electrostatic Embedding.
        // QM atoms cannot have charges in the standard location, to keep
        // NAMD from calculating electrostatic interactions between QM and MM atoms.
        // We handle electrostatics ourselves in ComputeQM.C and in special
        // modifications for PME.
        for (int i=0; i<qmNumQMAtoms; i++) {
            qmAtmChrg[i] = atoms[qmAtmIndx[i]].charge;
            atoms[qmAtmIndx[i]].charge = 0;
        }
    }
    
    
    if ( simParams->extraBondsOn) {
        // Lifted from Molecule::build_extra_bonds
        
        StringList *file = cfgList->find("extraBondsFile");
        
        char err_msg[512];
        int a1,a2; float k, ref;
        
        for ( ; file; file = file->next ) {  // loop over files
            FILE *f = fopen(file->data,"r");
//             if ( ! f ) {
//               sprintf(err_msg, "UNABLE TO OPEN EXTRA BONDS FILE %s", file->data);
//               NAMD_err(err_msg);
//             } else {
//               iout << iINFO << "READING EXTRA BONDS FILE " << file->data <<"\n"<<endi;
//             }
            
            while ( 1 ) {
              char buffer[512];
              int ret_code;
              do {
                ret_code = NAMD_read_line(f, buffer);
              } while ( (ret_code==0) && (NAMD_blank_string(buffer)) );
              if (ret_code!=0) break;

              char type[512];
              sscanf(buffer,"%s",type);
              
              int badline = 0;
              if ( ! strncasecmp(type,"bond",4) ) {
                if ( sscanf(buffer, "%s %d %d %f %f %s",
                    type, &a1, &a2, &k, &ref, err_msg) != 5 ) badline = 1;
                
                // If an extra bond is defined between QM atoms, we make
                // note so that it wont be deleted when we delete bonded 
                // interactions between QM atoms.
                if( qmAtomGroup[a1] > 0 && qmAtomGroup[a2]) {
                    Bond tmp;
                    tmp.bond_type = 0;
                    tmp.atom1 = a1;  tmp.atom2 = a2;
                    qmExtraBonds.add(tmp);
                }
                
                
              }else if ( ! strncasecmp(type,"#",1) ) {
                continue;  // comment
              }
              
            }
            fclose(f);
        }
    }
    
    return;
    
#endif // MEM_OPT_VERSION
}

void Molecule::print_atoms

(

Parameters *

params

)

Definition at line 5337 of file Molecule.C.

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

Referenced by NamdState::loadStructure().

{
#ifdef MEM_OPT_VERSION
    DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
#else   
  register int i;
  Real sigma;
  Real epsilon;
  Real sigma14;
  Real epsilon14;

  DebugM(2,"ATOM LIST\n" \
      << "******************************************\n" \
                  << "NUM  NAME TYPE RES  MASS    CHARGE CHARGE   FEP-CHARGE"  \
      << "SIGMA   EPSILON SIGMA14 EPSILON14\n" \
        << endi);

  for (i=0; i<numAtoms; i++)
  {
    params->get_vdw_params(&sigma, &epsilon, &sigma14, &epsilon14, 
        atoms[i].vdw_type);

    DebugM(2,i+1 << " " << atomNames[i].atomname  \
              << " " << atomNames[i].atomtype << " " \
              << atomNames[i].resname  << " " << atoms[i].mass  \
        << " " << atoms[i].charge << " " << sigma \
        << " " << epsilon << " " << sigma14 \
        << " " << epsilon14 << "\n" \
        << endi);
  }
#endif  
}

void Molecule::print_bonds

(

Parameters *

params

)

Definition at line 5380 of file Molecule.C.

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

Referenced by NamdState::loadStructure().

{
#ifdef MEM_OPT_VERSION
    DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
#else   
  register int i;
  Real k;
  Real x0;

  DebugM(2,"BOND LIST\n" << "********************************\n" \
      << "ATOM1 ATOM2 TYPE1 TYPE2      k        x0" \
      << endi);

  for (i=0; i<numBonds; i++)
  {
    params->get_bond_params(&k, &x0, bonds[i].bond_type);

    DebugM(2,bonds[i].atom1+1 << " " \
       << bonds[i].atom2+1 << " "   \
       << atomNames[bonds[i].atom1].atomtype << " "  \
       << atomNames[bonds[i].atom2].atomtype << " " << k \
       << " " << x0 << endi);
  }
  
#endif  
}

void Molecule::print_exclusions

(

)

Definition at line 5417 of file Molecule.C.

References DebugM, and endi().

Referenced by NamdState::loadStructure().

{
#ifdef MEM_OPT_VERSION
    DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
#else
  register int i;

  DebugM(2,"EXPLICIT EXCLUSION LIST\n" \
      << "********************************\n" \
            << "ATOM1 ATOM2 " \
      << endi);

  for (i=0; i<numExclusions; i++)
  {
    DebugM(2,exclusions[i].atom1+1 << "  " \
       << exclusions[i].atom2+1 << endi);
  }
#endif
}

void Molecule::print_go_params

(

)

Definition at line 548 of file GoMolecule.C.

References DebugM, go_array, MAX_GO_CHAINS, and NumGoChains.

{
  int i;
  int j;
  int index;

  DebugM(3,NumGoChains << " Go PARAMETERS 3\n" \
         << "*****************************************" << std::endl);

  for (i=0; i<NumGoChains; i++) {
    for (j=0; j<NumGoChains; j++) {
      index = (i * MAX_GO_CHAINS) + j;
      //  Real epsilon;    // Epsilon
      //  Real exp_a;      // First exponent for attractive L-J term
      //  Real exp_b;      // Second exponent for attractive L-J term
      //  Real sigmaRep;   // Sigma for repulsive term
      //  Real epsilonRep; // Epsilon for replusive term
      DebugM(3,"Go index=(" << i << "," << j << ") epsilon=" << go_array[index].epsilon \
             << " exp_a=" << go_array[index].exp_a << " exp_b=" << go_array[index].exp_b \
             << " exp_rep=" << go_array[index].exp_rep << " sigmaRep=" \
             << go_array[index].sigmaRep << " epsilonRep=" << go_array[index].epsilonRep \
             << " cutoff=" << go_array[index].cutoff << std::endl);
    }
  }

}

void Molecule::print_go_sigmas

(

)

Definition at line 1134 of file GoMolecule.C.

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

{
  int i;  //  Counter
  int j;  //  Counter
  Real sigma;

  DebugM(3,"GO SIGMA ARRAY\n" \
         << "***************************" << std::endl);
  DebugM(3, "numGoAtoms: " << numGoAtoms << std::endl);

  if (goSigmaIndices == NULL) {
    DebugM(3, "GO SIGMAS HAVE NOT BEEN BUILT" << std::endl);
    return;
  }

  for (i=0; i<numAtoms; i++) {
    for (j=0; j<numAtoms; j++) {
      if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 ) {
        //DebugM(3, "i: " << i << ", j: " << j << std::endl);
        sigma = goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]];
        if (sigma > 0.0) {
          DebugM(3, "(" << i << "," << j << ") - +" << sigma << " ");
        }
        else {
          DebugM(3, "(" << i << "," << j << ") - " << sigma << " ");
        }
      } else {
        //DebugM(3, "0 ");
      }
    }
    if ( goSigmaIndices[i] != -1 ) {
      DebugM(3, "-----------" << std::endl);
    }
  }
  return;
}

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

inline

Definition at line 1312 of file Molecule.h.

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

void Molecule::read_alch_unpert_angles

(

FILE *

fd

)

Definition at line 1783 of file Molecule.C.

References NAMD_die(), and NAMD_read_int().

                                               {
  int atom_nums[3];  //  Atom numbers for the three atoms
  register int j;    //  Loop counter
  int num_read=0;    //  Number of angles read so far

  alch_unpert_angles=new Angle[num_alch_unpert_Angles];

  if (alch_unpert_angles == NULL) {
    NAMD_die("memory allocation failed in Molecule::read_alch_unpert_angles");
  }

  while (num_read < num_alch_unpert_Angles) {
    for (j=0; j<3; j++) {
      atom_nums[j]=NAMD_read_int(fd, "ANGLES")-1;

      if (atom_nums[j] >= numAtoms) {
        char err_msg[128];
        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);
        NAMD_die(err_msg);
      }
    }

    alch_unpert_angles[num_read].atom1=atom_nums[0];
    alch_unpert_angles[num_read].atom2=atom_nums[1];
    alch_unpert_angles[num_read].atom3=atom_nums[2];

    ++num_read;
  }
  return;
}

void Molecule::read_alch_unpert_bonds

(

FILE *

fd

)

Definition at line 1660 of file Molecule.C.

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

                                              {
  int atom_nums[2];  // Atom indexes for the bonded atoms
  register int j;      // Loop counter
  int num_read=0;    // Number of bonds read so far

  alch_unpert_bonds=new Bond[num_alch_unpert_Bonds];

  if (alch_unpert_bonds == NULL) {
    NAMD_die("memory allocations failed in Molecule::read_alch_unpert_bonds");
  }

  while (num_read < num_alch_unpert_Bonds) {
    for (j=0; j<2; j++) {
      atom_nums[j]=NAMD_read_int(fd, "BONDS")-1;

      if (atom_nums[j] >= numAtoms) {
        char err_msg[128];

        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);
        NAMD_die(err_msg);
      }
    }

    Bond *b = &(alch_unpert_bonds[num_read]);
    b->atom1=atom_nums[0];
    b->atom2=atom_nums[1];

    ++num_read;
  }
  return;
}

void Molecule::read_alch_unpert_dihedrals

(

FILE *

fd

)

Definition at line 1931 of file Molecule.C.

References NAMD_die(), and NAMD_read_int().

                                                  {
  int atom_nums[4];  // The 4 atom indexes
  int num_read=0;    // number of dihedrals read so far
  register int j;    // loop counter

  alch_unpert_dihedrals = new Dihedral[num_alch_unpert_Dihedrals];

  if (alch_unpert_dihedrals == NULL) {
    NAMD_die("memory allocation failed in Molecule::read_alch_unpert_dihedrals");
  }

  while (num_read < num_alch_unpert_Dihedrals) {
    for (j=0; j<4; j++) {
      atom_nums[j]=NAMD_read_int(fd, "DIHEDRALS")-1;

      if (atom_nums[j] >= numAtoms) {
        char err_msg[128];

        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);
        NAMD_die(err_msg);
      }
    }

    alch_unpert_dihedrals[num_read].atom1=atom_nums[0];
    alch_unpert_dihedrals[num_read].atom2=atom_nums[1];
    alch_unpert_dihedrals[num_read].atom3=atom_nums[2];
    alch_unpert_dihedrals[num_read].atom4=atom_nums[3];

    num_read++;
  }
  return;
}

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

{

  int i;                   // Counter
  int j;                   // Counter
  int  par_type=0;         //  What type of parameter are we currently
                           //  dealing with? (vide infra)
  // JLai -- uncommented
  int  skipline;           //  skip this line?
  int  skipall = 0;        //  skip rest of file;
  char buffer[512];           //  Buffer to store each line of the file
  char first_word[512];           //  First word of the current line
  int read_count = 0;      //  Count of input parameters on a given line
  int chain1 = 0;          //  First chain type for pair interaction
  int chain2 = 0;          //  Second chain type for pair interaction
  int int1;                //  First parameter int
  int int2;                //  Second parameter int
  Real r1;                 //  Parameter Real
  char in2[512];           //  Second parameter word
  GoValue *goValue1 = NULL;    //  current GoValue for loading parameters
  GoValue *goValue2 = NULL;    //  other current GoValue for loading parameters
  Bool sameGoChain = FALSE;    //  whether the current GoValue is within a chain or between chains
  int restrictionCount = 0;    //  number of Go restrictions set for a given chain pair
  FILE *pfile;                 //  File descriptor for the parameter file

  /*  Check to make sure that we haven't previously been told     */
  /*  that all the files were read                                */
  /*if (AllFilesRead)
    {
    NAMD_die("Tried to read another parameter file after being told that all files were read!");
    }*/
  
  /*  Initialize go_indices  */
  for (i=0; i<MAX_GO_CHAINS+1; i++) {
    go_indices[i] = -1;
  }

  /*  Initialize go_array   */
  for (i=0; i<MAX_GO_CHAINS*MAX_GO_CHAINS; i++) {
    go_array[i].epsilon = 1.25;
    go_array[i].exp_a = 12;
    go_array[i].exp_b = 6;
    go_array[i].exp_rep = 12;
    go_array[i].sigmaRep = 2.25;
    go_array[i].epsilonRep = 0.03;
    go_array[i].cutoff = 4.0;
    for (j=0; j<MAX_RESTRICTIONS; j++) {
      go_array[i].restrictions[j] = -1;
    }
  }

  /*  Try and open the file                                        */
  if ( (pfile = fopen(fname, "r")) == NULL)
    {
      char err_msg[256];
      
      sprintf(err_msg, "UNABLE TO OPEN GO PARAMETER FILE %s\n", fname);
      NAMD_die(err_msg);
    }
  
  /*  Keep reading in lines until we hit the EOF                        */
  while (NAMD_read_line(pfile, buffer) != -1)
    {
      /*  Get the first word of the line                        */
      NAMD_find_first_word(buffer, first_word);
      skipline=0;
      
      /*  First, screen out things that we ignore.                   */   
      /*  blank lines, lines that start with '!' or '*', lines that  */
      /*  start with "END".                                          */
      if (!NAMD_blank_string(buffer) &&
          (strncmp(first_word, "!", 1) != 0) &&
          (strncmp(first_word, "*", 1) != 0) &&
          (strncasecmp(first_word, "END", 3) != 0))
        {
          if ( skipall ) {
            iout << iWARN << "SKIPPING PART OF GO PARAMETER FILE AFTER RETURN STATEMENT\n" << endi;
            break;
          }
          /*  Now, determine the apropriate parameter type.   */
          if (strncasecmp(first_word, "chaintypes", 10)==0)
            {
              read_count=sscanf(buffer, "%s %d %d\n", first_word, &int1, &int2);
              if (read_count != 3) {
                char err_msg[512];
                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);
                NAMD_die(err_msg);
              }
              chain1 = int1;
              chain2 = int2;
              if (chain1 < 1 || chain1 > MAX_GO_CHAINS ||
                  chain2 < 1 || chain2 > MAX_GO_CHAINS) {
                char err_msg[512];
                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);
                NAMD_die(err_msg);
              }
              if (go_indices[chain1] == -1) {
                go_indices[chain1] = NumGoChains;
                NumGoChains++;
              }
              if (go_indices[chain2] == -1) {
                go_indices[chain2] = NumGoChains;
                NumGoChains++;
              }
              if (chain1 == chain2) {
                sameGoChain = TRUE;
              } else {
                sameGoChain = FALSE;
              }
              //goValue = &go_array[(chain1 * MAX_GO_CHAINS) + chain2];
              goValue1 = &go_array[(chain1*MAX_GO_CHAINS) + chain2];
              goValue2 = &go_array[(chain2*MAX_GO_CHAINS) + chain1];
#if CODE_REDUNDANT
              goValue1 = &go_array[(go_indices[chain1]*MAX_GO_CHAINS) + go_indices[chain2]];
              goValue2 = &go_array[(go_indices[chain2]*MAX_GO_CHAINS) + go_indices[chain1]];
#endif
              restrictionCount = 0;    //  restrictionCount applies to each chain pair separately
            }
          else if (strncasecmp(first_word, "epsilonRep", 10)==0)
            {
              read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
              if (read_count != 2) {}
              goValue1->epsilonRep = r1;
              if (!sameGoChain) {
                goValue2->epsilonRep = r1;
              }
            }
          else if (strncasecmp(first_word, "epsilon", 7)==0)
            {
              // Read in epsilon
              read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
              if (read_count != 2) {}
              goValue1->epsilon = r1;
              if (!sameGoChain) {
                goValue2->epsilon = r1;
              }
            }
          else if (strncasecmp(first_word, "exp_a", 5)==0)
            {
              read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
              if (read_count != 2) {}
              goValue1->exp_a = int1;
              if (!sameGoChain) {
                goValue2->exp_a = int1;
              }
            }
          else if (strncasecmp(first_word, "exp_b", 5)==0)
            {
              read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
              if (read_count != 2) {}
              goValue1->exp_b = int1;
              if (!sameGoChain) {
                goValue2->exp_b = int1;
              }
            }
          else if (strncasecmp(first_word, "exp_rep", 5)==0)
            {
              read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
              if (read_count != 2) {}
              goValue1->exp_rep = int1;
              if (!sameGoChain) {
                goValue2->exp_rep = int1;
              }
            }
          else if (strncasecmp(first_word, "exp_Rep", 5)==0)
            {
              read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
              if (read_count != 2) {}
              goValue1->exp_rep = int1;
              if (!sameGoChain) {
              goValue2->exp_rep = int1;
              }
            }
          else if (strncasecmp(first_word, "sigmaRep", 8)==0)
            {
              read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
              if (read_count != 2) {}
              goValue1->sigmaRep = r1;
              if (!sameGoChain) {
                goValue2->sigmaRep = r1;
              }
            }
          else if (strncasecmp(first_word, "cutoff", 6)==0)
            {
              read_count=sscanf(buffer, "%s %f\n", first_word, &r1);
              if (read_count != 2) {}
              goValue1->cutoff = r1;
              if (!sameGoChain) {
                goValue2->cutoff = r1;
              }
            }
          else if (strncasecmp(first_word, "restriction", 10)==0)
            {
              read_count=sscanf(buffer, "%s %d\n", first_word, &int1);
              if (read_count != 2) {}
              if (int1 < 0) {
                DebugM(3, "ERROR: residue restriction value must be nonnegative.  int1=" << int1 << "\n");
              }
              if (!sameGoChain) {
                //goValue2->restrictions[restrictionCount] = int1;
                DebugM(3, "ERROR: residue restrictions should not be defined between two separate chains.  chain1=" << chain1 << ", chain2=" << chain2 << "\n");
              }
              else {
                goValue1->restrictions[restrictionCount] = int1;
              }
              restrictionCount++;
            }
          else if (strncasecmp(first_word, "return", 4)==0)
            {
              skipall=8;
              skipline=1;
            }        
          else // if (par_type == 0)
            {
              /*  This is an unknown paramter.        */
              /*  This is BAD                                */
              char err_msg[512];
              
              sprintf(err_msg, "UNKNOWN PARAMETER IN GO PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
              NAMD_die(err_msg);
            }
        }
      else
        {
          skipline=1;
        }
    }
  
  /*  Close the file  */
  fclose(pfile);
  
  return;
}

void Molecule::read_parm

(

Ambertoppar *

)

void Molecule::read_parm

(

AmberParm7Reader::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, MAX_GO_CHAINS, MAX_RESTRICTIONS, NAMD_die(), numAtoms, numGoAtoms, NumGoChains, pointerToGoBeg, pointerToGoEnd, go_val::restrictions, and go_val::sigmaRep.

                                               {
      // Ported by JLai -- Original by JE
      // JE - receive Go info
      Real *a1, *a2, *a3, *a4;
      int *i1, *i2, *i3, *i4;
      int maxGoChainsSqr = MAX_GO_CHAINS*MAX_GO_CHAINS;  // JE JLai Go code
      msg->get(NumGoChains);
      
      if (NumGoChains) {
        //go_indices = new int[MAX_GO_CHAINS+1];
        //go_array = new GoValue[MAX_GO_CHAINS*MAX_GO_CHAINS];
        
        //      int go_indices[MAX_GO_CHAINS+1];        //  Indices from chainIDs to go_array      
        //      GoValue go_array[MAX_GO_CHAINS*MAX_GO_CHAINS];   //  Array of Go params
        msg->get(MAX_GO_CHAINS+1,go_indices);
        
        a1 = new Real[maxGoChainsSqr];
        a2 = new Real[maxGoChainsSqr];
        a3 = new Real[maxGoChainsSqr];
        a4 = new Real[maxGoChainsSqr];
        i1 = new int[maxGoChainsSqr];
        i2 = new int[maxGoChainsSqr];
        i3 = new int[maxGoChainsSqr];
        i4 = new int[maxGoChainsSqr*MAX_RESTRICTIONS];
        
        if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) || 
             (i1 == NULL) || (i2 == NULL) || (i3 == NULL) || (i4 == NULL) )
          {
            NAMD_die("memory allocation failed in Molecule::send_Molecule");
          }
        
        msg->get(maxGoChainsSqr, a1);
        msg->get(maxGoChainsSqr, a2);
        msg->get(maxGoChainsSqr, a3);
        msg->get(maxGoChainsSqr, a4);
        msg->get(maxGoChainsSqr, i1);
        msg->get(maxGoChainsSqr, i2);
        msg->get(maxGoChainsSqr, i3);
        msg->get(maxGoChainsSqr*MAX_RESTRICTIONS, i4);
        
        for (int i=0; i<maxGoChainsSqr; i++) {
          go_array[i].epsilon = a1[i];
          go_array[i].sigmaRep = a2[i];
          go_array[i].epsilonRep = a3[i];
          go_array[i].cutoff = a4[i];
          go_array[i].exp_a = i1[i];
          go_array[i].exp_b = i2[i];
          go_array[i].exp_rep = i3[i];
          for (int j=0; j<MAX_RESTRICTIONS; j++) {
            go_array[i].restrictions[j] = i4[i*MAX_RESTRICTIONS + j];
          }
        }
        
        delete [] a1;
        delete [] a2;
        delete [] a3;
        delete [] a4;
        delete [] i1;
        delete [] i2;
        delete [] i3;
        delete [] i4;

        //msg->get(MAX_GO_CHAINS*MAX_GO_CHAINS, (char*)go_array);
        
        /*DebugM(3,"NumGoChains:" << NumGoChains << std::endl);
          for (int ii=0; ii<MAX_GO_CHAINS; ii++) {
          for (int jj=0; jj<MAX_GO_CHAINS; jj++) {
          DebugM(3,"go_array[" << ii << "][" << jj << "]:" << go_array[ii][jj] << std::endl);
          }
          }
          for (int ii=0; ii<MAX_GO_CHAINS+1; ii++) {
          DebugM(3,"go_indices[" << ii << "]:" << go_indices[ii] << std::endl);
          }*/
      }

      if (simParams->goForcesOn) {
        switch(simParams->goMethod) {
        case 1:
          msg->get(numGoAtoms);
          //printf("Deleting goSigmaIndiciesA\n");
          delete [] goSigmaIndices;
          goSigmaIndices = new int32[numAtoms];
          //printf("Deleting atomChainTypesA\n");
          delete [] atomChainTypes;
          atomChainTypes = new int32[numGoAtoms];
          //printf("Deleting goSigmasA\n");
          delete [] goSigmas;
          goSigmas = new Real[numGoAtoms*numGoAtoms];
          //printf("Deleting goWithinCutoffA\n"); 
          delete [] goWithinCutoff;
          goWithinCutoff = new bool[numGoAtoms*numGoAtoms];
          msg->get(numAtoms, goSigmaIndices);
          msg->get(numGoAtoms, atomChainTypes);
          msg->get(numGoAtoms*numGoAtoms, goSigmas);
          msg->get(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
          break;          
        case 2: //GSR
          msg->get(numGoAtoms);
          delete [] pointerToGoBeg;
          pointerToGoBeg = new int[numAtoms];
          msg->get(numAtoms,pointerToGoBeg);
          delete [] pointerToGoEnd;
          pointerToGoEnd = new int[numAtoms];
          msg->get(numAtoms,pointerToGoEnd);
          delete [] goSigmaIndices;
          goSigmaIndices = new int32[numAtoms];
          msg->get(numAtoms,goSigmaIndices);
          delete [] goResidIndices;
          goResidIndices = new int32[numAtoms];
          msg->get(numAtoms,goResidIndices);      
          delete [] atomChainTypes;
          atomChainTypes = new int32[numGoAtoms];
          msg->get(numGoAtoms,atomChainTypes);
          msg->get(goNumLJPair);
          delete [] goIndxLJA;
          goIndxLJA = new int[goNumLJPair];
          msg->get(goNumLJPair,goIndxLJA);
          delete [] goIndxLJB;
          goIndxLJB = new int[goNumLJPair];
          msg->get(goNumLJPair,goIndxLJB);
          delete [] goSigmaPairA;
          goSigmaPairA = new double[goNumLJPair];
          msg->get(goNumLJPair,goSigmaPairA);
          delete [] goSigmaPairB;
          goSigmaPairB = new double[goNumLJPair];
          msg->get(goNumLJPair,goSigmaPairB);  
          break;
        case 3:
          msg->get(numGoAtoms);
          //printf("Deleting goSigmaIndiciesB\n");
          delete [] goSigmaIndices;
          goSigmaIndices = new int32[numAtoms];
          //printf("Deleting atomChainTypesB\n");
          delete [] atomChainTypes;
          atomChainTypes = new int32[numGoAtoms];
          //delete [] goSigmas;
          //goSigmas = new Real[numGoAtoms*numGoAtoms];
          //delete [] goWithinCutoff;
          //goWithinCutoff = new bool[numGoAtoms*numGoAtoms];
          //printf("Deleting goCoordinatesB\n");
          delete [] goCoordinates;
          goCoordinates = new Real[numGoAtoms*3];
          //printf("Deleting goResidsB\n");
          delete [] goResids;
          goResids = new int[numGoAtoms];
          msg->get(numAtoms, goSigmaIndices);
          msg->get(numGoAtoms, atomChainTypes);
          //msg->get(numGoAtoms*numGoAtoms, goSigmas);
          //msg->get(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
          msg->get(numGoAtoms*3, goCoordinates);
          msg->get(numGoAtoms, goResids);
          break;
        }
      }

      delete msg;

}

void Molecule::receive_Molecule

(

MIStream *

msg

)

Definition at line 5806 of file Molecule.C.

References SimParameters::alchOn, gromacsPair::atom1, gromacsPair::atom2, atomNamePool, atoms, atomSigPool, SimParameters::consForceOn, SimParameters::consTorqueOn, SimParameters::constraintsOn, DebugM, SimParameters::excludeFromPressure, exclusions, SimParameters::fixedAtomsOn, MIStream::get(), Parameters::get_num_vdw_params(), SimParameters::goGroPair, gromacsPair::gromacsPair_type, SimParameters::langevinOn, SimParameters::LCPOOn, SimParameters::lesOn, SimParameters::mgridforceOn, SimParameters::movDragOn, numGaussPair, numLJPair, gromacsPair::pairC12, gromacsPair::pairC6, SimParameters::pairInteractionOn, SimParameters::qmForcesOn, SimParameters::rotDragOn, SimParameters::sdScaling, SimParameters::soluteScalingOn, SimParameters::stirOn, SimParameters::tCoupleOn, and GridforceGrid::unpack_grid().

Referenced by Node::resendMolecule().

                                            {
  //  Get the atom information
  msg->get(numAtoms);

#ifdef MEM_OPT_VERSION
//in the memory optimized version, only the atom signatures are recved
//from the master Node. --Chao Mei

  msg->get(massPoolSize);
  if(atomMassPool) delete [] atomMassPool;
  atomMassPool = new Real[massPoolSize];
  msg->get(massPoolSize, atomMassPool);

  msg->get(chargePoolSize);
  if(atomChargePool) delete [] atomChargePool;
  atomChargePool = new Real[chargePoolSize];
  msg->get(chargePoolSize, atomChargePool);

  //get atoms' signatures
  msg->get(atomSigPoolSize);
  if(atomSigPool) delete [] atomSigPool;
  atomSigPool = new AtomSignature[atomSigPoolSize];
  for(int i=0; i<atomSigPoolSize; i++)
      atomSigPool[i].unpack(msg);

  //get exclusions' signatures
  msg->get(exclSigPoolSize);
  if(exclSigPool) delete [] exclSigPool;
  exclSigPool = new ExclusionSignature[exclSigPoolSize];
  for(int i=0; i<exclSigPoolSize; i++)
      exclSigPool[i].unpack(msg);
 
  msg->get(numHydrogenGroups);      
  msg->get(maxHydrogenGroupSize);      
  msg->get(numMigrationGroups);      
  msg->get(maxMigrationGroupSize);      
  msg->get(isOccupancyValid);
  msg->get(isBFactorValid);

   //get names for atoms
  msg->get(atomNamePoolSize);
  atomNamePool = new char *[atomNamePoolSize];
  for(int i=0; i<atomNamePoolSize;i++) {
    int len;
    msg->get(len);
    atomNamePool[i] = nameArena->getNewArray(len+1);
    msg->get(len, atomNamePool[i]);
  }
  
  if(simParams->fixedAtomsOn){
    int numFixedAtomsSet;
    msg->get(numFixedAtoms);
    msg->get(numFixedAtomsSet);
    fixedAtomsSet = new AtomSetList(numFixedAtomsSet);
    msg->get(numFixedAtomsSet*sizeof(AtomSet), (char *)(fixedAtomsSet->begin()));
  } 

  if(simParams->constraintsOn){
    int numConstrainedAtomsSet;
    msg->get(numConstraints);
    msg->get(numConstrainedAtomsSet);
    constrainedAtomsSet = new AtomSetList(numConstrainedAtomsSet);
    msg->get(numConstrainedAtomsSet*sizeof(AtomSet), (char *)(constrainedAtomsSet->begin()));
  } 

#else
  delete [] atoms;
  atoms= new Atom[numAtoms];  
  msg->get(numAtoms*sizeof(Atom), (char*)atoms);

  //  Get the bond information
  msg->get(numRealBonds);
  msg->get(numBonds);    
  if (numBonds)
  {
    delete [] bonds;
    bonds=new Bond[numBonds]; 
    msg->get(numBonds*sizeof(Bond), (char*)bonds);
  }  

  //  Get the angle information
  msg->get(numAngles);  
  if (numAngles)
  {
    delete [] angles;
    angles=new Angle[numAngles];  
    msg->get(numAngles*sizeof(Angle), (char*)angles);
  }  

  //  Get the dihedral information
  msg->get(numDihedrals);    
  if (numDihedrals)
  {
    delete [] dihedrals;
    dihedrals=new Dihedral[numDihedrals];  
    msg->get(numDihedrals*sizeof(Dihedral), (char*)dihedrals);
  }  

  if (simParams->sdScaling) {
    msg->get(num_alch_unpert_Bonds);
    alch_unpert_bonds=new Bond[num_alch_unpert_Bonds];
    msg->get(num_alch_unpert_Bonds*sizeof(Bond), (char*)alch_unpert_bonds);

    msg->get(num_alch_unpert_Angles);
    alch_unpert_angles=new Angle[num_alch_unpert_Angles];
    msg->get(num_alch_unpert_Angles*sizeof(Angle), (char*)alch_unpert_angles);

    msg->get(num_alch_unpert_Dihedrals);
    alch_unpert_dihedrals=new Dihedral[num_alch_unpert_Dihedrals];
    msg->get(num_alch_unpert_Dihedrals*sizeof(Dihedral), (char*)alch_unpert_dihedrals);
  }
  
  //  Get the improper information
  msg->get(numImpropers);
  if (numImpropers)
  {
    delete [] impropers;
    impropers=new Improper[numImpropers];  
    msg->get(numImpropers*sizeof(Improper), (char*)impropers);
  }
  
  //  Get the crossterm information
  msg->get(numCrossterms);
  if (numCrossterms)
  {
    delete [] crossterms;
    crossterms=new Crossterm[numCrossterms];  
    msg->get(numCrossterms*sizeof(Crossterm), (char*)crossterms);
  }
  
  //  Get the hydrogen bond donors
  msg->get(numDonors);  
  if (numDonors)
  {
    delete [] donors;
    donors=new Bond[numDonors];  
    msg->get(numDonors*sizeof(Bond), (char*)donors);
  }
  
  //  Get the hydrogen bond acceptors
  msg->get(numAcceptors);  
  if (numAcceptors)
  {
    delete [] acceptors;
    acceptors=new Bond[numAcceptors];  
    msg->get(numAcceptors*sizeof(Bond), (char*)acceptors);
  }
  
  //  Get the exclusion information 
  msg->get(numExclusions);  
  if (numExclusions)
  {
    delete [] exclusions;
    exclusions=new Exclusion[numExclusions];  
    msg->get(numExclusions*sizeof(Exclusion), (char*)exclusions);
  }
        
      //  Get the constraint information, if they are active
      if (simParams->constraintsOn)
      {
         msg->get(numConstraints);

         delete [] consIndexes;
         consIndexes = new int32[numAtoms];
         
         msg->get(numAtoms, consIndexes);
         
         if (numConstraints)
         {
           delete [] consParams;
           consParams = new ConstraintParams[numConstraints];
      
           msg->get(numConstraints*sizeof(ConstraintParams), (char*)consParams);
         }
      }
#endif

      /* BEGIN gf */
      if (simParams->mgridforceOn)
      {
         DebugM(3, "Receiving gridforce info\n");
         
         msg->get(numGridforceGrids);
         
         DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n");
         
         delete [] numGridforces;
         numGridforces = new int[numGridforceGrids];
         
         delete [] gridfrcIndexes;      // Should I be deleting elements of these first?
         delete [] gridfrcParams;
         delete [] gridfrcGrid;
         gridfrcIndexes = new int32*[numGridforceGrids];
         gridfrcParams = new GridforceParams*[numGridforceGrids];
         gridfrcGrid = new GridforceGrid*[numGridforceGrids];
         
         int grandTotalGrids = 0;
         for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
             msg->get(numGridforces[gridnum]);
             
             gridfrcIndexes[gridnum] = new int32[numAtoms];
             msg->get(numAtoms, gridfrcIndexes[gridnum]);
         
             if (numGridforces[gridnum])
             {
                 gridfrcParams[gridnum] = new GridforceParams[numGridforces[gridnum]];
                 msg->get(numGridforces[gridnum]*sizeof(GridforceParams), (char*)gridfrcParams[gridnum]);
             }
             
             gridfrcGrid[gridnum] = GridforceGrid::unpack_grid(gridnum, msg);
             
             grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
         }
      }
      /* END gf */
      
      //  Get the stirring information, if stirring is  active
      if (simParams->stirOn)
      {
         msg->get(numStirredAtoms);

         delete [] stirIndexes;
         stirIndexes = new int32[numAtoms];
         
         msg->get(numAtoms, stirIndexes);
         
         if (numStirredAtoms)
         {
           delete [] stirParams;
           stirParams = new StirParams[numStirredAtoms];
      
           msg->get(numStirredAtoms*sizeof(StirParams), (char*)stirParams);
         }
      }
      
      //  Get the moving drag information, if it is active
      if (simParams->movDragOn) {
         msg->get(numMovDrag);
         delete [] movDragIndexes;
         movDragIndexes = new int32[numAtoms];
         msg->get(numAtoms, movDragIndexes);
         if (numMovDrag)
         {
           delete [] movDragParams;
           movDragParams = new MovDragParams[numMovDrag];
           msg->get(numMovDrag*sizeof(MovDragParams), (char*)movDragParams);
         }
      }
      
      //  Get the rotating drag information, if it is active
      if (simParams->rotDragOn) {
         msg->get(numRotDrag);
         delete [] rotDragIndexes;
         rotDragIndexes = new int32[numAtoms];
         msg->get(numAtoms, rotDragIndexes);
         if (numRotDrag)
         {
           delete [] rotDragParams;
           rotDragParams = new RotDragParams[numRotDrag];
           msg->get(numRotDrag*sizeof(RotDragParams), (char*)rotDragParams);
         }
      }
      
      //  Get the "constant" torque information, if it is active
      if (simParams->consTorqueOn) {
         msg->get(numConsTorque);
         delete [] consTorqueIndexes;
         consTorqueIndexes = new int32[numAtoms];
         msg->get(numAtoms, consTorqueIndexes);
         if (numConsTorque)
         {
           delete [] consTorqueParams;
           consTorqueParams = new ConsTorqueParams[numConsTorque];
           msg->get(numConsTorque*sizeof(ConsTorqueParams), (char*)consTorqueParams);
         }
      }
      
      // Get the constant force information, if it's active
      if (simParams->consForceOn)
      { msg->get(numConsForce);
        delete [] consForceIndexes;
        consForceIndexes = new int32[numAtoms];
        msg->get(numAtoms, consForceIndexes);
        if (numConsForce)
        { delete [] consForce;
          consForce = new Vector[numConsForce];
          msg->get(numConsForce*sizeof(Vector), (char*)consForce);
        }
      }

      if (simParams->excludeFromPressure) {
        exPressureAtomFlags = new int32[numAtoms];
        msg->get(numExPressureAtoms);
        msg->get(numAtoms, exPressureAtomFlags);
      }

#ifndef MEM_OPT_VERSION
      //  Get the langevin parameters, if they are active
      if (simParams->langevinOn || simParams->tCoupleOn)
      {
        delete [] langevinParams;
        langevinParams = new Real[numAtoms];

        msg->get(numAtoms, langevinParams);
      }

      //  Get the fixed atoms, if they are active
      if (simParams->fixedAtomsOn)
      {
        delete [] fixedAtomFlags;
        fixedAtomFlags = new int32[numAtoms];

        msg->get(numFixedAtoms);
        msg->get(numAtoms, fixedAtomFlags);
        msg->get(numFixedRigidBonds);
      }

      if (simParams->qmForcesOn)
      {
        if( qmAtomGroup != 0)
            delete [] qmAtomGroup;
        qmAtomGroup = new Real[numAtoms];
        
        msg->get(numAtoms, qmAtomGroup);
        
        msg->get(qmNumQMAtoms);
        
        if( qmAtmChrg != 0)
            delete [] qmAtmChrg;
        qmAtmChrg = new Real[qmNumQMAtoms];
        
        msg->get(qmNumQMAtoms, qmAtmChrg);
        
        if( qmAtmIndx != 0)
            delete [] qmAtmIndx;
        qmAtmIndx = new int[qmNumQMAtoms];
        
        msg->get(qmNumQMAtoms, qmAtmIndx);
        
        msg->get(qmNoPC);
        
        msg->get(qmNumBonds);
        
        msg->get(qmMeNumBonds);
        
        if( qmMeMMindx != 0)
            delete [] qmMeMMindx;
        qmMeMMindx = new int[qmMeNumBonds];
        
        msg->get(qmMeNumBonds, qmMeMMindx);
        
        if( qmMeQMGrp != 0)
            delete [] qmMeQMGrp;
        qmMeQMGrp = new Real[qmMeNumBonds];
        
        msg->get(qmMeNumBonds, qmMeQMGrp);
        
        msg->get(qmPCFreq);
        
        msg->get(qmNumGrps);
        
        if( qmGrpID != 0)
            delete [] qmGrpID;
        qmGrpID = new Real[qmNumGrps];
        msg->get(qmNumGrps, qmGrpID);
        
        if( qmCustPCSizes != 0)
            delete [] qmCustPCSizes;
        qmCustPCSizes = new int[qmNumGrps];
        msg->get(qmNumGrps, qmCustPCSizes);
        
        msg->get(qmTotCustPCs);
        
        if( qmCustomPCIdxs != 0)
            delete [] qmCustomPCIdxs;
        qmCustomPCIdxs = new int[qmTotCustPCs];
        msg->get(qmTotCustPCs, qmCustomPCIdxs);
      }
    
//fepb
      //receive fep atom info
      if (simParams->alchOn || simParams->lesOn || simParams->pairInteractionOn) {
        delete [] fepAtomFlags;
        fepAtomFlags = new unsigned char[numAtoms];

        msg->get(numFepInitial);
        msg->get(numFepFinal);
        msg->get(numAtoms*sizeof(unsigned char), (char*)fepAtomFlags);
      }
//fepe

//soluteScaling
      if (simParams->soluteScalingOn) {
        delete [] ssAtomFlags;
        delete [] ss_vdw_type;
        delete [] ss_index;
        ssAtomFlags = new unsigned char[numAtoms];
        ss_vdw_type = new int [params->get_num_vdw_params()];
        ss_index = new int [numAtoms];
        msg->get(numAtoms*sizeof(unsigned char), (char*)ssAtomFlags);
        msg->get(ss_num_vdw_params);
        msg->get(params->get_num_vdw_params()*sizeof(int), (char*)ss_vdw_type);
        msg->get(numAtoms*sizeof(int), (char*)ss_index);
      }
//soluteScaling
#ifdef OPENATOM_VERSION
      // This needs to be refactored into its own version
      if (simParams->openatomOn) {
        delete [] fepAtomFlags;
        fepAtomFlags = new unsigned char[numAtoms];

        msg->get(numFepInitial);
        msg->get(numAtoms*sizeof(unsigned char), (char*)fepAtomFlags);
#endif //OPENATOM_VERSION

      // DRUDE: receive data read from PSF
      msg->get(is_lonepairs_psf);
      if (is_lonepairs_psf) {
        msg->get(numLphosts);
        delete[] lphosts;
        lphosts = new Lphost[numLphosts];
        msg->get(numLphosts*sizeof(Lphost), (char*)lphosts);
      }
      msg->get(is_drude_psf);
      if (is_drude_psf) {
        delete[] drudeConsts;
        drudeConsts = new DrudeConst[numAtoms];
        msg->get(numAtoms*sizeof(DrudeConst), (char*)drudeConsts);
        msg->get(numTholes);
        delete[] tholes;
        tholes = new Thole[numTholes];
        msg->get(numTholes*sizeof(Thole), (char*)tholes);
        msg->get(numAnisos);
        delete[] anisos;
        anisos = new Aniso[numAnisos];
        msg->get(numAnisos*sizeof(Aniso), (char*)anisos);
      }
      msg->get(numZeroMassAtoms);
      // DRUDE

  //LCPO
  if (simParams->LCPOOn) {
    delete [] lcpoParamType;
    lcpoParamType = new int[numAtoms];
    msg->get(numAtoms, (int*)lcpoParamType);
  }

  //Receive GromacsPairStuff -- JLai

  if (simParams->goGroPair) {
    msg->get(numLJPair);
    delete [] indxLJA;
    indxLJA = new int[numLJPair];
    msg->get(numLJPair,indxLJA);
    delete [] indxLJB;
    indxLJB = new int[numLJPair];
    msg->get(numLJPair,indxLJB);
    delete [] pairC6;
    pairC6 = new Real[numLJPair];    
    msg->get(numLJPair,pairC6);
    delete [] pairC12;
    pairC12 = new Real[numLJPair];
    msg->get(numLJPair,pairC12);
    delete [] gromacsPair_type;
    gromacsPair_type = new int[numLJPair];
    msg->get(numLJPair,gromacsPair_type);
    delete [] pointerToLJBeg;
    pointerToLJBeg = new int[numAtoms];
    msg->get((numAtoms),pointerToLJBeg);
    delete [] pointerToLJEnd;
    pointerToLJEnd = new int[numAtoms];
    msg->get((numAtoms),pointerToLJEnd);
    // JLai
    delete [] gromacsPair;
    gromacsPair = new GromacsPair[numLJPair];
    for(int i=0; i < numLJPair; i++) {
        gromacsPair[i].atom1 = indxLJA[i];
        gromacsPair[i].atom2 = indxLJB[i];
        gromacsPair[i].pairC6  = pairC6[i];
        gromacsPair[i].pairC12 = pairC12[i];
        gromacsPair[i].gromacsPair_type = gromacsPair_type[i];
    }
    //
    msg->get(numGaussPair);
    delete [] indxGaussA;
    indxGaussA = new int[numGaussPair];
    msg->get(numGaussPair,indxGaussA);
    delete [] indxGaussB;
    indxGaussB = new int[numGaussPair];
    msg->get(numGaussPair,indxGaussB);
    delete [] gA;
    gA = new Real[numGaussPair];
    msg->get(numGaussPair,gA);
    delete [] gMu1;
    gMu1 = new Real[numGaussPair];
    msg->get(numGaussPair,gMu1);
    delete [] giSigma1;
    giSigma1 = new Real[numGaussPair];
    msg->get(numGaussPair,giSigma1);
    delete [] gMu2;
    gMu2 = new Real[numGaussPair];
    msg->get(numGaussPair,gMu2);
    delete [] giSigma2;
    giSigma2 = new Real[numGaussPair];
    msg->get(numGaussPair,giSigma2);
    delete [] gRepulsive;
    gRepulsive = new Real[numGaussPair];
    msg->get(numGaussPair,gRepulsive);
    delete [] pointerToGaussBeg;
    pointerToGaussBeg = new int[numAtoms];
    msg->get((numAtoms),pointerToGaussBeg);
    delete [] pointerToGaussEnd;
    pointerToGaussEnd = new int[numAtoms];
    msg->get((numAtoms),pointerToGaussEnd);
    //
  }
#endif

      //  Now free the message 
      delete msg;

#ifdef MEM_OPT_VERSION

      build_excl_check_signatures();

    //set num{Calc}Tuples(Bonds,...,Impropers) to 0
    numBonds = numCalcBonds = 0;
    numAngles = numCalcAngles = 0;
    numDihedrals = numCalcDihedrals = 0;
    numImpropers = numCalcImpropers = 0;
    numCrossterms = numCalcCrossterms = 0;
    numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;  
    // JLai
    numLJPair = numCalcLJPair = 0;
    // End of JLai

#else

      //  analyze the data and find the status of each atom
      build_atom_status();
      build_lists_by_atom();      

      
#endif
}

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

Referenced by Node::reloadCharges().

Real Molecule::rigid_bond_length ( int  atomnum ) const

inline

Definition at line 1457 of file Molecule.h.

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

  {
    return(rigidBondLengths[atomnum]);
  }

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, MAX_GO_CHAINS, MAX_RESTRICTIONS, NAMD_die(), numAtoms, numGoAtoms, NumGoChains, pointerToGoBeg, pointerToGoEnd, MOStream::put(), go_val::restrictions, and go_val::sigmaRep.

                                            {
  Real *a1, *a2, *a3, *a4;
  int *i1, *i2, *i3, *i4;
  int maxGoChainsSqr = MAX_GO_CHAINS*MAX_GO_CHAINS;  // JE JLai Go code
  msg->put(NumGoChains);
  
  if (NumGoChains) {
    //      int go_indices[MAX_GO_CHAINS+1];        //  Indices from chainIDs to go_array      
    //      GoValue go_array[MAX_GO_CHAINS*MAX_GO_CHAINS];   //  Array of Go params
    msg->put(MAX_GO_CHAINS+1,go_indices);

    a1 = new Real[maxGoChainsSqr];
    a2 = new Real[maxGoChainsSqr];
    a3 = new Real[maxGoChainsSqr];
    a4 = new Real[maxGoChainsSqr];
    i1 = new int[maxGoChainsSqr];
    i2 = new int[maxGoChainsSqr];
    i3 = new int[maxGoChainsSqr];
    i4 = new int[maxGoChainsSqr*MAX_RESTRICTIONS];

    if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) || 
         (i1 == NULL) || (i2 == NULL) || (i3 == NULL) || (i4 == NULL) )
    {
      NAMD_die("memory allocation failed in Molecules::send_Molecules");
    }

    for (int i=0; i<maxGoChainsSqr; i++) {
      a1[i] = go_array[i].epsilon;
      a2[i] = go_array[i].sigmaRep;
      a3[i] = go_array[i].epsilonRep;
      a4[i] = go_array[i].cutoff;
      i1[i] = go_array[i].exp_a;
      i2[i] = go_array[i].exp_b;
      i3[i] = go_array[i].exp_rep;
      for (int j=0; j<MAX_RESTRICTIONS; j++) {
        i4[i*MAX_RESTRICTIONS + j] = go_array[i].restrictions[j];
      }
    }

    msg->put(maxGoChainsSqr, a1);
    msg->put(maxGoChainsSqr, a2);
    msg->put(maxGoChainsSqr, a3);
    msg->put(maxGoChainsSqr, a4);
    msg->put(maxGoChainsSqr, i1);
    msg->put(maxGoChainsSqr, i2);
    msg->put(maxGoChainsSqr, i3);
    msg->put(maxGoChainsSqr*MAX_RESTRICTIONS, i4);

    delete [] a1;
    delete [] a2;
    delete [] a3;
    delete [] a4;
    delete [] i1;
    delete [] i2;
    delete [] i3;
    delete [] i4;
  }

  //Ported JLai
  if (simParams->goForcesOn) {
    switch(simParams->goMethod) {
    case 1:
      msg->put(numGoAtoms);
      msg->put(numAtoms, goSigmaIndices);
      msg->put(numGoAtoms, atomChainTypes);
      msg->put(numGoAtoms*numGoAtoms, goSigmas);
      msg->put(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
      // printf("Molecule.C sending atomChainTypes %d %d \n", numGoAtoms, atomChainTypes);
      break;
    case 2: //GSS
      msg->put(numGoAtoms);
      msg->put(numAtoms,pointerToGoBeg);
      msg->put(numAtoms,pointerToGoEnd);
      msg->put(numAtoms,goSigmaIndices);
      msg->put(numAtoms,goResidIndices);
      msg->put(numGoAtoms,atomChainTypes);
      msg->put(goNumLJPair);
      msg->put(goNumLJPair,goIndxLJA);
      msg->put(goNumLJPair,goIndxLJB);
      msg->put(goNumLJPair,goSigmaPairA);
      msg->put(goNumLJPair,goSigmaPairB);
      break;
    case 3:
      msg->put(numGoAtoms);
      msg->put(numAtoms, goSigmaIndices);
      msg->put(numGoAtoms, atomChainTypes);
      //msg->put(numGoAtoms*numGoAtoms, goSigmas);
      //msg->put(numGoAtoms*numGoAtoms*sizeof(bool), (char*)goWithinCutoff);
      msg->put(numGoAtoms*3, goCoordinates);
      msg->put(numGoAtoms, goResids);
      break;
    }
  } 

  msg->end();
  delete msg;
}

void Molecule::send_Molecule

(

MOStream *

msg

)

Definition at line 5448 of file Molecule.C.

References SimParameters::alchOn, atomNamePool, atomSigPool, SimParameters::consForceOn, SimParameters::consTorqueOn, SimParameters::constraintsOn, DebugM, MOStream::end(), endi(), SimParameters::excludeFromPressure, SimParameters::fixedAtomsOn, Parameters::get_num_vdw_params(), SimParameters::goGroPair, SimParameters::langevinOn, SimParameters::LCPOOn, SimParameters::lesOn, SimParameters::mgridforceOn, SimParameters::movDragOn, numGaussPair, GridforceGrid::pack_grid(), SimParameters::pairInteractionOn, MOStream::put(), SimParameters::qmForcesOn, SimParameters::rotDragOn, SimParameters::sdScaling, SimParameters::soluteScalingOn, SimParameters::stirOn, and SimParameters::tCoupleOn.

Referenced by Node::resendMolecule().

                                         {
#ifdef MEM_OPT_VERSION
//in the memory optimized version, only the atom signatures are broadcast
//to other Nodes. --Chao Mei

  msg->put(numAtoms);

  msg->put(massPoolSize);
  msg->put(massPoolSize, atomMassPool);

  msg->put(chargePoolSize);
  msg->put(chargePoolSize, atomChargePool); 

  //put atoms' signatures
  msg->put(atomSigPoolSize);
  for(int i=0; i<atomSigPoolSize; i++)
      atomSigPool[i].pack(msg);

  //put atom's exclusion signatures
  msg->put(exclSigPoolSize);
  for(int i=0; i<exclSigPoolSize; i++)
      exclSigPool[i].pack(msg);

  msg->put(numHydrogenGroups);      
  msg->put(maxHydrogenGroupSize);      
  msg->put(numMigrationGroups);      
  msg->put(maxMigrationGroupSize);            
  msg->put(isOccupancyValid);
  msg->put(isBFactorValid);
  
  //put names for atoms
  msg->put(atomNamePoolSize);
  for(int i=0; i<atomNamePoolSize;i++) {
    int len = strlen(atomNamePool[i]);
    msg->put(len);
    msg->put(len*sizeof(char), atomNamePool[i]);
  } 
  
  if(simParams->fixedAtomsOn){
    int numFixedAtomsSet = fixedAtomsSet->size();
    msg->put(numFixedAtoms);
    msg->put(numFixedAtomsSet);
    msg->put(numFixedAtomsSet*sizeof(AtomSet), (char *)(fixedAtomsSet->begin()));
  }

  if (simParams->constraintsOn) {
    int numConstrainedAtomsSet = constrainedAtomsSet->size();
    msg->put(numConstraints);
    msg->put(numConstrainedAtomsSet);
    msg->put(numConstrainedAtomsSet*sizeof(AtomSet), (char *)(constrainedAtomsSet->begin()));
  }
    
#else
  msg->put(numAtoms);
  msg->put(numAtoms*sizeof(Atom), (char*)atoms);
  
  //  Send the bond information
  msg->put(numRealBonds);
  msg->put(numBonds);
 
  if (numBonds)
  {
    msg->put(numBonds*sizeof(Bond), (char*)bonds);
  }

  //  Send the angle information
  msg->put(numAngles);  
  if (numAngles)
  {
    msg->put(numAngles*sizeof(Angle), (char*)angles);
  }  

  //  Send the dihedral information
  msg->put(numDihedrals);
  if (numDihedrals)
  {
    msg->put(numDihedrals*sizeof(Dihedral), (char*)dihedrals);
  }  

  if (simParams->sdScaling) {
    msg->put(num_alch_unpert_Bonds);
    msg->put(num_alch_unpert_Bonds*sizeof(Bond), (char*)alch_unpert_bonds);

    msg->put(num_alch_unpert_Angles);
    msg->put(num_alch_unpert_Angles*sizeof(Angle), (char*)alch_unpert_angles);

    msg->put(num_alch_unpert_Dihedrals);
    msg->put(num_alch_unpert_Dihedrals*sizeof(Dihedral), (char*)alch_unpert_dihedrals);
  }

  //  Send the improper information
  msg->put(numImpropers);  
  if (numImpropers)
  {
    msg->put(numImpropers*sizeof(Improper), (char*)impropers);
  }

  //  Send the crossterm information
  msg->put(numCrossterms);
  if (numCrossterms)
  {
    msg->put(numCrossterms*sizeof(Crossterm), (char*)crossterms);
  }

  // send the hydrogen bond donor information
  msg->put(numDonors);
  if(numDonors)
  {
    msg->put(numDonors*sizeof(Bond), (char*)donors);
  }

  // send the hydrogen bond acceptor information
  msg->put(numAcceptors);
  if(numAcceptors)
  {
    msg->put(numAcceptors*sizeof(Bond), (char*)acceptors);
  }

  //  Send the exclusion information  
  msg->put(numExclusions);
  if (numExclusions)
  {
    msg->put(numExclusions*sizeof(Exclusion), (char*)exclusions);
  }      
  //  Send the constraint information, if used
  if (simParams->constraintsOn)
  {
     msg->put(numConstraints);
     
     msg->put(numAtoms, consIndexes);
     
     if (numConstraints)
     {
       msg->put(numConstraints*sizeof(ConstraintParams), (char*)consParams);
     }
  }
#endif
  
  /* BEGIN gf */
  // Send the gridforce information, if used
  if (simParams->mgridforceOn)
  {
    DebugM(3, "Sending gridforce info\n" << endi);
    msg->put(numGridforceGrids);
    
    for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
      msg->put(numGridforces[gridnum]);
      msg->put(numAtoms, gridfrcIndexes[gridnum]);
      if (numGridforces[gridnum])
      {
       msg->put(numGridforces[gridnum]*sizeof(GridforceParams), (char*)gridfrcParams[gridnum]);
      }
      GridforceGrid::pack_grid(gridfrcGrid[gridnum], msg);
    }
  }
  /* END gf */
  
  //  Send the stirring information, if used
  if (simParams->stirOn)
  {
     //CkPrintf ("DEBUG: putting numStirredAtoms..\n");
     msg->put(numStirredAtoms);
     //CkPrintf ("DEBUG: putting numAtoms,stirIndexes.. numAtoms=%d\n",numStirredAtoms);
     msg->put(numAtoms, stirIndexes);
     //CkPrintf ("DEBUG: if numStirredAtoms..\n");
     if (numStirredAtoms)
     {
       //CkPrintf ("DEBUG: big put, with (char*)stirParams\n");
       msg->put(numStirredAtoms*sizeof(StirParams), (char*)stirParams);
     }
  }
  
  
  //  Send the moving drag information, if used
  if (simParams->movDragOn) {
     msg->put(numMovDrag);
     msg->put(numAtoms, movDragIndexes);
     if (numMovDrag)
     {
       msg->put(numMovDrag*sizeof(MovDragParams), (char*)movDragParams);
     }
  }
  
  //  Send the rotating drag information, if used
  if (simParams->rotDragOn) {
     msg->put(numRotDrag);
     msg->put(numAtoms, rotDragIndexes);
     if (numRotDrag)
     {
       msg->put(numRotDrag*sizeof(RotDragParams), (char*)rotDragParams);
     }
  }
  
  //  Send the "constant" torque information, if used
  if (simParams->consTorqueOn) {
     msg->put(numConsTorque);
     msg->put(numAtoms, consTorqueIndexes);
     if (numConsTorque)
     {
       msg->put(numConsTorque*sizeof(ConsTorqueParams), (char*)consTorqueParams);
     }
  }
  
  // Send the constant force information, if used
  if (simParams->consForceOn)
  { msg->put(numConsForce);
    msg->put(numAtoms, consForceIndexes);
    if (numConsForce)
      msg->put(numConsForce*sizeof(Vector), (char*)consForce);
  }
  
  if (simParams->excludeFromPressure) {
    msg->put(numExPressureAtoms);
    msg->put(numAtoms, exPressureAtomFlags);
  }
  
#ifndef MEM_OPT_VERSION
  //  Send the langevin parameters, if active
  if (simParams->langevinOn || simParams->tCoupleOn)
  {
    msg->put(numAtoms, langevinParams);
  }
  
  //  Send fixed atoms, if active
  if (simParams->fixedAtomsOn)
  {
    msg->put(numFixedAtoms);
    msg->put(numAtoms, fixedAtomFlags);
  msg->put(numFixedRigidBonds);
  }
  
  if (simParams->qmForcesOn)
  {
    msg->put(numAtoms, qmAtomGroup);
    msg->put(qmNumQMAtoms);
    msg->put(qmNumQMAtoms, qmAtmChrg);
    msg->put(qmNumQMAtoms, qmAtmIndx);
    msg->put(qmNoPC);
    msg->put(qmNumBonds);
    msg->put(qmMeNumBonds);
    msg->put(qmMeNumBonds, qmMeMMindx);
    msg->put(qmMeNumBonds, qmMeQMGrp);
    msg->put(qmPCFreq);
    msg->put(qmNumGrps);
    msg->put(qmNumGrps, qmGrpID);
    msg->put(qmNumGrps, qmCustPCSizes);
    msg->put(qmTotCustPCs);
    msg->put(qmTotCustPCs, qmCustomPCIdxs);
  }
  
  //fepb
  // send fep atom info
  if (simParams->alchOn || simParams->lesOn || simParams->pairInteractionOn) {
    msg->put(numFepInitial);
    msg->put(numFepFinal);
    msg->put(numAtoms*sizeof(char), (char*)fepAtomFlags);
  }
  //fepe

  if (simParams->soluteScalingOn) {
    msg->put(numAtoms*sizeof(char), (char*)ssAtomFlags);
    msg->put(ss_num_vdw_params);
    msg->put(params->get_num_vdw_params()*sizeof(int), (char*)ss_vdw_type);
    msg->put(numAtoms*sizeof(int), (char*)ss_index);
  }

  #ifdef OPENATOM_VERSION
  // needs to be refactored into its own openatom version
  if (simParams->openatomOn ) {
    msg->put(numFepInitial);
    msg->put(numAtoms*sizeof(char), (char*)fepAtomFlags);
  }
  #endif //OPENATOM_VERSION
  
  // DRUDE: send data read from PSF
  msg->put(is_lonepairs_psf);
  if (is_lonepairs_psf) {
    msg->put(numLphosts);
    msg->put(numLphosts*sizeof(Lphost), (char*)lphosts);
  }
  msg->put(is_drude_psf);
  if (is_drude_psf) {
    msg->put(numAtoms*sizeof(DrudeConst), (char*)drudeConsts);
    msg->put(numTholes);
    msg->put(numTholes*sizeof(Thole), (char*)tholes);
    msg->put(numAnisos);
    msg->put(numAnisos*sizeof(Aniso), (char*)anisos);
  }
  msg->put(numZeroMassAtoms);
  // DRUDE

  //LCPO
  if (simParams->LCPOOn) {
    msg->put(numAtoms, (int*)lcpoParamType);
  }
  
  //Send GromacsPairStuff -- JLai
  if (simParams->goGroPair) {
    msg->put(numLJPair);
    msg->put(numLJPair,indxLJA);
    msg->put(numLJPair,indxLJB);
    msg->put(numLJPair,pairC6);
    msg->put(numLJPair,pairC12);
    msg->put(numLJPair,gromacsPair_type);
    msg->put((numAtoms),pointerToLJBeg);
    msg->put((numAtoms),pointerToLJEnd);
    msg->put(numGaussPair);
    msg->put(numGaussPair,indxGaussA);
    msg->put(numGaussPair,indxGaussB);
    msg->put(numGaussPair,gA);
    msg->put(numGaussPair,gMu1);
    msg->put(numGaussPair,giSigma1);
    msg->put(numGaussPair,gMu2);
    msg->put(numGaussPair,giSigma2);
    msg->put(numGaussPair,gRepulsive);
    msg->put((numAtoms),pointerToGaussBeg);
    msg->put((numAtoms),pointerToGaussEnd);
  }
#endif

  // Broadcast the message to the other nodes
  msg->end();
  delete msg;

#ifdef MEM_OPT_VERSION

  build_excl_check_signatures();

  //set num{Calc}Tuples(Bonds,...,Impropers) to 0
  numBonds = numCalcBonds = 0;
  numAngles = numCalcAngles = 0;
  numDihedrals = numCalcDihedrals = 0;
  numImpropers = numCalcImpropers = 0;
  numCrossterms = numCalcCrossterms = 0;
  numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;  
  // JLai
  numLJPair = numCalcLJPair = 0;
  // End of JLai

#else

  //  Now build arrays of indexes into these arrays by atom      
  build_lists_by_atom();

#endif
}

int Molecule::set_gridfrc_grid ( int  gridnum, GridforceGrid *  grid  )

inline

Definition at line 1289 of file Molecule.h.

References numGridforceGrids.

Referenced by Node::reloadGridforceGrid().

  {
      if (grid && gridnum >= 0 && gridnum < numGridforceGrids) {
          gridfrcGrid[gridnum] = grid;
          return 0;
      } else {
          return -1;
      }
  }

void Molecule::set_qm_replaceAll ( Bool  newReplaceAll )

inline

Definition at line 818 of file Molecule.h.

Referenced by NamdState::loadStructure().

{ qmReplaceAll = newReplaceAll; } ;

void Molecule::setBFactorData

(

molfile_atom_t *

atomarray

)

Definition at line 3243 of file Molecule.C.

                                                      {
    bfactor = new float[numAtoms];
    for(int i=0; i<numAtoms; i++) {
        bfactor[i] = atomarray[i].bfactor;
    }
}

void Molecule::setOccupancyData

(

molfile_atom_t *

atomarray

)

Definition at line 3236 of file Molecule.C.

                                                        {
    occupancy = new float[numAtoms];
    for(int i=0; i<numAtoms; i++) {
        occupancy[i] = atomarray[i].occupancy;
    }
}

Friends And Related Function Documentation

friend class AngleElem

friend

Definition at line 216 of file Molecule.h.

friend class AnisoElem

friend

Definition at line 220 of file Molecule.h.

friend class BondElem

friend

Definition at line 215 of file Molecule.h.

friend class CrosstermElem

friend

Definition at line 221 of file Molecule.h.

friend class DihedralElem

friend

Definition at line 217 of file Molecule.h.

friend class ExclElem

friend

Definition at line 214 of file Molecule.h.

friend class GromacsPairElem

friend

Definition at line 223 of file Molecule.h.

friend class ImproperElem

friend

Definition at line 218 of file Molecule.h.

friend class TholeElem

friend

Definition at line 219 of file Molecule.h.

friend class WorkDistrib

friend

Definition at line 225 of file Molecule.h.

Member Data Documentation

Angle* Molecule::alch_unpert_angles

Definition at line 577 of file Molecule.h.

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

Bond* Molecule::alch_unpert_bonds

Definition at line 576 of file Molecule.h.

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

Dihedral* Molecule::alch_unpert_dihedrals

Definition at line 578 of file Molecule.h.

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

int Molecule::alchDroppedAngles

Definition at line 564 of file Molecule.h.

Referenced by NamdState::loadStructure().

int Molecule::alchDroppedDihedrals

Definition at line 565 of file Molecule.h.

Referenced by NamdState::loadStructure().

int Molecule::alchDroppedImpropers

Definition at line 566 of file Molecule.h.

Referenced by NamdState::loadStructure().

int32* Molecule::atomChainTypes

Definition at line 644 of file Molecule.h.

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

Vector* Molecule::consForce

Definition at line 614 of file Molecule.h.

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

int32* Molecule::consForceIndexes

Definition at line 613 of file Molecule.h.

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

int32* Molecule::consTorqueIndexes

Definition at line 616 of file Molecule.h.

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

ConsTorqueParams* Molecule::consTorqueParams

Definition at line 617 of file Molecule.h.

Referenced by get_constorque_params().

BigReal Molecule::energyNative

Definition at line 685 of file Molecule.h.

Referenced by build_go_arrays(), and goInit().

BigReal Molecule::energyNonnative

Definition at line 686 of file Molecule.h.

Referenced by build_go_arrays(), and goInit().

Real* Molecule::gA

Definition at line 677 of file Molecule.h.

Referenced by get_gro_force2().

Real* Molecule::giSigma1

Definition at line 679 of file Molecule.h.

Referenced by get_gro_force2().

Real* Molecule::giSigma2

Definition at line 681 of file Molecule.h.

Referenced by get_gro_force2().

Real* Molecule::gMu1

Definition at line 678 of file Molecule.h.

Referenced by get_gro_force2().

Real* Molecule::gMu2

Definition at line 680 of file Molecule.h.

Referenced by get_gro_force2().

GoValue Molecule::go_array[MAX_GO_CHAINS *MAX_GO_CHAINS]

Definition at line 1588 of file Molecule.h.

Referenced by 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(), print_go_params(), read_go_file(), receive_GoMolecule(), and send_GoMolecule().

int Molecule::go_indices[MAX_GO_CHAINS+1]

Definition at line 1589 of file Molecule.h.

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

Real* Molecule::goCoordinates

Definition at line 649 of file Molecule.h.

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

int* Molecule::goIndxLJA

Definition at line 654 of file Molecule.h.

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

int* Molecule::goIndxLJB

Definition at line 655 of file Molecule.h.

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

int Molecule::goNumLJPair

Definition at line 653 of file Molecule.h.

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

PDB* Molecule::goPDB

Definition at line 651 of file Molecule.h.

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

int32* Molecule::goResidIndices

Definition at line 646 of file Molecule.h.

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

int* Molecule::goResids

Definition at line 650 of file Molecule.h.

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

int32* Molecule::goSigmaIndices

Definition at line 645 of file Molecule.h.

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

double* Molecule::goSigmaPairA

Definition at line 656 of file Molecule.h.

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

double* Molecule::goSigmaPairB

Definition at line 657 of file Molecule.h.

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

Real* Molecule::goSigmas

Definition at line 647 of file Molecule.h.

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

bool* Molecule::goWithinCutoff

Definition at line 648 of file Molecule.h.

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

Real* Molecule::gRepulsive

Definition at line 682 of file Molecule.h.

Referenced by get_gro_force2().

int* Molecule::gromacsPair_type

Definition at line 670 of file Molecule.h.

HydrogenGroup Molecule::hydrogenGroup

Definition at line 640 of file Molecule.h.

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

int* Molecule::indxGaussA

Definition at line 675 of file Molecule.h.

int* Molecule::indxGaussB

Definition at line 676 of file Molecule.h.

Referenced by get_gro_force2().

int* Molecule::indxLJA

Definition at line 666 of file Molecule.h.

int* Molecule::indxLJB

Definition at line 667 of file Molecule.h.

Referenced by get_gro_force2().

int Molecule::is_drude_psf

Definition at line 462 of file Molecule.h.

int Molecule::is_lonepairs_psf

Definition at line 463 of file Molecule.h.

int Molecule::isBFactorValid

Definition at line 1469 of file Molecule.h.

int Molecule::isOccupancyValid

Definition at line 1469 of file Molecule.h.

int Molecule::maxHydrogenGroupSize

Definition at line 601 of file Molecule.h.

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

int Molecule::maxMigrationGroupSize

Definition at line 603 of file Molecule.h.

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

int Molecule::num_alch_unpert_Angles

Definition at line 574 of file Molecule.h.

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

int Molecule::num_alch_unpert_Bonds

Definition at line 573 of file Molecule.h.

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

int Molecule::num_alch_unpert_Dihedrals

Definition at line 575 of file Molecule.h.

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

int Molecule::numAcceptors

Definition at line 570 of file Molecule.h.

int Molecule::numAngles

Definition at line 561 of file Molecule.h.

Referenced by buildAngleData(), delete_qm_bonded(), dumpbench(), AngleElem::getMoleculePointers(), and NamdState::loadStructure().

int Molecule::numAnisos

Number of anisotropic terms.

Definition at line 584 of file Molecule.h.

Referenced by AnisoElem::getMoleculePointers().

int Molecule::numAtoms

Definition at line 557 of file Molecule.h.

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

int Molecule::numBonds

Definition at line 560 of file Molecule.h.

Referenced by buildBondData(), delete_qm_bonded(), dumpbench(), BondElem::getMoleculePointers(), NamdState::loadStructure(), and prepare_qm().

int Molecule::numCalcAngles

Definition at line 623 of file Molecule.h.

Referenced by Controller::compareChecksums(), and dumpbench().

int Molecule::numCalcAnisos

Definition at line 632 of file Molecule.h.

Referenced by Controller::compareChecksums().

int Molecule::numCalcBonds

Definition at line 622 of file Molecule.h.

Referenced by Controller::compareChecksums(), and dumpbench().

int Molecule::numCalcCrossterms

Definition at line 626 of file Molecule.h.

Referenced by Controller::compareChecksums().

int Molecule::numCalcDihedrals

Definition at line 624 of file Molecule.h.

Referenced by Controller::compareChecksums(), and dumpbench().

int64 Molecule::numCalcExclusions

Definition at line 627 of file Molecule.h.

Referenced by Controller::compareChecksums(), and dumpbench().

int64 Molecule::numCalcFullExclusions

Definition at line 628 of file Molecule.h.

Referenced by Controller::compareChecksums().

int Molecule::numCalcImpropers

Definition at line 625 of file Molecule.h.

Referenced by Controller::compareChecksums(), and dumpbench().

int Molecule::numCalcLJPair

Definition at line 663 of file Molecule.h.

int Molecule::numCalcTholes

Definition at line 631 of file Molecule.h.

Referenced by Controller::compareChecksums().

int Molecule::numConsForce

Definition at line 612 of file Molecule.h.

Referenced by NamdState::loadStructure().

int Molecule::numConsTorque

Definition at line 596 of file Molecule.h.

Referenced by is_atom_constorqued().

int Molecule::numConstraints

Definition at line 589 of file Molecule.h.

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

int Molecule::numCrossterms

Definition at line 568 of file Molecule.h.

Referenced by buildCrosstermData(), delete_qm_bonded(), CrosstermElem::getMoleculePointers(), and NamdState::loadStructure().

int Molecule::numDihedrals

Definition at line 562 of file Molecule.h.

Referenced by buildDihedralData(), delete_qm_bonded(), dumpbench(), DihedralElem::getMoleculePointers(), and NamdState::loadStructure().

int Molecule::numDonors

Definition at line 569 of file Molecule.h.

int Molecule::numDrudeAtoms

Number of Drude particles.

Definition at line 582 of file Molecule.h.

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

int Molecule::numExclusions

Definition at line 571 of file Molecule.h.

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

int Molecule::numExPressureAtoms

Definition at line 599 of file Molecule.h.

Referenced by is_atom_exPressure().

int Molecule::numFepFinal

Definition at line 609 of file Molecule.h.

Referenced by NamdState::loadStructure().

int Molecule::numFepInitial

Definition at line 608 of file Molecule.h.

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

int Molecule::numFixedAtoms

Definition at line 597 of file Molecule.h.

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

int Molecule::numFixedGroups

Definition at line 604 of file Molecule.h.

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

int Molecule::numFixedRigidBonds

Definition at line 606 of file Molecule.h.

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

int Molecule::numGaussPair

Definition at line 674 of file Molecule.h.

int Molecule::numGoAtoms

Definition at line 643 of file Molecule.h.

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

int Molecule::NumGoChains

Definition at line 1590 of file Molecule.h.

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

int Molecule::numGridforceGrids

Definition at line 591 of file Molecule.h.

Referenced by ComputeGridForce::doForce(), get_gridfrc_grid(), colvarproxy_namd::init_volmap(), is_atom_gridforced(), NamdState::loadStructure(), and set_gridfrc_grid().

int* Molecule::numGridforces

Definition at line 592 of file Molecule.h.

int Molecule::numHydrogenGroups

Definition at line 600 of file Molecule.h.

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

int Molecule::numImpropers

Definition at line 567 of file Molecule.h.

Referenced by buildImproperData(), delete_qm_bonded(), dumpbench(), ImproperElem::getMoleculePointers(), and NamdState::loadStructure().

int Molecule::numLJPair

Definition at line 662 of file Molecule.h.

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

int Molecule::numLonepairs

Number of lone pairs.

Definition at line 581 of file Molecule.h.

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

int Molecule::numLphosts

Number of lone pair host records in PSF.

Definition at line 585 of file Molecule.h.

int Molecule::numMigrationGroups

Definition at line 602 of file Molecule.h.

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

int Molecule::numMovDrag

Definition at line 594 of file Molecule.h.

Referenced by is_atom_movdragged().

int Molecule::numMultipleDihedrals

Definition at line 636 of file Molecule.h.

Referenced by NamdState::loadStructure().

int Molecule::numMultipleImpropers

Definition at line 638 of file Molecule.h.

Referenced by NamdState::loadStructure().

int Molecule::numPair

Definition at line 661 of file Molecule.h.

int Molecule::numRealBonds

Definition at line 559 of file Molecule.h.

Referenced by buildBondData(), delete_qm_bonded(), and prepare_qm().

int Molecule::numRigidBonds

Definition at line 605 of file Molecule.h.

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

int Molecule::numRotDrag

Definition at line 595 of file Molecule.h.

Referenced by is_atom_rotdragged().

int Molecule::numStirredAtoms

Definition at line 598 of file Molecule.h.

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

int Molecule::numTholes

Number of Thole terms.

Definition at line 583 of file Molecule.h.

Referenced by TholeElem::getMoleculePointers().

int64 Molecule::numTotalExclusions

Definition at line 572 of file Molecule.h.

int Molecule::numZeroMassAtoms

Number of atoms with zero mass.

Definition at line 586 of file Molecule.h.

Real* Molecule::pairC12

Definition at line 669 of file Molecule.h.

Referenced by get_gro_force2().

Real* Molecule::pairC6

Definition at line 668 of file Molecule.h.

Referenced by get_gro_force2().

int* Molecule::pointerToGaussBeg

Definition at line 672 of file Molecule.h.

Referenced by get_gro_force2().

int* Molecule::pointerToGaussEnd

Definition at line 673 of file Molecule.h.

Referenced by get_gro_force2().

int* Molecule::pointerToGoBeg

Definition at line 658 of file Molecule.h.

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

int* Molecule::pointerToGoEnd

Definition at line 659 of file Molecule.h.

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

int* Molecule::pointerToLJBeg

Definition at line 664 of file Molecule.h.

Referenced by get_gro_force2().

int* Molecule::pointerToLJEnd

Definition at line 665 of file Molecule.h.

Referenced by get_gro_force2().

Real Molecule::r_ohc

Definition at line 468 of file Molecule.h.

Real Molecule::r_om

Definition at line 467 of file Molecule.h.

int* Molecule::ss_index

Definition at line 459 of file Molecule.h.

int Molecule::ss_num_vdw_params

Definition at line 457 of file Molecule.h.

Referenced by LJTable::LJTable().

int* Molecule::ss_vdw_type

Definition at line 458 of file Molecule.h.

int Molecule::suspiciousAlchBonds

Definition at line 563 of file Molecule.h.

Referenced by NamdState::loadStructure().

BigReal Molecule::tail_corr_dUdl_1

Definition at line 473 of file Molecule.h.

BigReal Molecule::tail_corr_dUdl_2

Definition at line 473 of file Molecule.h.

BigReal Molecule::tail_corr_ener

Definition at line 471 of file Molecule.h.

Referenced by Controller::rescaleaccelMD().

BigReal Molecule::tail_corr_virial

Definition at line 472 of file Molecule.h.

BigReal Molecule::tail_corr_virial_1

Definition at line 474 of file Molecule.h.

BigReal Molecule::tail_corr_virial_2

Definition at line 474 of file Molecule.h.

---

The documentation for this class was generated from the following files:

• Molecule.h
• GoMolecule.C
• Molecule.C
• MoleculeQM.C