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 *, 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 150 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 *	,
		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 1048 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 1038 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 1058 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 1034 of file Molecule.h.

Referenced by NamdState::loadStructure().

{ delete [] bfactor; bfactor=NULL; }

void Molecule::freeOccupancyData ( )

inline

Definition at line 1030 of file Molecule.h.

Referenced by NamdState::loadStructure().

{ delete [] occupancy; occupancy=NULL; }

Bond* Molecule::get_acceptor ( int  dnum ) const

inline

Definition at line 1105 of file Molecule.h.

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

Angle* Molecule::get_angle ( int  anum ) const

inline

Definition at line 1068 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

int32* Molecule::get_angles_for_atom ( int  anum )

inline

Definition at line 1147 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

Referenced by atoms_1to4(), and dumpbench().

      { return bondsByAtom[anum]; } 

int Molecule::get_cluster ( int  anum ) const

inline

Definition at line 1023 of file Molecule.h.

Referenced by wrap_coor_int().

{ return cluster[anum]; }

int Molecule::get_clusterSize ( int  anum ) const

inline

Definition at line 1024 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 1262 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 1336 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 1077 of file Molecule.h.

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

int32* Molecule::get_crossterms_for_atom ( int  anum )

inline

Definition at line 1153 of file Molecule.h.

      { return crosstermsByAtom[anum]; }  

const Real* Molecule::get_cSMDcoffs ( )

inline

Definition at line 874 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return cSMDcoffs;} ;

int const* const* Molecule::get_cSMDindex ( )

inline

Definition at line 870 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return cSMDindex;} ;

int const* Molecule::get_cSMDindxLen ( )

inline

Definition at line 869 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return cSMDindxLen;} ;

const Real* Molecule::get_cSMDKs ( )

inline

Definition at line 872 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return cSMDKs;} ;

int Molecule::get_cSMDnumInst ( )

inline

Definition at line 868 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return cSMDnumInst;} ;

int const* const* Molecule::get_cSMDpairs ( )

inline

Definition at line 871 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return cSMDpairs;} ;

const Real* Molecule::get_cSMDVels ( )

inline

Definition at line 873 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return cSMDVels;} ;

Dihedral* Molecule::get_dihedral ( int  dnum ) const

inline

Definition at line 1074 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 1149 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

      { return dihedralsByAtom[anum]; }

Bond* Molecule::get_donor ( int  dnum ) const

inline

Definition at line 1102 of file Molecule.h.

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

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

inline

Definition at line 1173 of file Molecule.h.

Referenced by dumpbench().

                                                               {      
      return &all_exclusions[anum];             
  }

Exclusion* Molecule::get_exclusion ( int  ex ) const

inline

Definition at line 1112 of file Molecule.h.

References exclusions.

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

int32* Molecule::get_exclusions_for_atom ( int  anum )

inline

Definition at line 1155 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 1385 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 1345 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 1157 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 1601 of file Molecule.h.

References go_val::cutoff, go_array, and MAX_GO_CHAINS.

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

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

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

inline

Definition at line 1610 of file Molecule.h.

References go_val::epsilon, go_array, and MAX_GO_CHAINS.

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

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

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

inline

Definition at line 1604 of file Molecule.h.

References go_val::epsilonRep, go_array, and MAX_GO_CHAINS.

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

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

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

inline

Definition at line 1613 of file Molecule.h.

References go_val::exp_a, go_array, and MAX_GO_CHAINS.

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

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

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

inline

Definition at line 1616 of file Molecule.h.

References go_val::exp_b, go_array, and MAX_GO_CHAINS.

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

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

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

inline

Definition at line 1619 of file Molecule.h.

References go_val::exp_rep, go_array, and MAX_GO_CHAINS.

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

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

References go_array, MAX_GO_CHAINS, and go_val::sigmaRep.

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

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

GridforceGrid* Molecule::get_gridfrc_grid ( int  gridnum ) const

inline

Definition at line 1276 of file Molecule.h.

References numGridforceGrids.

Referenced by ComputeGridForce::doForce(), 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 1270 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 1071 of file Molecule.h.

Referenced by dumpbench().

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

int32* Molecule::get_impropers_for_atom ( int  anum )

inline

Definition at line 1151 of file Molecule.h.

Referenced by dumpbench().

      { return impropersByAtom[anum]; }  

Lphost* Molecule::get_lphost ( int  atomid ) const

inline

Definition at line 1081 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 1159 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 1321 of file Molecule.h.

Referenced by Sequencer::addMovDragToPosition().

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

Bool Molecule::get_noPC ( )

inline

Definition at line 856 of file Molecule.h.

Referenced by ComputeQM::initialize().

{ return qmNoPC; } ;

int Molecule::get_numQMAtoms ( )

inline

Definition at line 825 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 822 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 823 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 819 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 820 of file Molecule.h.

{return qmAtomGroup[indx]; } ;

Bool Molecule::get_qmcSMD ( )

inline

Definition at line 867 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmcSMD;} ;

int* Molecule::get_qmCustomPCIdxs ( )

inline

Definition at line 865 of file Molecule.h.

Referenced by ComputeQM::initialize().

{ return qmCustomPCIdxs; } ;

int* Molecule::get_qmCustPCSizes ( )

inline

Definition at line 864 of file Molecule.h.

Referenced by ComputeQM::initialize().

{ return qmCustPCSizes; } ;

const BigReal* Molecule::get_qmDummyBondVal ( )

inline

Definition at line 834 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmDummyBondVal; } ;

const char* const* Molecule::get_qmDummyElement ( )

inline

Definition at line 839 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmDummyElement; } ;

const char* const* Molecule::get_qmElements ( )

inline

Definition at line 826 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{return qmElementArray;} ;

const int* const* Molecule::get_qmGrpBonds ( )

inline

Definition at line 837 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmGrpBonds; } ;

Real* Molecule::get_qmGrpChrg ( )

inline

Definition at line 830 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return qmGrpChrg; } ;

Real* Molecule::get_qmGrpID ( )

inline

Definition at line 829 of file Molecule.h.

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

{ return qmGrpID; } ;

Real* Molecule::get_qmGrpMult ( )

inline

Definition at line 831 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return qmGrpMult; } ;

const int* Molecule::get_qmGrpNumBonds ( )

inline

Definition at line 835 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmGrpNumBonds; } ;

const int* Molecule::get_qmGrpSizes ( )

inline

Definition at line 828 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{return qmGrpSizes; } ;

int Molecule::get_qmLSSFreq ( )

inline

Definition at line 847 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSFreq; } ;

int* Molecule::get_qmLSSIdxs ( )

inline

Definition at line 845 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSIdxs;} ;

Mass* Molecule::get_qmLSSMass ( )

inline

Definition at line 846 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSMass; } ;

int* Molecule::get_qmLSSRefIDs ( )

inline

Definition at line 849 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSRefIDs ; } ;

Mass* Molecule::get_qmLSSRefMass ( )

inline

Definition at line 851 of file Molecule.h.

{return qmLSSRefMass; } ;

int* Molecule::get_qmLSSRefSize ( )

inline

Definition at line 850 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSRefSize ; } ;

int Molecule::get_qmLSSResSize ( )

inline

Definition at line 848 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSResidueSize; } ;

int* Molecule::get_qmLSSSize ( )

inline

Definition at line 844 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmLSSSize;} ;

int* Molecule::get_qmMeMMindx ( )

inline

Definition at line 858 of file Molecule.h.

Referenced by ComputeQM::initialize().

{ return qmMeMMindx; } ;

int Molecule::get_qmMeNumBonds ( )

inline

Definition at line 857 of file Molecule.h.

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

{ return qmMeNumBonds; } ;

Real* Molecule::get_qmMeQMGrp ( )

inline

Definition at line 859 of file Molecule.h.

Referenced by ComputeQM::initialize().

{ return qmMeQMGrp; } ;

const int* const* Molecule::get_qmMMBond ( )

inline

Definition at line 836 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmMMBond; } ;

const int* const* Molecule::get_qmMMBondedIndx ( )

inline

Definition at line 838 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmMMBondedIndx; } ;

const int* const* Molecule::get_qmMMChargeTarget ( )

inline

Definition at line 841 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmMMChargeTarget; } ;

const int* Molecule::get_qmMMNumTargs ( )

inline

Definition at line 842 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

{ return qmMMNumTargs; } ;

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

inline

Definition at line 852 of file Molecule.h.

{ return qmClassicSolv;} ;

int Molecule::get_qmNumBonds ( )

inline

Definition at line 833 of file Molecule.h.

{ return qmNumBonds; } ;

int Molecule::get_qmNumGrps ( ) const

inline

Definition at line 827 of file Molecule.h.

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

{return qmNumGrps; } ;

int Molecule::get_qmPCFreq ( )

inline

Definition at line 861 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{ return qmPCFreq; } ;

Bool Molecule::get_qmReplaceAll ( )

inline

Definition at line 854 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

{return qmReplaceAll; } ;

int Molecule::get_qmTotCustPCs ( )

inline

Definition at line 863 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 1327 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 1351 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 1302 of file Molecule.h.

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

Real Molecule::get_stir_startTheta ( int  atomnum ) const

inline

Definition at line 1314 of file Molecule.h.

Referenced by ComputeStir::doForce().

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

Bond* Molecule::getAllAcceptors ( ) const

inline

Definition at line 1108 of file Molecule.h.

{return acceptors;}

Angle* Molecule::getAllAngles ( ) const

inline

Definition at line 1091 of file Molecule.h.

Referenced by buildAngleData().

{return angles;}

Bond* Molecule::getAllBonds ( ) const

inline

Definition at line 1090 of file Molecule.h.

Referenced by buildBondData().

{return bonds;}

Crossterm* Molecule::getAllCrossterms ( ) const

inline

Definition at line 1094 of file Molecule.h.

Referenced by buildCrosstermData().

{return crossterms;}

Dihedral* Molecule::getAllDihedrals ( ) const

inline

Definition at line 1093 of file Molecule.h.

Referenced by buildDihedralData().

{return dihedrals;}

Bond* Molecule::getAllDonors ( ) const

inline

Definition at line 1107 of file Molecule.h.

{return donors;}

Improper* Molecule::getAllImpropers ( ) const

inline

Definition at line 1092 of file Molecule.h.

Referenced by buildImproperData().

{return impropers;}

Lphost* Molecule::getAllLphosts ( ) const

inline

Definition at line 1098 of file Molecule.h.

{ return lphosts; }

BigReal Molecule::GetAtomAlpha ( int  i ) const

inline

Definition at line 482 of file Molecule.h.

References drude_constants::alpha.

                                    {
    return drudeConsts[i].alpha;
  }

AtomNameInfo* Molecule::getAtomNames ( ) const

inline

Definition at line 491 of file Molecule.h.

Referenced by buildAtomData().

{ return atomNames; }

Atom* Molecule::getAtoms ( ) const

inline

Definition at line 490 of file Molecule.h.

References atoms.

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

{ return atoms; }

AtomSegResInfo* Molecule::getAtomSegResInfo ( ) const

inline

Definition at line 494 of file Molecule.h.

Referenced by buildAtomData().

{ return atomSegResids; }

const float* Molecule::getBFactorData ( )

inline

Definition at line 1032 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 478 of file Molecule.h.

Referenced by HomePatch::setLcpoType().

                                 {
    return lcpoParamType;
  }

const float* Molecule::getOccupancyData ( )

inline

Definition at line 1028 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 1246 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 1197 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 1447 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 1407 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 1181 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 1214 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 1230 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 1428 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 1443 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 1296 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 523 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 497 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 503 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 510 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.");
            }