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

Public Attributes

int ss_num_vdw_params
 
int * ss_vdw_type
 
int * ss_index
 
int is_drude_psf
 
int is_lonepairs_psf
 
Real r_om
 
Real r_ohc
 
BigReal tail_corr_ener
 
BigReal tail_corr_virial
 
BigReal tail_corr_dUdl_1
 
BigReal tail_corr_dUdl_2
 
BigReal tail_corr_virial_1
 
BigReal tail_corr_virial_2
 
int numAtoms
 
int numRealBonds
 
int numBonds
 
int numAngles
 
int numDihedrals
 
int suspiciousAlchBonds
 
int alchDroppedAngles
 
int alchDroppedDihedrals
 
int alchDroppedImpropers
 
int numImpropers
 
int numCrossterms
 
int numDonors
 
int numAcceptors
 
int numExclusions
 
int64 numTotalExclusions
 
int num_alch_unpert_Bonds
 
int num_alch_unpert_Angles
 
int num_alch_unpert_Dihedrals
 
Bondalch_unpert_bonds
 
Anglealch_unpert_angles
 
Dihedralalch_unpert_dihedrals
 
int numLonepairs
 Number of lone pairs. 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
 
int32consForceIndexes
 
VectorconsForce
 
int32consTorqueIndexes
 
ConsTorqueParams * consTorqueParams
 
int numCalcBonds
 
int numCalcAngles
 
int numCalcDihedrals
 
int numCalcImpropers
 
int numCalcCrossterms
 
int64 numCalcExclusions
 
int64 numCalcFullExclusions
 
int numCalcTholes
 
int numCalcAnisos
 
int numMultipleDihedrals
 
int numMultipleImpropers
 
HydrogenGroup hydrogenGroup
 
int numGoAtoms
 
int32atomChainTypes
 
int32goSigmaIndices
 
int32goResidIndices
 
RealgoSigmas
 
bool * goWithinCutoff
 
RealgoCoordinates
 
int * goResids
 
PDBgoPDB
 
int goNumLJPair
 
int * goIndxLJA
 
int * goIndxLJB
 
double * goSigmaPairA
 
double * goSigmaPairB
 
int * pointerToGoBeg
 
int * pointerToGoEnd
 
int numPair
 
int numLJPair
 
int numCalcLJPair
 
int * pointerToLJBeg
 
int * pointerToLJEnd
 
int * indxLJA
 
int * indxLJB
 
RealpairC6
 
RealpairC12
 
int * gromacsPair_type
 
int * pointerToGaussBeg
 
int * pointerToGaussEnd
 
int numGaussPair
 
int * indxGaussA
 
int * indxGaussB
 
RealgA
 
RealgMu1
 
RealgiSigma1
 
RealgMu2
 
RealgiSigma2
 
RealgRepulsive
 
BigReal energyNative
 
BigReal energyNonnative
 
int isOccupancyValid
 
int isBFactorValid
 
GoValue go_array [MAX_GO_CHAINS *MAX_GO_CHAINS]
 
int go_indices [MAX_GO_CHAINS+1]
 
int NumGoChains
 

Friends

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

Detailed Description

Definition at line 150 of file Molecule.h.

Constructor & Destructor Documentation

Molecule::Molecule ( SimParameters simParams,
Parameters param 
)

Definition at line 429 of file Molecule.C.

430 {
431  initialize(simParams,param);
432 }
void initialize()
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.

443 {
444  initialize(simParams,param);
445 
446 #ifdef MEM_OPT_VERSION
447  if(simParams->useCompressedPsf)
448  read_mol_signatures(filename, param, cfgList);
449 #else
450  read_psf_file(filename, param);
451  //LCPO
452  if (simParams->LCPOOn)
453  assignLCPOTypes( 0 );
454 #endif
455  }
void initialize()
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().

465 {
466 #ifdef MEM_OPT_VERSION
467  NAMD_die("Sorry, plugin IO is not supported in the memory optimized version.");
468 #else
469  initialize(simParams, param);
470  numAtoms = natoms;
471  int optflags = MOLFILE_BADOPTIONS;
472  molfile_atom_t *atomarray = (molfile_atom_t *) malloc(natoms*sizeof(molfile_atom_t));
473  memset(atomarray, 0, natoms*sizeof(molfile_atom_t));
474 
475  //1a. read basic atoms information
476  int rc = pIOHdl->read_structure(pIOFileHdl, &optflags, atomarray);
477  if (rc != MOLFILE_SUCCESS && rc != MOLFILE_NOSTRUCTUREDATA) {
478  free(atomarray);
479  NAMD_die("ERROR: plugin failed reading structure data");
480  }
481  if(optflags == MOLFILE_BADOPTIONS) {
482  free(atomarray);
483  NAMD_die("ERROR: plugin didn't initialize optional data flags");
484  }
485  if(optflags & MOLFILE_OCCUPANCY) {
486  setOccupancyData(atomarray);
487  }
488  if(optflags & MOLFILE_BFACTOR) {
489  setBFactorData(atomarray);
490  }
491  //1b. load basic atoms information to the molecule object
492  plgLoadAtomBasics(atomarray);
493  free(atomarray);
494 
495  //2a. read bonds
496  //indices are one-based in read_bonds
497  int *from, *to;
498  float *bondorder;
499  int *bondtype, nbondtypes;
500  char **bondtypename;
501  if(pIOHdl->read_bonds!=NULL) {
502  if(pIOHdl->read_bonds(pIOFileHdl, &numBonds, &from, &to, &bondorder,
503  &bondtype, &nbondtypes, &bondtypename)){
504  NAMD_die("ERROR: failed reading bond information.");
505  }
506  }
507  //2b. load bonds information to the molecule object
508  if(numBonds!=0) {
509  plgLoadBonds(from,to);
510  }
511 
512  //3a. read other bonded structures
513  int *plgAngles, *plgDihedrals, *plgImpropers, *plgCterms;
514  int ctermcols, ctermrows;
515  int *angletypes, numangletypes, *dihedraltypes, numdihedraltypes;
516  int *impropertypes, numimpropertypes;
517  char **angletypenames, **dihedraltypenames, **impropertypenames;
518 
519  plgAngles=plgDihedrals=plgImpropers=plgCterms=NULL;
520  if(pIOHdl->read_angles!=NULL) {
521  if(pIOHdl->read_angles(pIOFileHdl,
522  &numAngles, &plgAngles,
523  &angletypes, &numangletypes, &angletypenames,
524  &numDihedrals, &plgDihedrals,
525  &dihedraltypes, &numdihedraltypes, &dihedraltypenames,
526  &numImpropers, &plgImpropers,
527  &impropertypes, &numimpropertypes, &impropertypenames,
528  &numCrossterms, &plgCterms, &ctermcols, &ctermrows)) {
529  NAMD_die("ERROR: failed reading angle information.");
530  }
531  }
532  //3b. load other bonded structures to the molecule object
533  if(numAngles!=0) plgLoadAngles(plgAngles);
534  if(numDihedrals!=0) plgLoadDihedrals(plgDihedrals);
535  if(numImpropers!=0) plgLoadImpropers(plgImpropers);
536  if(numCrossterms!=0) plgLoadCrossterms(plgCterms);
537 
539  build_atom_status();
540  //LCPO
541  if (simParams->LCPOOn)
542  assignLCPOTypes( 2 );
543 #endif
544 }
int numBonds
Definition: Molecule.h:559
void setOccupancyData(molfile_atom_t *atomarray)
Definition: Molecule.C:3236
void setBFactorData(molfile_atom_t *atomarray)
Definition: Molecule.C:3243
int numRealBonds
Definition: Molecule.h:558
int numAngles
Definition: Molecule.h:560
int numAtoms
Definition: Molecule.h:556
int numCrossterms
Definition: Molecule.h:567
void NAMD_die(const char *err_msg)
Definition: common.C:83
int numDihedrals
Definition: Molecule.h:561
int numImpropers
Definition: Molecule.h:566
void initialize()
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.

559 {
560  /* Check to see if each array was ever allocated. If it was */
561  /* then free it */
562  if (atoms != NULL)
563  delete [] atoms;
564 
565  if (atomNames != NULL)
566  {
567  // subarrarys allocated from arena - automatically deleted
568  delete [] atomNames;
569  }
570  delete nameArena;
571 
572  if (resLookup != NULL)
573  delete resLookup;
574 
575  // DRUDE: free arrays read from PSF
576  if (drudeConsts != NULL) delete [] drudeConsts;
577  if (lphosts != NULL) delete [] lphosts;
578  if (anisos != NULL) delete [] anisos;
579  if (tholes != NULL) delete [] tholes;
580  if (lphostIndexes != NULL) delete [] lphostIndexes;
581  // DRUDE
582 
583  //LCPO
584  if (lcpoParamType != NULL) delete [] lcpoParamType;
585 
586  #ifdef MEM_OPT_VERSION
587  if(eachAtomSig) delete [] eachAtomSig;
588  if(atomSigPool) delete [] atomSigPool;
589  #else
590  if (bonds != NULL)
591  delete [] bonds;
592 
593  if (angles != NULL)
594  delete [] angles;
595 
596  if (dihedrals != NULL)
597  delete [] dihedrals;
598 
599  if (impropers != NULL)
600  delete [] impropers;
601 
602  if (crossterms != NULL)
603  delete [] crossterms;
604 
605  if (exclusions != NULL)
606  delete [] exclusions;
607  #endif
608 
609  if (donors != NULL)
610  delete [] donors;
611 
612  if (acceptors != NULL)
613  delete [] acceptors;
614 
615  #ifdef MEM_OPT_VERSION
616  if(exclSigPool) delete [] exclSigPool;
617  if(exclChkSigPool) delete [] exclChkSigPool;
618  if(eachAtomExclSig) delete [] eachAtomExclSig;
619  if(fixedAtomsSet) delete fixedAtomsSet;
620  if(constrainedAtomsSet) delete constrainedAtomsSet;
621  #else
622  if (bondsByAtom != NULL)
623  delete [] bondsByAtom;
624 
625  if (anglesByAtom != NULL)
626  delete [] anglesByAtom;
627 
628  if (dihedralsByAtom != NULL)
629  delete [] dihedralsByAtom;
630 
631  if (impropersByAtom != NULL)
632  delete [] impropersByAtom;
633 
634  if (crosstermsByAtom != NULL)
635  delete [] crosstermsByAtom;
636 
637  if (exclusionsByAtom != NULL)
638  delete [] exclusionsByAtom;
639 
640  if (fullExclusionsByAtom != NULL)
641  delete [] fullExclusionsByAtom;
642 
643  if (modExclusionsByAtom != NULL)
644  delete [] modExclusionsByAtom;
645 
646  if (all_exclusions != NULL)
647  delete [] all_exclusions;
648 
649  // JLai
650  if (gromacsPairByAtom != NULL)
651  delete [] gromacsPairByAtom;
652  // End of JLai
653 
654  // DRUDE
655  if (tholesByAtom != NULL)
656  delete [] tholesByAtom;
657  if (anisosByAtom != NULL)
658  delete [] anisosByAtom;
659  // DRUDE
660  #endif
661 
662  //LCPO
663  if (lcpoParamType != NULL)
664  delete [] lcpoParamType;
665 
666  if (fixedAtomFlags != NULL)
667  delete [] fixedAtomFlags;
668 
669  if (stirIndexes != NULL)
670  delete [] stirIndexes;
671 
672 
673  #ifdef MEM_OPT_VERSION
674  if(clusterSigs != NULL){
675  delete [] clusterSigs;
676  }
677  #else
678  if (cluster != NULL)
679  delete [] cluster;
680  #endif
681  if (clusterSize != NULL)
682  delete [] clusterSize;
683 
684  if (exPressureAtomFlags != NULL)
685  delete [] exPressureAtomFlags;
686 
687  if (rigidBondLengths != NULL)
688  delete [] rigidBondLengths;
689 
690 //fepb
691  if (fepAtomFlags != NULL)
692  delete [] fepAtomFlags;
693  if (alch_unpert_bonds != NULL)
694  delete [] alch_unpert_bonds;
695  if (alch_unpert_angles != NULL)
696  delete [] alch_unpert_angles;
697  if (alch_unpert_dihedrals != NULL)
698  delete [] alch_unpert_dihedrals;
699 //fepe
700 
701 //soluteScaling
702  if (ssAtomFlags != NULL)
703  delete [] ssAtomFlags;
704  if (ss_vdw_type != NULL)
705  delete [] ss_vdw_type;
706  if (ss_index != NULL)
707  delete [] ss_index;
708 //soluteScaling
709 
710  if (qmAtomGroup != NULL)
711  delete [] qmAtomGroup;
712 
713  if (qmAtmIndx != NULL)
714  delete [] qmAtmIndx;
715 
716  if (qmAtmChrg != NULL)
717  delete [] qmAtmChrg;
718 
719 
720  if (qmGrpNumBonds != NULL)
721  delete [] qmGrpNumBonds;
722 
723  if (qmGrpSizes != NULL)
724  delete [] qmGrpSizes;
725 
726  if (qmDummyBondVal != NULL)
727  delete [] qmDummyBondVal;
728 
729  if (qmMMNumTargs != NULL)
730  delete [] qmMMNumTargs;
731 
732  if (qmGrpID != NULL)
733  delete [] qmGrpID;
734 
735  if (qmGrpChrg != NULL)
736  delete [] qmGrpChrg;
737 
738  if (qmGrpMult != NULL)
739  delete [] qmGrpMult;
740 
741  if (qmMeMMindx != NULL)
742  delete [] qmMeMMindx;
743 
744  if (qmMeQMGrp != NULL)
745  delete [] qmMeQMGrp;
746 
747  if (qmLSSSize != NULL)
748  delete [] qmLSSSize;
749 
750  if (qmLSSIdxs != NULL)
751  delete [] qmLSSIdxs;
752 
753  if (qmLSSMass != NULL)
754  delete [] qmLSSMass;
755 
756  if (qmLSSRefSize != NULL)
757  delete [] qmLSSRefSize;
758 
759  if (qmLSSRefIDs != NULL)
760  delete [] qmLSSRefIDs;
761 
762  if (qmLSSRefMass != NULL)
763  delete [] qmLSSRefMass;
764 
765  if (qmMMBond != NULL) {
766  for(int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
767  if (qmMMBond[grpIndx] != NULL)
768  delete [] qmMMBond[grpIndx];
769  }
770  delete [] qmMMBond;
771  }
772 
773  if (qmGrpBonds != NULL) {
774  for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
775  if (qmGrpBonds[grpIndx] != NULL)
776  delete [] qmGrpBonds[grpIndx];
777  }
778  delete [] qmGrpBonds;
779  }
780 
781  if (qmMMBondedIndx != NULL) {
782  for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
783  if (qmMMBondedIndx[grpIndx] != NULL)
784  delete [] qmMMBondedIndx[grpIndx];
785  }
786  delete [] qmMMBondedIndx;
787  }
788 
789  if (qmMMChargeTarget != NULL) {
790  for(int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
791  if (qmMMChargeTarget[grpIndx] != NULL)
792  delete [] qmMMChargeTarget[grpIndx];
793  }
794  delete [] qmMMChargeTarget;
795  }
796 
797  if (qmElementArray != NULL){
798  for(int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
799  if (qmElementArray[atmInd] != NULL)
800  delete [] qmElementArray[atmInd];
801  }
802  delete [] qmElementArray;
803  }
804 
805  if (qmDummyElement != NULL){
806  for(int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
807  if (qmDummyElement[atmInd] != NULL)
808  delete [] qmDummyElement[atmInd];
809  }
810  delete [] qmDummyElement;
811  }
812 
813  if (qmCustPCSizes != NULL){
814  delete [] qmCustPCSizes;
815  }
816 
817  if (qmCustomPCIdxs != NULL){
818  delete [] qmCustomPCIdxs;
819  }
820 
821  if (cSMDindex != NULL) {
822  for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
823  if (cSMDindex[grpIndx] != NULL)
824  delete [] cSMDindex[grpIndx];
825  }
826  delete [] cSMDindex;
827  }
828 
829  if (cSMDpairs != NULL) {
830  for(int instIndx = 0 ; instIndx < cSMDnumInst; instIndx++) {
831  if (cSMDpairs[instIndx] != NULL)
832  delete [] cSMDpairs[instIndx];
833  }
834  delete [] cSMDpairs;
835  }
836 
837  if (cSMDindxLen != NULL)
838  delete [] cSMDindxLen;
839  if (cSMDKs != NULL)
840  delete [] cSMDKs;
841  if (cSMDVels != NULL)
842  delete [] cSMDVels;
843  if (cSMDcoffs != NULL)
844  delete [] cSMDcoffs;
845 
846  #ifndef MEM_OPT_VERSION
847  delete arena;
848  delete exclArena;
849  #endif
850 }
Angle * alch_unpert_angles
Definition: Molecule.h:576
Dihedral * alch_unpert_dihedrals
Definition: Molecule.h:577
Bond * alch_unpert_bonds
Definition: Molecule.h:575
int numAtoms
Definition: Molecule.h:556
int * ss_index
Definition: Molecule.h:458
int * ss_vdw_type
Definition: Molecule.h:457
HashPool< AtomSigInfo > atomSigPool
Definition: CompressPsf.C:313

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

1049  {
1050  #ifdef MEM_OPT_VERSION
1051  return atomChargePool[eachAtomCharge[anum]];
1052  #else
1053  return(atoms[anum].charge);
1054  #endif
1055  }
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
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().

1039  {
1040  #ifdef MEM_OPT_VERSION
1041  return atomMassPool[eachAtomMass[anum]];
1042  #else
1043  return(atoms[anum].mass);
1044  #endif
1045  }
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().

1539 {
1540  int bondNum; // Bonds in bonded list
1541  int angleNum; // Angles in angle list
1542  int dihedralNum; // Dihedrals in dihedral list
1543  int *bonds;
1544  int *angles;
1545  int *dihedrals;
1546  Bond *bond; // Temporary bond structure
1547  Angle *angle; // Temporary angle structure
1548  Dihedral *dihedral; // Temporary dihedral structure
1549 
1550  DebugM(2,"atoms_1to4(" << atom1 << "," << atom2 << ")" << std::endl);
1551 
1552  bonds = get_bonds_for_atom(atom1);
1553  bondNum = *bonds;
1554  while(bondNum != -1) {
1555  bond = get_bond(bondNum);
1556  DebugM(2,"bond atom1:" << bond->atom1 << ", atom2:" << bond->atom2 << std::endl);
1557  if (atom2 == bond->atom1 || atom2 == bond->atom2) {
1558  return TRUE;
1559  }
1560  bondNum = *(++bonds);
1561  }
1562 
1563  bonds = get_bonds_for_atom(atom2);
1564  bondNum = *bonds;
1565  while(bondNum != -1) {
1566  bond = get_bond(bondNum);
1567  DebugM(2,"bond atom1:" << bond->atom1 << ", atom2:" << bond->atom2 << std::endl);
1568  if (atom1 == bond->atom1 || atom1 == bond->atom2) {
1569  return TRUE;
1570  }
1571  bondNum = *(++bonds);
1572  }
1573 
1574  angles = get_angles_for_atom(atom1);
1575  angleNum = *angles;
1576  while(angleNum != -1) {
1577  angle = get_angle(angleNum);
1578  DebugM(2,"angle atom1:" << angle->atom1 << ", atom2:" << angle->atom2 << ", atom3:" << angle->atom3 << std::endl);
1579  if (atom2 == angle->atom1 || atom2 == angle->atom2 || atom2 == angle->atom3) {
1580  return TRUE;
1581  }
1582  angleNum = *(++angles);
1583  }
1584 
1585  angles = get_angles_for_atom(atom2);
1586  angleNum = *angles;
1587  while(angleNum != -1) {
1588  angle = get_angle(angleNum);
1589  DebugM(2,"angle atom1:" << angle->atom1 << ", atom2:" << angle->atom2 << ", atom3:" << angle->atom3 << std::endl);
1590  if (atom1 == angle->atom1 || atom1 == angle->atom2 || atom1 == angle->atom3) {
1591  return TRUE;
1592  }
1593  angleNum = *(++angles);
1594  }
1595 
1596  dihedrals = get_dihedrals_for_atom(atom1);
1597  dihedralNum = *dihedrals;
1598  while(dihedralNum != -1) {
1599  dihedral = get_dihedral(dihedralNum);
1600  DebugM(2,"dihedral atom1:" << dihedral->atom1 << ", atom2:" << dihedral->atom2 << ", atom3:" << dihedral->atom3 << ", atom4:" << dihedral->atom4 << std::endl);
1601  if (atom2 == dihedral->atom1 || atom2 == dihedral->atom2 \
1602  || atom2 == dihedral->atom3 || atom2 == dihedral->atom4) {
1603  return TRUE;
1604  }
1605  dihedralNum = *(++dihedrals);
1606  }
1607 
1608  dihedrals = get_dihedrals_for_atom(atom2);
1609  dihedralNum = *dihedrals;
1610  while(dihedralNum != -1) {
1611  dihedral = get_dihedral(dihedralNum);
1612  DebugM(2,"dihedral atom1:" << dihedral->atom1 << ", atom2:" << dihedral->atom2 << ", atom3:" << dihedral->atom3 << ", atom4:" << dihedral->atom4 << std::endl);
1613  if (atom1 == dihedral->atom1 || atom1 == dihedral->atom2 \
1614  || atom1 == dihedral->atom3 || atom1 == dihedral->atom4) {
1615  return TRUE;
1616  }
1617  dihedralNum = *(++dihedrals);
1618  }
1619 
1620  return FALSE;
1621 }
int32 * get_dihedrals_for_atom(int anum)
Definition: Molecule.h:1149
int32 atom3
Definition: structures.h:65
#define DebugM(x, y)
Definition: Debug.h:59
#define FALSE
Definition: common.h:116
int32 * get_angles_for_atom(int anum)
Definition: Molecule.h:1147
int32 atom4
Definition: structures.h:66
int32 atom1
Definition: structures.h:48
Dihedral * get_dihedral(int dnum) const
Definition: Molecule.h:1074
int32 atom2
Definition: structures.h:56
int32 atom1
Definition: structures.h:55
int32 * get_bonds_for_atom(int anum)
Definition: Molecule.h:1145
int32 atom3
Definition: structures.h:57
Bond * get_bond(int bnum) const
Definition: Molecule.h:1065
int32 atom2
Definition: structures.h:64
Angle * get_angle(int anum) const
Definition: Molecule.h:1068
#define TRUE
Definition: common.h:117
int32 atom1
Definition: structures.h:63
int32 atom2
Definition: structures.h:49
Index Molecule::atomvdwtype ( int  anum) const
inline

Definition at line 1058 of file Molecule.h.

References atoms.

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

1059  {
1060  return(atoms[anum].vdw_type);
1061  }
void Molecule::build_alch_unpert_bond_lists ( char *  )
void Molecule::build_constant_forces ( char *  )
void Molecule::build_constorque_params ( StringList ,
StringList ,
StringList ,
StringList ,
StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_constraint_params ( StringList ,
StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_exPressure_atoms ( StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_extra_bonds ( Parameters parameters,
StringList file 
)
void Molecule::build_fep_flags ( StringList ,
StringList ,
PDB ,
char *  ,
const char *   
)
void Molecule::build_fixed_atoms ( StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_go_arrays ( StringList goCoordFile,
char *  cwd 
)

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

Definition at line 950 of file GoMolecule.C.

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

Referenced by NamdState::loadStructure().

952 {
953  DebugM(3,"->build_go_arrays" << std::endl);
954  //PDB *goPDB; // Pointer to PDB object to use
955  int bcol = 4; // Column that data is in
956  int32 chainType = 0; // b value from PDB file
957  int i; // Loop counter
958  int j; // Loop counter
959  BigReal atomAtomDist; // Distance between two atoms -- JLai put back
960  int resid1; // Residue ID for first atom
961  int resid2; // Residue ID for second atom
962  int residDiff; // Difference between resid1 and resid2
963  int goIndex; // Index into the goCoordinates array
964  int goIndx; // Index into the goCoordinates array
965  char filename[129]; // Filename
966 
967  //JLai
968  BigReal nativeEnergy, *native;
969  BigReal nonnativeEnergy, *nonnative;
970  nativeEnergy = 0;
971  nonnativeEnergy = 0;
972  native = &nativeEnergy;
973  nonnative = &nonnativeEnergy;
974 
975  long nativeContacts = 0;
976  long nonnativeContacts = 0;
977 
978  // Get the PDB object that contains the Go coordinates. If
979  // the user gave another file name, use it. Otherwise, just use
980  // the PDB file that has the initial coordinates.
981  if (goCoordFile == NULL)
982  {
983  //goPDB = initial_pdb;
984  NAMD_die("Error: goCoordFile is NULL - build_go_arrays");
985  }
986  else
987  {
988  if (goCoordFile->next != NULL)
989  {
990  NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
991  }
992 
993  if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
994  {
995  strcpy(filename, goCoordFile->data);
996  }
997  else
998  {
999  strcpy(filename, cwd);
1000  strcat(filename, goCoordFile->data);
1001  }
1002 
1003  goPDB = new PDB(filename);
1004  if ( goPDB == NULL )
1005  {
1006  NAMD_die("goPDB memory allocation failed in Molecule::build_go_arrays");
1007  }
1008 
1009  if (goPDB->num_atoms() != numAtoms)
1010  {
1011  NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
1012  }
1013  }
1014 
1015  // Allocate the array to hold Go atom indices into the sigma array
1016  goSigmaIndices = new int32[numAtoms];
1017 
1018  if (goSigmaIndices == NULL) {
1019  NAMD_die("goSigmaIndices memory allocation failed in Molecule::build_go_arrays");
1020  }
1021 
1022  numGoAtoms = 0;
1023 
1024  // Loop through all the atoms and get the Go chain types
1025  for (i=0; i<numAtoms; i++) {
1026  chainType = (int32)(goPDB->atom(i))->occupancy();
1027  if ( chainType != 0 ) {
1028  DebugM(3,"build_go_arrays - atom:" << i << std::endl);
1030  numGoAtoms++;
1031  }
1032  else {
1033  goSigmaIndices[i] = -1;
1034  }
1035  }
1036 
1037  // Allocate the array to hold the sigma values for each Go atom pair
1049  // Allocate the array to hold the chain types
1051 
1052  if (atomChainTypes == NULL) {
1053  NAMD_die("atomChainTypes memory allocation failed in Molecule::build_go_arrays");
1054  }
1055 
1056  // Allocate the array to hold (x,y,z) coordinates for all of the Go atoms
1057  goCoordinates = new Real[numGoAtoms*3];
1058 
1059  if (goCoordinates == NULL) {
1060  NAMD_die("goCoordinates memory allocation failed in Molecule::build_go_arrays");
1061  }
1062 
1063  goResids = new int[numGoAtoms];
1064 
1065  // Allocate the array to hold PDB residu IDs for all of the Go atoms
1066  if (goResids == NULL) {
1067  NAMD_die("goResids memory allocation failed in Molecule::build_go_arrays");
1068  }
1069 
1070  for (i=0; i<numAtoms; i++) {
1071  goIndex = goSigmaIndices[i];
1072  if (goIndex != -1) {
1073  // Assign the chainType value!
1074  // Get the chainType from the occupancy field
1075  atomChainTypes[goIndex] = (int32)(goPDB->atom(i))->occupancy();
1076  goCoordinates[goIndex*3] = goPDB->atom(i)->xcoor();
1077  goCoordinates[goIndex*3 + 1] = goPDB->atom(i)->ycoor();
1078  goCoordinates[goIndex*3 + 2] = goPDB->atom(i)->zcoor();
1079  goResids[goIndex] = goPDB->atom(i)->residueseq();
1080  }
1081  }
1082  // JLai
1083  energyNative = 0;
1084  energyNonnative = 0;
1085  //printf("INIT ENERGY: (N) %f (NN) %f\n", energyNative, energyNonnative);
1086  for (i=0; i<numAtoms-1; i++) {
1087  goIndex = goSigmaIndices[i];
1088  if (goIndex != -1) {
1089  for (j=i+1; j<numAtoms;j++) {
1090  goIndx = goSigmaIndices[j];
1091  if (goIndx != -1) {
1092  resid1 = (goPDB->atom(i))->residueseq();
1093  resid2 = (goPDB->atom(j))->residueseq();
1094  residDiff = resid2 - resid1;
1095  if (residDiff < 0) residDiff = -residDiff;
1096  if (atomChainTypes[goIndex] && atomChainTypes[goIndx] &&
1097  !(this->go_restricted(atomChainTypes[goIndex],atomChainTypes[goIndx],residDiff)) &&
1098  !atoms_1to4(i,j)) {
1099  atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
1100  pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
1101  pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
1102  if (atomAtomDist <= this->get_go_cutoff(atomChainTypes[goIndex],atomChainTypes[goIndx]) ) {
1103  nativeContacts++;
1104  } else {
1105  nonnativeContacts++;
1106  }
1107  }
1108  }
1109  }
1110  }
1111  }
1112  iout << iINFO << "Number of UNIQUE native contacts: " << nativeContacts << "\n" << endi;
1113  iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts << "\n" << endi;
1114 
1115  // If we had to create a PDB object, delete it now
1116  if (goCoordFile != NULL) {
1117  delete goPDB;
1118  }
1119 
1120  return;
1121 }
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
Real get_go_cutoff(int chain1, int chain2)
Definition: Molecule.h:1601
Definition: PDB.h:35
PDB * goPDB
Definition: Molecule.h:650
short int32
Definition: dumpdcd.c:24
Real * goCoordinates
Definition: Molecule.h:648
float Real
Definition: common.h:107
#define DebugM(x, y)
Definition: Debug.h:59
BigReal zcoor(void)
Definition: PDBData.C:433
#define iout
Definition: InfoStream.h:87
int num_atoms(void)
Definition: PDB.C:323
int32 * atomChainTypes
Definition: Molecule.h:643
PDBAtom * atom(int place)
Definition: PDB.C:394
BigReal ycoor(void)
Definition: PDBData.C:429
int residueseq(void)
Definition: PDBData.C:411
int * goResids
Definition: Molecule.h:649
BigReal energyNonnative
Definition: Molecule.h:685
Bool go_restricted(int, int, int)
Definition: GoMolecule.C:525
int numAtoms
Definition: Molecule.h:556
void NAMD_die(const char *err_msg)
Definition: common.C:83
BigReal xcoor(void)
Definition: PDBData.C:425
StringList * next
Definition: ConfigList.h:49
char * data
Definition: ConfigList.h:48
BigReal energyNative
Definition: Molecule.h:684
infostream & endi(infostream &s)
Definition: InfoStream.C:38
Bool atoms_1to4(unsigned int, unsigned int)
Definition: GoMolecule.C:1537
int32 * goSigmaIndices
Definition: Molecule.h:644
int numGoAtoms
Definition: Molecule.h:642
double BigReal
Definition: common.h:112
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().

79  {
80  iout << iINFO << "Building Go Parameters" << "\n" << endi;
81 #ifdef MEM_OPT_VERSION
82  NAMD_die("Go forces are not supported in memory-optimized builds.");
83 #else
84  build_lists_by_atom();
85 #endif
86  int iterator = 0;
87  do
88  {
89  iout << iINFO << "Reading Go File: " << iterator << "\n" << endi;
90  read_go_file(g->data);
91  g = g->next;
92  iterator++;
93  } while ( g != NULL && iterator < 100);
94 }
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
#define iout
Definition: InfoStream.h:87
void read_go_file(char *)
Definition: GoMolecule.C:113
void NAMD_die(const char *err_msg)
Definition: common.C:83
StringList * next
Definition: ConfigList.h:49
char * data
Definition: ConfigList.h:48
infostream & endi(infostream &s)
Definition: InfoStream.C:38
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(), StringList::next, PDB::num_atoms(), numAtoms, and numGoAtoms.

Referenced by NamdState::loadStructure().

579 {
580  DebugM(3,"->build_go_sigmas" << std::endl);
581  PDB *goPDB; // Pointer to PDB object to use
582  int bcol = 4; // Column that data is in
583  int32 chainType = 0; // b value from PDB file
584  int i; // Loop counter
585  int j; // Loop counter
586  int resid1; // Residue ID for first atom
587  int resid2; // Residue ID for second atom
588  int residDiff; // Difference between resid1 and resid2
589  Real sigma; // Sigma calculated for a Go atom pair
590  Real atomAtomDist; // Distance between two atoms
591  Real exp_a; // First exponent in L-J like formula
592  Real exp_b; // Second exponent in L-J like formula
593  char filename[129]; // Filename
594 
595  //JLai
596  BigReal nativeEnergy, *native;
597  BigReal nonnativeEnergy, *nonnative;
598  nativeEnergy = 0;
599  nonnativeEnergy = 0;
600  native = &nativeEnergy;
601  nonnative = &nonnativeEnergy;
602 
603  long nativeContacts = 0;
604  long nonnativeContacts = 0;
605 
606  // Get the PDB object that contains the Go coordinates. If
607  // the user gave another file name, use it. Otherwise, just use
608  // the PDB file that has the initial coordinates.
609  if (goCoordFile == NULL)
610  {
611  //goPDB = initial_pdb;
612  NAMD_die("Error: goCoordFile is NULL - build_go_sigmas");
613  }
614  else
615  {
616  if (goCoordFile->next != NULL)
617  {
618  NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
619  }
620 
621  if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
622  {
623  strcpy(filename, goCoordFile->data);
624  }
625  else
626  {
627  strcpy(filename, cwd);
628  strcat(filename, goCoordFile->data);
629  }
630 
631  goPDB = new PDB(filename);
632  if ( goPDB == NULL )
633  {
634  NAMD_die("Memory allocation failed in Molecule::build_go_sigmas");
635  }
636 
637  if (goPDB->num_atoms() != numAtoms)
638  {
639  NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
640  }
641  }
642  // Allocate the array to hold the chain types
644  // Allocate the array to hold Go atom indices into the sigma array
646 
647  if (atomChainTypes == NULL) {
648  NAMD_die("memory allocation failed in Molecule::build_go_sigmas");
649  }
650 
651  numGoAtoms = 0;
652 
653  // Loop through all the atoms and get the Go chain types
654  for (i=0; i<numAtoms; i++) {
655  // Get the chainType from the occupancy field
656  chainType = (int32)(goPDB->atom(i))->occupancy();
657  // Assign the chainType value
658  if ( chainType != 0 ) {
659  //DebugM(3,"build_go_sigmas - atom:" << i << ", chainType:" << chainType << std::endl);
660  atomChainTypes[i] = chainType;
662  numGoAtoms++;
663  }
664  else {
665  atomChainTypes[i] = 0;
666  goSigmaIndices[i] = -1;
667  }
668  //printf("CT: %d %d %d %d\n",i,numGoAtoms,atomChainTypes[i],goSigmaIndices[i]);
669  }
670 
671  // Allocate the array to hold the sigma values for each Go atom pair
674  for (i=0; i<numGoAtoms; i++) {
675  for (j=0; j<numGoAtoms; j++) {
676  goSigmas[i*numGoAtoms + j] = -1.0;
677  goWithinCutoff[i*numGoAtoms + j] = false;
678  }
679  }
680  // Loop through all atom pairs and calculate sigmas
681  DebugM(3," numAtoms=" << numAtoms << std::endl);
682  for (i=0; i<numAtoms; i++) {
683  //DebugM(3," i=" << i << std::endl);
684  resid1 = (goPDB->atom(i))->residueseq();
685  //DebugM(3," resid1=" << resid1 << std::endl);
686  //if ( goSigmaIndices[i] != -1) {
687  // goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[i]] = 0.0;
688  //}
689  for (j=i+1; j<numAtoms; j++) {
690  //DebugM(3," j=" << j << std::endl);
691  resid2 = (goPDB->atom(j))->residueseq();
692  //printf("GSIi %d %d %d\n",i,numAtoms,goSigmaIndices[i]);
693  //printf("SAN CHECK: %d\n",goSigmaIndices[37]);
694  //printf("GSIj %d %d %d\n",j,numAtoms,goSigmaIndices[j]);
695  //printf("ATOMS_1to4 %d\n",!atoms_1to4(i,j));
696  //DebugM(3," resid2=" << resid2 << std::endl);
697  // if goSigmaIndices aren't defined, don't set anything in goSigmas
698  if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 && !atoms_1to4(i,j) ) {
699  //printf("TAKING DIFFERENCE\n");
700  residDiff = resid2 - resid1;
701  //printf("RESIDDIFF %d\n",residDiff);
702  if (residDiff < 0) residDiff = -residDiff;
703  //printf("RESIDDIFF2 %d\n",residDiff);
704  // if this is a Go atom pair that is not restricted
705  // calculate sigma
706  // sigmas are initially set to -1.0 if the atom pair fails go_restricted
707  //printf("CHECKING LOOPING\n");
708  if ( atomChainTypes[i] && atomChainTypes[j] &&
709  !(this->go_restricted(atomChainTypes[i],atomChainTypes[j],residDiff)) ) {
710  atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
711  pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
712  pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
713  if ( atomAtomDist <= this->get_go_cutoff(atomChainTypes[i],atomChainTypes[j]) ) {
714  exp_a = this->get_go_exp_a(atomChainTypes[i],atomChainTypes[j]);
715  exp_b = this->get_go_exp_b(atomChainTypes[i],atomChainTypes[j]);
716  sigma = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
717  goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = sigma;
718  goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = sigma;
719  goWithinCutoff[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = true;
720  goWithinCutoff[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = true;
721  nativeContacts++;
722  } else {
723  goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = 0.0;
724  goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = 0.0;
725  nonnativeContacts++;
726  }
727  } else {
728  goSigmas[goSigmaIndices[i]*numGoAtoms + goSigmaIndices[j]] = -1.0;
729  goSigmas[goSigmaIndices[j]*numGoAtoms + goSigmaIndices[i]] = -1.0;
730  }
731  }
732  }
733  }
734 
735  iout << iINFO << "Number of UNIQUE native contacts: " << nativeContacts << "\n" << endi;
736  iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts << "\n" << endi;
737 
738  // If we had to create a PDB object, delete it now
739  if (goCoordFile != NULL) {
740  delete goPDB;
741  }
742 
743  return;
744 }
int get_go_exp_a(int chain1, int chain2)
Definition: Molecule.h:1613
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
Real get_go_cutoff(int chain1, int chain2)
Definition: Molecule.h:1601
Definition: PDB.h:35
PDB * goPDB
Definition: Molecule.h:650
short int32
Definition: dumpdcd.c:24
float Real
Definition: common.h:107
#define DebugM(x, y)
Definition: Debug.h:59
bool * goWithinCutoff
Definition: Molecule.h:647
#define iout
Definition: InfoStream.h:87
int num_atoms(void)
Definition: PDB.C:323
int32 * atomChainTypes
Definition: Molecule.h:643
PDBAtom * atom(int place)
Definition: PDB.C:394
Real * goSigmas
Definition: Molecule.h:646
Bool go_restricted(int, int, int)
Definition: GoMolecule.C:525
int numAtoms
Definition: Molecule.h:556
void NAMD_die(const char *err_msg)
Definition: common.C:83
int get_go_exp_b(int chain1, int chain2)
Definition: Molecule.h:1616
StringList * next
Definition: ConfigList.h:49
char * data
Definition: ConfigList.h:48
infostream & endi(infostream &s)
Definition: InfoStream.C:38
Bool atoms_1to4(unsigned int, unsigned int)
Definition: GoMolecule.C:1537
int32 * goSigmaIndices
Definition: Molecule.h:644
int numGoAtoms
Definition: Molecule.h:642
double BigReal
Definition: common.h:112
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(), StringList::next, PDB::num_atoms(), numAtoms, numGoAtoms, numLJPair, pointerToGoBeg, pointerToGoEnd, and sort.

Referenced by NamdState::loadStructure().

749 {
750  DebugM(3,"->build_go_sigmas2" << std::endl);
751  PDB *goPDB; // Pointer to PDB object to use
752  int bcol = 4; // Column that data is in
753  int32 chainType = 0; // b value from PDB file
754  int32 residType = 0; // resid value from PDB file
755  int i; // Loop counter
756  int j; // Loop counter
757  int resid1; // Residue ID for first atom
758  int resid2; // Residue ID for second atom
759  int residDiff; // Difference between resid1 and resid2
760  Real sigma; // Sigma calculated for a Go atom pair
761  Real atomAtomDist; // Distance between two atoms
762  Real exp_a; // First exponent in L-J like formula
763  Real exp_b; // Second exponent in L-J like formula
764  char filename[129]; // Filename
765 
766  //JLai
767  int numLJPair = 0;
768  long nativeContacts = 0;
769  long nonnativeContacts = 0;
770 
771  // Get the PDB object that contains the Go coordinates. If
772  // the user gave another file name, use it. Otherwise, just use
773  // the PDB file that has the initial coordinates.
774  if (goCoordFile == NULL)
775  {
776  //goPDB = initial_pdb;
777  NAMD_die("Error: goCoordFile is NULL - build_go_sigmas2");
778  }
779  else
780  {
781  if (goCoordFile->next != NULL)
782  {
783  NAMD_die("Multiple definitions of Go atoms PDB file in configuration file");
784  }
785 
786  if ( (cwd == NULL) || (goCoordFile->data[0] == '/') )
787  {
788  strcpy(filename, goCoordFile->data);
789  }
790  else
791  {
792  strcpy(filename, cwd);
793  strcat(filename, goCoordFile->data);
794  }
795 
796  goPDB = new PDB(filename);
797  if ( goPDB == NULL )
798  {
799  NAMD_die("Memory allocation failed in Molecule::build_go_sigmas2");
800  }
801 
802  if (goPDB->num_atoms() != numAtoms)
803  {
804  NAMD_die("Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
805  }
806  }
807  // Allocate the array to hold the chain types
809  // Allocate the array to hold Go atom indices into the sigma array
811  // Allocate the array to hold resid
813 
814  if (atomChainTypes == NULL) {
815  NAMD_die("memory allocation failed in Molecule::build_go_sigmas2");
816  }
817 
818  numGoAtoms = 0;
819 
820  // Loop through all the atoms and get the Go chain types
821  for (i=0; i<numAtoms; i++) {
822  // Get the chainType from the occupancy field
823  chainType = (int32)(goPDB->atom(i))->occupancy();
824  // Get the resid from the resid field
825  residType = (int32)(goPDB->atom(i))->residueseq();
826  // Assign the chainType value
827  if ( chainType != 0 ) {
828  //DebugM(3,"build_go_sigmas2 - atom:" << i << ", chainType:" << chainType << std::endl);
829  atomChainTypes[i] = chainType;
831  goResidIndices[i] = residType;
832  numGoAtoms++;
833  }
834  else {
835  atomChainTypes[i] = 0;
836  goSigmaIndices[i] = -1;
837  goResidIndices[i] = -1;
838  }
839  }
840 
841  //Define:
842  ResizeArray<GoPair> tmpGoPair;
843 
844  // Loop through all atom pairs and calculate sigmas
845  // Store sigmas into sorted array
846  DebugM(3," numAtoms=" << numAtoms << std::endl);
847  for (i=0; i<numAtoms; i++) {
848  resid1 = (goPDB->atom(i))->residueseq();
849  for (j=i+1; j<numAtoms; j++) {
850  resid2 = (goPDB->atom(j))->residueseq();
851  if ( goSigmaIndices[i] != -1 && goSigmaIndices[j] != -1 && !atoms_1to4(i,j) ) {
852  residDiff = resid2 - resid1;
853  if (residDiff < 0) residDiff = -residDiff;
854  if ( atomChainTypes[i] && atomChainTypes[j] &&
855  !(this->go_restricted(atomChainTypes[i],atomChainTypes[j],residDiff)) ) {
856  atomAtomDist = sqrt(pow((goPDB->atom(i))->xcoor() - (goPDB->atom(j))->xcoor(), 2.0) +
857  pow((goPDB->atom(i))->ycoor() - (goPDB->atom(j))->ycoor(), 2.0) +
858  pow((goPDB->atom(i))->zcoor() - (goPDB->atom(j))->zcoor(), 2.0));
859  if ( atomAtomDist <= this->get_go_cutoff(atomChainTypes[i],atomChainTypes[j]) ) {
860  exp_a = this->get_go_exp_a(atomChainTypes[i],atomChainTypes[j]);
861  exp_b = this->get_go_exp_b(atomChainTypes[i],atomChainTypes[j]);
862  sigma = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
863  double tmpA = pow(sigma,exp_a);
864  double tmpB = pow(sigma,exp_b);
865  GoPair gp;
866  GoPair gp2;
867  gp.goIndxA = i;
868  gp.goIndxB = j;
869  gp.A = tmpA;
870  gp.B = tmpB;
871  tmpGoPair.add(gp);
872  gp2.goIndxA = j;
873  gp2.goIndxB = i;
874  gp2.A = tmpA;
875  gp2.B = tmpB;
876  tmpGoPair.add(gp2);
877  nativeContacts++;
878  } else {
879  nonnativeContacts++;
880  }
881  }
882  }
883  }
884  }
885 
886  iout << iINFO << "Number of UNIQUE native contacts: " << nativeContacts << "\n" << endi;
887  iout << iINFO << "Number of UNIQUE nonnative contacts: " << nonnativeContacts << "\n" << endi;
888 
889  // Copies the resizeArray into a block of continuous memory
890  std::sort(tmpGoPair.begin(),tmpGoPair.end(),goPairCompare);
891  goNumLJPair = 2*nativeContacts;
892  goIndxLJA = new int[goNumLJPair];
893  goIndxLJB = new int[goNumLJPair];
894  goSigmaPairA = new double[goNumLJPair];
895  goSigmaPairB = new double[goNumLJPair];
896  for(i=0; i< goNumLJPair; i++) {
897  goIndxLJA[i] = tmpGoPair[i].goIndxA;
898  goIndxLJB[i] = tmpGoPair[i].goIndxB;
899  goSigmaPairA[i] = tmpGoPair[i].A;
900  goSigmaPairB[i] = tmpGoPair[i].B;
901  }
902 
903  pointerToGoBeg = new int[numAtoms];
904  pointerToGoEnd = new int[numAtoms];
905  int oldIndex = -1;
906  for(i=0; i<numAtoms; i++) {
907  pointerToGoBeg[i] = -1;
908  pointerToGoEnd[i] = -2;
909  }
910  for(i=0; i < goNumLJPair; i++) {
911  if(pointerToGoBeg[goIndxLJA[i]] == -1) {
912  pointerToGoBeg[goIndxLJA[i]] = i;
913  oldIndex = goIndxLJA[i];
914  }
915  pointerToGoEnd[oldIndex] = i;
916  }
917 
918  // If we had to create a PDB object, delete it now
919  if (goCoordFile != NULL) {
920  delete goPDB;
921  }
922 
923  return;
924 }
int get_go_exp_a(int chain1, int chain2)
Definition: Molecule.h:1613
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
Real get_go_cutoff(int chain1, int chain2)
Definition: Molecule.h:1601
Definition: PDB.h:35
double A
Definition: Molecule.h:120
int * pointerToGoBeg
Definition: Molecule.h:657
PDB * goPDB
Definition: Molecule.h:650
short int32
Definition: dumpdcd.c:24
int goIndxB
Definition: Molecule.h:119
float Real
Definition: common.h:107
#define DebugM(x, y)
Definition: Debug.h:59
double * goSigmaPairA
Definition: Molecule.h:655
#define iout
Definition: InfoStream.h:87
int num_atoms(void)
Definition: PDB.C:323
int32 * atomChainTypes
Definition: Molecule.h:643
double B
Definition: Molecule.h:121
int goIndxA
Definition: Molecule.h:118
int * pointerToGoEnd
Definition: Molecule.h:658
int * goIndxLJA
Definition: Molecule.h:653
PDBAtom * atom(int place)
Definition: PDB.C:394
int goNumLJPair
Definition: Molecule.h:652
iterator end(void)
Definition: ResizeArray.h:37
int32 * goResidIndices
Definition: Molecule.h:645
Bool go_restricted(int, int, int)
Definition: GoMolecule.C:525
int numAtoms
Definition: Molecule.h:556
void NAMD_die(const char *err_msg)
Definition: common.C:83
double * goSigmaPairB
Definition: Molecule.h:656
int add(const Elem &elem)
Definition: ResizeArray.h:97
BlockRadixSort::TempStorage sort
int get_go_exp_b(int chain1, int chain2)
Definition: Molecule.h:1616
StringList * next
Definition: ConfigList.h:49
char * data
Definition: ConfigList.h:48
infostream & endi(infostream &s)
Definition: InfoStream.C:38
Bool atoms_1to4(unsigned int, unsigned int)
Definition: GoMolecule.C:1537
int32 * goSigmaIndices
Definition: Molecule.h:644
int * goIndxLJB
Definition: Molecule.h:654
int numGoAtoms
Definition: Molecule.h:642
int numLJPair
Definition: Molecule.h:661
iterator begin(void)
Definition: ResizeArray.h:36
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(), GridforceGrid::new_grid(), MGridforceParams::next, and PDB::num_atoms().

Referenced by NamdState::loadStructure().

6379 {
6380  PDB *kPDB;
6381  register int i; // Loop counters
6382  register int j;
6383  register int k;
6384 
6385  DebugM(3, "Entered build_gridforce_params multi...\n");
6386 // DebugM(3, "\tgridfrcfile = " << gridfrcfile->data << endi);
6387 // DebugM(3, "\tgridfrccol = " << gridfrccol->data << endi);
6388 
6389  MGridforceParams* mgridParams = simParams->mgridforcelist.get_first();
6390  numGridforceGrids = 0;
6391  while (mgridParams != NULL) {
6393  mgridParams = mgridParams->next;
6394  }
6395 
6396  DebugM(3, "numGridforceGrids = " << numGridforceGrids << "\n");
6397  gridfrcIndexes = new int32*[numGridforceGrids];
6398  gridfrcParams = new GridforceParams*[numGridforceGrids];
6399  gridfrcGrid = new GridforceGrid*[numGridforceGrids];
6400  numGridforces = new int[numGridforceGrids];
6401 
6402  int grandTotalGrids = 0; // including all subgrids
6403 
6404  mgridParams = simParams->mgridforcelist.get_first();
6405  for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
6406  int current_index=0; // Index into values used
6407  int kcol = 5; // Column to look for force constant in
6408  int qcol = 0; // Column for charge (default 0: use electric charge)
6409  Real kval = 0; // Force constant value retreived
6410  char filename[129]; // PDB filename
6411  char potfilename[129]; // Potential file name
6412 
6413  if (mgridParams == NULL) {
6414  NAMD_die("Problem with mgridParams!");
6415  }
6416 
6417  // Now load values from mgridforcelist object
6418  if (mgridParams->gridforceFile == NULL)
6419  {
6420  if ( ! initial_pdb ) NAMD_die("Initial PDB file unavailable, gridforceFile required.");
6421  kPDB = initial_pdb;
6422  }
6423  else
6424  {
6425  DebugM(4, "mgridParams->gridforceFile = " << mgridParams->gridforceFile << "\n" << endi);
6426 
6427  if ( (cwd == NULL) || (mgridParams->gridforceFile[0] == '/') )
6428  {
6429  strcpy(filename, mgridParams->gridforceFile);
6430  }
6431  else
6432  {
6433  strcpy(filename, cwd);
6434  strcat(filename, mgridParams->gridforceFile);
6435  }
6436 
6437  kPDB = new PDB(filename);
6438  if ( kPDB == NULL )
6439  {
6440  NAMD_die("Memory allocation failed in Molecule::build_gridforce_params");
6441  }
6442 
6443  if (kPDB->num_atoms() != numAtoms)
6444  {
6445  NAMD_die("Number of atoms in grid force PDB doesn't match coordinate PDB");
6446  }
6447  }
6448 
6449  // Get the column that the force constant is going to be in. It
6450  // can be in any of the 5 floating point fields in the PDB, according
6451  // to what the user wants. The allowable fields are X, Y, Z, O, or
6452  // B which correspond to the 1st, 2nd, ... 5th floating point fields.
6453  // The default is the 5th field, which is beta (temperature factor)
6454  if (mgridParams->gridforceCol == NULL)
6455  {
6456  kcol = 5;
6457  }
6458  else
6459  {
6460  if (strcasecmp(mgridParams->gridforceCol, "X") == 0)
6461  {
6462  kcol=1;
6463  }
6464  else if (strcasecmp(mgridParams->gridforceCol, "Y") == 0)
6465  {
6466  kcol=2;
6467  }
6468  else if (strcasecmp(mgridParams->gridforceCol, "Z") == 0)
6469  {
6470  kcol=3;
6471  }
6472  else if (strcasecmp(mgridParams->gridforceCol, "O") == 0)
6473  {
6474  kcol=4;
6475  }
6476  else if (strcasecmp(mgridParams->gridforceCol, "B") == 0)
6477  {
6478  kcol=5;
6479  }
6480  else
6481  {
6482  NAMD_die("gridforcecol must have value of X, Y, Z, O, or B");
6483  }
6484  }
6485 
6486  // Get the column that the charge is going to be in.
6487  if (mgridParams->gridforceQcol == NULL)
6488  {
6489  qcol = 0; // Default: don't read charge from file, use electric charge
6490  }
6491  else
6492  {
6493  if (strcasecmp(mgridParams->gridforceQcol, "X") == 0)
6494  {
6495  qcol=1;
6496  }
6497  else if (strcasecmp(mgridParams->gridforceQcol, "Y") == 0)
6498  {
6499  qcol=2;
6500  }
6501  else if (strcasecmp(mgridParams->gridforceQcol, "Z") == 0)
6502  {
6503  qcol=3;
6504  }
6505  else if (strcasecmp(mgridParams->gridforceQcol, "O") == 0)
6506  {
6507  qcol=4;
6508  }
6509  else if (strcasecmp(mgridParams->gridforceQcol, "B") == 0)
6510  {
6511  qcol=5;
6512  }
6513  else
6514  {
6515  NAMD_die("gridforcechargecol must have value of X, Y, Z, O, or B");
6516  }
6517  }
6518 
6519  if (kcol == qcol) {
6520  NAMD_die("gridforcecol and gridforcechargecol cannot have same value");
6521  }
6522 
6523 
6524  // Allocate an array that will store an index into the constraint
6525  // parameters for each atom. If the atom is not constrained, its
6526  // value will be set to -1 in this array.
6527  gridfrcIndexes[gridnum] = new int32[numAtoms];
6528 
6529  if (gridfrcIndexes[gridnum] == NULL)
6530  {
6531  NAMD_die("memory allocation failed in Molecule::build_gridforce_params()");
6532  }
6533 
6534  // Loop through all the atoms and find out which ones are constrained
6535  for (i=0; i<numAtoms; i++)
6536  {
6537  // Get the k value based on where we were told to find it
6538  switch (kcol)
6539  {
6540  case 1:
6541  kval = (kPDB->atom(i))->xcoor();
6542  break;
6543  case 2:
6544  kval = (kPDB->atom(i))->ycoor();
6545  break;
6546  case 3:
6547  kval = (kPDB->atom(i))->zcoor();
6548  break;
6549  case 4:
6550  kval = (kPDB->atom(i))->occupancy();
6551  break;
6552  case 5:
6553  kval = (kPDB->atom(i))->temperaturefactor();
6554  break;
6555  }
6556 
6557  if (kval > 0.0)
6558  {
6559  // This atom is constrained
6560  gridfrcIndexes[gridnum][i] = current_index;
6561  current_index++;
6562  }
6563  else
6564  {
6565  // This atom is not constrained
6566  gridfrcIndexes[gridnum][i] = -1;
6567  }
6568  }
6569 
6570  if (current_index == 0)
6571  {
6572  // Constraints were turned on, but there weren't really any constrained
6573  iout << iWARN << "NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" << endi;
6574  }
6575  else
6576  {
6577  // Allocate an array to hold the constraint parameters
6578  gridfrcParams[gridnum] = new GridforceParams[current_index];
6579  if (gridfrcParams[gridnum] == NULL)
6580  {
6581  NAMD_die("memory allocation failed in Molecule::build_gridforce_params");
6582  }
6583  }
6584 
6585  numGridforces[gridnum] = current_index;
6586 
6587  // Loop through all the atoms and assign the parameters for those
6588  // that are constrained
6589  for (i=0; i<numAtoms; i++)
6590  {
6591  if (gridfrcIndexes[gridnum][i] != -1)
6592  {
6593  // This atom has grid force, so get the k value again
6594  switch (kcol)
6595  {
6596  case 1:
6597  gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->xcoor();
6598  break;
6599  case 2:
6600  gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->ycoor();
6601  break;
6602  case 3:
6603  gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->zcoor();
6604  break;
6605  case 4:
6606  gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->occupancy();
6607  break;
6608  case 5:
6609  gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->temperaturefactor();
6610  break;
6611  }
6612 
6613  // Also get charge column
6614  switch (qcol)
6615  {
6616  case 0:
6617 #ifdef MEM_OPT_VERSION
6618  gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]];
6619 #else
6620  gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge;
6621 #endif
6622  break;
6623  case 1:
6624  gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->xcoor();
6625  break;
6626  case 2:
6627  gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->ycoor();
6628  break;
6629  case 3:
6630  gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->zcoor();
6631  break;
6632  case 4:
6633  gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->occupancy();
6634  break;
6635  case 5:
6636  gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->temperaturefactor();
6637  break;
6638  }
6639  }
6640  }
6641 
6642  // If we had to create new PDB objects, delete them now
6643  if (mgridParams->gridforceFile != NULL)
6644  {
6645  delete kPDB;
6646  }
6647 
6648  // Now we fill in our grid information
6649 
6650  // Open potential file
6651  if ( (cwd == NULL) || (mgridParams->gridforceVfile[0] == '/') )
6652  {
6653  strcpy(potfilename, mgridParams->gridforceVfile);
6654  }
6655  else
6656  {
6657  strcpy(potfilename, cwd);
6658  strcat(potfilename, mgridParams->gridforceVfile);
6659  }
6660 
6661 // iout << iINFO << "Allocating grid " << gridnum
6662 // << "\n" << endi;
6663 
6664  DebugM(3, "allocating GridforceGrid(" << gridnum << ")\n");
6665  gridfrcGrid[gridnum] = GridforceGrid::new_grid(gridnum, potfilename, simParams, mgridParams);
6666 
6667  grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
6668  DebugM(4, "grandTotalGrids = " << grandTotalGrids << "\n" << endi);
6669 
6670  // Finally, get next mgridParams pointer
6671  mgridParams = mgridParams->next;
6672  }
6673 }
Definition: PDB.h:35
short int32
Definition: dumpdcd.c:24
static GridforceGrid * new_grid(int gridnum, char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)
Definition: GridForceGrid.C:34
float Real
Definition: common.h:107
#define DebugM(x, y)
Definition: Debug.h:59
int numGridforceGrids
Definition: Molecule.h:590
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
#define iout
Definition: InfoStream.h:87
int num_atoms(void)
Definition: PDB.C:323
MGridforceParams * get_first()
PDBAtom * atom(int place)
Definition: PDB.C:394
int numAtoms
Definition: Molecule.h:556
void NAMD_die(const char *err_msg)
Definition: common.C:83
MGridforceParamsList mgridforcelist
MGridforceParams * next
int * numGridforces
Definition: Molecule.h:591
infostream & endi(infostream &s)
Definition: InfoStream.C:38
virtual int get_total_grids(void) const =0
void Molecule::build_gro_pair ( )
void Molecule::build_langevin_params ( BigReal  coupling,
BigReal  drudeCoupling,
Bool  doHydrogen 
)
void Molecule::build_langevin_params ( StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_movdrag_params ( StringList ,
StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_rotdrag_params ( StringList ,
StringList ,
StringList ,
StringList ,
StringList ,
StringList ,
PDB ,
char *   
)
void Molecule::build_ss_flags ( const StringList ssfile,
const StringList sscol,
PDB initial_pdb,
const char *  cwd 
)

Build the flags needed for solute scaling.

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

Parameters
ssfileconfig "soluteScalingFile = my.pdb" for PDB filename
sscolconfig "soluteScalingCol = O" for column of PDB ATOM records
initial_pdbthe initial PDB file
cwdcurrent working directory

Referenced by NamdState::loadStructure().

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

Definition at line 1547 of file MoleculeQM.C.

References 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.

1547  {
1548 
1549 #ifdef MEM_OPT_VERSION
1550  NAMD_die("QMMM interface is not supported in memory optimized builds.");
1551 #else
1552 
1553  DebugM(3,"Cleaning QM bonds, angles, dihedrals, impropers and crossterms.\n")
1554 
1555  // Bonds
1556  Bond *nonQMBonds;
1557  nonQMBonds = new Bond[numBonds] ;
1558  int nonQMBondCount = 0;
1559  qmDroppedBonds = 0;
1560  for (int i = 0; i < numBonds; i++) {
1561 
1562  int part1 = qmAtomGroup[bonds[i].atom1];
1563  int part2 = qmAtomGroup[bonds[i].atom2];
1564  if (part1 > 0 && part2 > 0 ) {
1565 
1566  // If the user defined extra bonds, we check if the QM-QM bond is an "extra"
1567  // bond, and do not delete it.
1568  if (simParams->extraBondsOn) {
1569 // std::cout << "Checking Bond: " << bonds[i].atom1 << "," << bonds[i].atom2 << std::endl ;
1570 
1571  for (int ebi=0; ebi < qmExtraBonds.size() ; ebi++) {
1572 
1573  if( (qmExtraBonds[ebi].atom1 == bonds[i].atom1 ||
1574  qmExtraBonds[ebi].atom1 == bonds[i].atom2) &&
1575  (qmExtraBonds[ebi].atom2 == bonds[i].atom1 ||
1576  qmExtraBonds[ebi].atom2 == bonds[i].atom2) ) {
1577 // std::cout << "This is an extra bond! We will keep it." << std::endl;
1578  nonQMBonds[nonQMBondCount++] = bonds[i];
1579  break;
1580  }
1581  }
1582  }
1583 
1584  qmDroppedBonds++;
1585  } else {
1586  // Just a simple sanity check.
1587  // If there are no QM bonds defined by the user but we find a
1588  // bond between a QM and an MM atom, error out.
1589  if ((part1 > 0 || part2 > 0) && qmNumBonds == 0) {
1590  iout << iERROR << " Atoms " << bonds[i].atom1
1591  << " and " << bonds[i].atom2 << " form a QM-MM bond!\n" << endi ;
1592  NAMD_die("No QM-MM bonds were defined by the user, but one was found.");
1593  }
1594  nonQMBonds[nonQMBondCount++] = bonds[i];
1595  }
1596  }
1597  numBonds = nonQMBondCount;
1598  delete [] bonds;
1599  bonds = new Bond[numBonds];
1600  for (int i = 0; i < nonQMBondCount; i++) {
1601  bonds[i]=nonQMBonds[i];
1602  }
1603  delete [] nonQMBonds;
1605 
1606  // Angles
1607  Angle *nonQMAngles;
1608  nonQMAngles = new Angle[numAngles];
1609  int nonQMAngleCount = 0;
1610  qmDroppedAngles = 0;
1611  for (int i = 0; i < numAngles; i++) {
1612  int part1 = qmAtomGroup[angles[i].atom1];
1613  int part2 = qmAtomGroup[angles[i].atom2];
1614  int part3 = qmAtomGroup[angles[i].atom3];
1615  if (part1 > 0 && part2 > 0 && part3 > 0) {
1616  //CkPrintf("-----ANGLE ATOMS %i %i %i partitions %i %i %i\n",angles[i].atom1, angles[i].atom2, angles[i].atom3, part1, part2, part3);
1617  qmDroppedAngles++;
1618  }
1619  else {
1620  nonQMAngles[nonQMAngleCount++] = angles[i];
1621  }
1622  }
1623  numAngles = nonQMAngleCount;
1624  delete [] angles;
1625  angles = new Angle[numAngles];
1626  for (int i = 0; i < nonQMAngleCount; i++) {
1627  angles[i]=nonQMAngles[i];
1628  }
1629  delete [] nonQMAngles;
1630 
1631 
1632  // Dihedrals
1633  Dihedral *nonQMDihedrals;
1634  nonQMDihedrals = new Dihedral[numDihedrals];
1635  int nonQMDihedralCount = 0;
1636  qmDroppedDihedrals = 0;
1637  for (int i = 0; i < numDihedrals; i++) {
1638  int part1 = qmAtomGroup[dihedrals[i].atom1];
1639  int part2 = qmAtomGroup[dihedrals[i].atom2];
1640  int part3 = qmAtomGroup[dihedrals[i].atom3];
1641  int part4 = qmAtomGroup[dihedrals[i].atom4];
1642  if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0) {
1643  //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);
1644  qmDroppedDihedrals++;
1645  }
1646  else {
1647  nonQMDihedrals[nonQMDihedralCount++] = dihedrals[i];
1648  }
1649  }
1650  numDihedrals = nonQMDihedralCount;
1651  delete [] dihedrals;
1652  dihedrals = new Dihedral[numDihedrals];
1653  for (int i = 0; i < numDihedrals; i++) {
1654  dihedrals[i]=nonQMDihedrals[i];
1655  }
1656  delete [] nonQMDihedrals;
1657 
1658  // Impropers
1659  Improper *nonQMImpropers;
1660  nonQMImpropers = new Improper[numImpropers];
1661  int nonQMImproperCount = 0;
1662  qmDroppedImpropers = 0;
1663  for (int i = 0; i < numImpropers; i++) {
1664  int part1 = qmAtomGroup[impropers[i].atom1];
1665  int part2 = qmAtomGroup[impropers[i].atom2];
1666  int part3 = qmAtomGroup[impropers[i].atom3];
1667  int part4 = qmAtomGroup[impropers[i].atom4];
1668  if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0) {
1669  //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);
1670  qmDroppedImpropers++;
1671  }
1672  else {
1673  nonQMImpropers[nonQMImproperCount++] = impropers[i];
1674  }
1675  }
1676  numImpropers = nonQMImproperCount;
1677  delete [] impropers;
1678  impropers = new Improper[numImpropers];
1679  for (int i = 0; i < numImpropers; i++) {
1680  impropers[i]=nonQMImpropers[i];
1681  }
1682  delete [] nonQMImpropers;
1683 
1684  // Crossterms
1685  Crossterm *nonQMCrossterms;
1686  nonQMCrossterms = new Crossterm[numCrossterms];
1687  int nonQMCrosstermCount = 0;
1688  qmDroppedCrossterms = 0 ;
1689  for (int i = 0; i < numCrossterms; i++) {
1690  int part1 = qmAtomGroup[crossterms[i].atom1];
1691  int part2 = qmAtomGroup[crossterms[i].atom2];
1692  int part3 = qmAtomGroup[crossterms[i].atom3];
1693  int part4 = qmAtomGroup[crossterms[i].atom4];
1694  int part5 = qmAtomGroup[crossterms[i].atom5];
1695  int part6 = qmAtomGroup[crossterms[i].atom6];
1696  int part7 = qmAtomGroup[crossterms[i].atom7];
1697  int part8 = qmAtomGroup[crossterms[i].atom8];
1698  if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0 &&
1699  part5 > 0 && part6 > 0 && part7 > 0 && part8 > 0) {
1700  //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);
1701  qmDroppedCrossterms++;
1702  }
1703  else {
1704  nonQMCrossterms[nonQMCrosstermCount++] = crossterms[i];
1705  }
1706  }
1707  numCrossterms = nonQMCrosstermCount;
1708  delete [] crossterms;
1709  crossterms = new Crossterm[numCrossterms] ;
1710  for (int i=0; i < numCrossterms; i++){
1711  crossterms[i] = nonQMCrossterms[i] ;
1712  }
1713  delete [] nonQMCrossterms ;
1714 
1715  if (!CkMyPe()) {
1716  iout << iINFO << "The QM region will remove " << qmDroppedBonds << " bonds, " <<
1717  qmDroppedAngles << " angles, " << qmDroppedDihedrals << " dihedrals, "
1718  << qmDroppedImpropers << " impropers and " << qmDroppedCrossterms
1719  << " crossterms.\n" << endi ;
1720  }
1721 
1722 #endif
1723 }
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
int numBonds
Definition: Molecule.h:559
#define DebugM(x, y)
Definition: Debug.h:59
int numRealBonds
Definition: Molecule.h:558
int32 atom1
Definition: structures.h:81
#define iout
Definition: InfoStream.h:87
int32 atom1
Definition: structures.h:48
int32 atom1
Definition: structures.h:55
int numAngles
Definition: Molecule.h:560
int numCrossterms
Definition: Molecule.h:567
void NAMD_die(const char *err_msg)
Definition: common.C:83
int numDihedrals
Definition: Molecule.h:561
int numImpropers
Definition: Molecule.h:566
int32 atom1
Definition: structures.h:72
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:109
infostream & endi(infostream &s)
Definition: InfoStream.C:38
int32 atom1
Definition: structures.h:63
for(int i=0;i< n1;++i)
void Molecule::freeBFactorData ( )
inline

Definition at line 1034 of file Molecule.h.

Referenced by NamdState::loadStructure().

1034 { delete [] bfactor; bfactor=NULL; }
void Molecule::freeOccupancyData ( )
inline

Definition at line 1030 of file Molecule.h.

Referenced by NamdState::loadStructure().

1030 { delete [] occupancy; occupancy=NULL; }
Bond* Molecule::get_acceptor ( int  dnum) const
inline

Definition at line 1105 of file Molecule.h.

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

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

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

159  {
160 
161  if (atomNames == NULL || resLookup == NULL)
162  {
163  NAMD_die("Tried to find atom from name on node other than node 0");
164  }
165  int i = 0;
166  int end = 0;
167  if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
168  if ( index >= 0 && index < ( end - i ) ) return ( index + i );
169  return -1;
170 }
int lookup(const char *segid, int resid, int *begin, int *end) const
Definition: Molecule.C:76
void NAMD_die(const char *err_msg)
Definition: common.C:83
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().

122  {
123 
124  if (atomNames == NULL || resLookup == NULL)
125  {
126  NAMD_die("Tried to find atom from name on node other than node 0");
127  }
128 
129  int i = 0;
130  int end = 0;
131  if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
132  for ( ; i < end; ++i ) {
133  #ifdef MEM_OPT_VERSION
134  Index idx = atomNames[i].atomnameIdx;
135  if(!strcasecmp(aname, atomNamePool[idx])) return i;
136  #else
137  if ( ! strcasecmp(aname,atomNames[i].atomname) ) return i;
138  #endif
139  }
140  return -1;
141 }
int Index
Definition: structures.h:26
int lookup(const char *segid, int resid, int *begin, int *end) const
Definition: Molecule.C:76
void NAMD_die(const char *err_msg)
Definition: common.C:83
HashPool< HashString > atomNamePool
Definition: CompressPsf.C:309
const char* Molecule::get_atomtype ( int  anum) const
inline

Definition at line 1116 of file Molecule.h.

References atomTypePool, and NAMD_die().

1117  {
1118  if (atomNames == NULL)
1119  {
1120  NAMD_die("Tried to find atom type on node other than node 0");
1121  }
1122 
1123  #ifdef MEM_OPT_VERSION
1124  return atomTypePool[atomNames[anum].atomtypeIdx];
1125  #else
1126  return(atomNames[anum].atomtype);
1127  #endif
1128  }
void NAMD_die(const char *err_msg)
Definition: common.C:83
HashPool< HashString > atomTypePool
Definition: CompressPsf.C:310
Bond* Molecule::get_bond ( int  bnum) const
inline

Definition at line 1065 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

1146  { return bondsByAtom[anum]; }
int Molecule::get_cluster ( int  anum) const
inline

Definition at line 1023 of file Molecule.h.

Referenced by wrap_coor_int().

1023 { return cluster[anum]; }
int Molecule::get_clusterSize ( int  anum) const
inline

Definition at line 1024 of file Molecule.h.

Referenced by wrap_coor_int().

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

1263  {
1264  k = consParams[consIndexes[atomnum]].k;
1265  refPos = consParams[consIndexes[atomnum]].refPos;
1266  }
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().

1338  {
1339  v = consTorqueParams[consTorqueIndexes[atomnum]].v;
1340  a = consTorqueParams[consTorqueIndexes[atomnum]].a;
1341  p = consTorqueParams[consTorqueIndexes[atomnum]].p;
1342  }
int32 * consTorqueIndexes
Definition: Molecule.h:615
ConsTorqueParams * consTorqueParams
Definition: Molecule.h:616
Crossterm* Molecule::get_crossterm ( int  inum) const
inline

Definition at line 1077 of file Molecule.h.

1077 {return (&(crossterms[inum]));}
int32* Molecule::get_crossterms_for_atom ( int  anum)
inline

Definition at line 1153 of file Molecule.h.

1154  { return crosstermsByAtom[anum]; }
const Real* Molecule::get_cSMDcoffs ( )
inline

Definition at line 874 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

874 {return cSMDcoffs;} ;
int const* const* Molecule::get_cSMDindex ( )
inline

Definition at line 870 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

870 {return cSMDindex;} ;
int const* Molecule::get_cSMDindxLen ( )
inline

Definition at line 869 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

869 { return cSMDindxLen;} ;
const Real* Molecule::get_cSMDKs ( )
inline

Definition at line 872 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

872 {return cSMDKs;} ;
int Molecule::get_cSMDnumInst ( )
inline

Definition at line 868 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

868 { return cSMDnumInst;} ;
int const* const* Molecule::get_cSMDpairs ( )
inline

Definition at line 871 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

871 {return cSMDpairs;} ;
const Real* Molecule::get_cSMDVels ( )
inline

Definition at line 873 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

873 {return cSMDVels;} ;
Dihedral* Molecule::get_dihedral ( int  dnum) const
inline

Definition at line 1074 of file Molecule.h.

Referenced by atoms_1to4(), and dumpbench().

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

1150  { return dihedralsByAtom[anum]; }
Bond* Molecule::get_donor ( int  dnum) const
inline

Definition at line 1102 of file Molecule.h.

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

1173  {
1174  return &all_exclusions[anum];
1175  }
Exclusion* Molecule::get_exclusion ( int  ex) const
inline

Definition at line 1112 of file Molecule.h.

References exclusions.

1112 {return (&(exclusions[ex]));}
int32* Molecule::get_exclusions_for_atom ( int  anum)
inline

Definition at line 1155 of file Molecule.h.

References exclusionsByAtom.

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

1386  {
1387  int typeSum = 0;
1388  for ( int i=0; i < order; ++i ) {
1389  typeSum += (fepAtomFlags[atomID[i]] == 2 ? -1 : fepAtomFlags[atomID[i]]);
1390  }
1391  // Increase the cutoff if scaling purely alchemical bonds.
1392  // This really only applies when order = 2.
1393  if ( simParams->alchBondDecouple ) order++;
1394 
1395  // The majority of interactions are type 0, so catch those first.
1396  if ( typeSum == 0 || abs(typeSum) == order ) return 0;
1397  else if ( 0 < typeSum && typeSum < order ) return 1;
1398  else if ( -order < typeSum && typeSum < 0 ) return 2;
1399 
1400  // Alchemify should always keep this from bombing, but just in case...
1401  NAMD_die("Unexpected alchemical bonded interaction!");
1402  return 0;
1403  }
#define order
Definition: PmeRealSpace.C:235
void NAMD_die(const char *err_msg)
Definition: common.C:83
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().

1346  {
1347  return(fepAtomFlags[anum]);
1348  }
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().

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

1601 { return go_array[MAX_GO_CHAINS*chain1 + chain2].cutoff; };
#define MAX_GO_CHAINS
Definition: Molecule.h:36
GoValue go_array[MAX_GO_CHAINS *MAX_GO_CHAINS]
Definition: Molecule.h:1584
Real cutoff
Definition: Molecule.h:112
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().

1610 { return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilon; };
#define MAX_GO_CHAINS
Definition: Molecule.h:36
GoValue go_array[MAX_GO_CHAINS *MAX_GO_CHAINS]
Definition: Molecule.h:1584
Real epsilon
Definition: Molecule.h:106
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().

1604 { return go_array[MAX_GO_CHAINS*chain1 + chain2].epsilonRep; };
Real epsilonRep
Definition: Molecule.h:111
#define MAX_GO_CHAINS
Definition: Molecule.h:36
GoValue go_array[MAX_GO_CHAINS *MAX_GO_CHAINS]
Definition: Molecule.h:1584
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().

1613 { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_a; };
#define MAX_GO_CHAINS
Definition: Molecule.h:36
int exp_a
Definition: Molecule.h:107
GoValue go_array[MAX_GO_CHAINS *MAX_GO_CHAINS]
Definition: Molecule.h:1584
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().

1616 { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_b; };
int exp_b
Definition: Molecule.h:108
#define MAX_GO_CHAINS
Definition: Molecule.h:36
GoValue go_array[MAX_GO_CHAINS *MAX_GO_CHAINS]
Definition: Molecule.h:1584
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().

1619 { return go_array[MAX_GO_CHAINS*chain1 + chain2].exp_rep; };
#define MAX_GO_CHAINS
Definition: Molecule.h:36
GoValue go_array[MAX_GO_CHAINS *MAX_GO_CHAINS]
Definition: Molecule.h:1584
int exp_rep
Definition: Molecule.h:109
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.

1265 {
1266  BigReal goForce = 0.0;
1267  Real pow1;
1268  Real pow2;
1269  // determine which Go chain pair we are working with
1270  //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ")" << std::endl);
1271  int32 chain1 = atomChainTypes[atom1];
1272  int32 chain2 = atomChainTypes[atom2];
1273 
1274  //DebugM(3," chain1:" << chain1 << ", chain2:" << chain2 << std::endl);
1275  if (chain1 == 0 || chain2 == 0) return 0.0;
1276 
1277  // retrieve Go cutoff for this chain pair
1278  //TMP// JLai -- I'm going to replace this with a constant accessor. This is just a temporary thing
1279  Real goCutoff = const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2);
1280  //DebugM(3," goCutoff:" << goCutoff << std::endl);
1281  if (goCutoff == 0) return 0.0;
1282  // if repulsive then calculate repulsive
1283  // sigmas are initially set to -1.0 if the atom pair fails go_restricted
1284  if (goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]] != -1.0) {
1285  if (!goWithinCutoff[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]]) {
1286  Real epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1,chain2);
1287  Real sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1,chain2);
1288  int exp_rep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
1289  pow1 = pow(sigmaRep/r,exp_rep);
1290  goForce = 4*((exp_rep/(r*r)) * epsilonRep * pow1);
1291  *goNative = 0.0;
1292  *goNonnative = (4 * epsilonRep * pow1 );
1293  //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", epsilonRep:" << epsilonRep << ", sigmaRep:" << sigmaRep << ", r:" << r << ", goForce:" << goForce << std::endl);
1294  }
1295  // if attractive then calculate attractive
1296  else {
1297  int goSigmaIndex1 = goSigmaIndices[atom1];
1298  int goSigmaIndex2 = goSigmaIndices[atom2];
1299  if (goSigmaIndex1 != -1 && goSigmaIndex2 != -1) {
1300  Real epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
1301  int exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
1302  int exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
1303  Real sigma_ij = goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]];
1304  // Positive gradient of potential, not negative gradient of potential
1305  pow1 = pow(sigma_ij/r,exp_a);
1306  pow2 = pow(sigma_ij/r,exp_b);
1307  goForce = ((4/(r*r)) * epsilon * (exp_a * pow1 - exp_b * pow2));
1308  //DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", sigma_ij:" << sigma_ij << ", r:" << r << ", goForce:" << goForce << std::endl);
1309  *goNative = (4 * epsilon * ( pow1 - pow2 ) );
1310  *goNonnative = 0.0;
1311  }
1312  }
1313  }
1314  //DebugM(3,"goForce:" << goForce << std::endl);
1315  return goForce;
1316 }
int get_go_exp_a(int chain1, int chain2)
Definition: Molecule.h:1613
Real get_go_cutoff(int chain1, int chain2)
Definition: Molecule.h:1601
short int32
Definition: dumpdcd.c:24
Real get_go_epsilonRep(int chain1, int chain2)
Definition: Molecule.h:1604
float Real
Definition: common.h:107
bool * goWithinCutoff
Definition: Molecule.h:647
Real get_go_epsilon(int chain1, int chain2)
Definition: Molecule.h:1610
int32 * atomChainTypes
Definition: Molecule.h:643
Real get_go_sigmaRep(int chain1, int chain2)
Definition: Molecule.h:1607
int get_go_exp_rep(int chain1, int chain2)
Definition: Molecule.h:1619
Real * goSigmas
Definition: Molecule.h:646
BigReal goNonnative
int get_go_exp_b(int chain1, int chain2)
Definition: Molecule.h:1616
BigReal goForce
int32 * goSigmaIndices
Definition: Molecule.h:644
int numGoAtoms
Definition: Molecule.h:642
double BigReal
Definition: common.h:112
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.

1463 {
1464 
1465  // Check to see if restricted. If so, escape function early
1466  int32 chain1 = atomChainTypes[atom1];
1467  int32 chain2 = atomChainTypes[atom2];
1468 
1469  if(chain1 == 0 || chain2 == 0) return 0.0;
1470  Molecule *mol = const_cast<Molecule*>(this);
1471  Real goCutoff = mol->get_go_cutoff(chain1,chain2);
1472  if(goCutoff == 0) return 0.0;
1473 
1474  int resid1 = goResidIndices[atom1];
1475  int resid2 = goResidIndices[atom2];
1476  int residDiff = abs(resid1 - resid2);
1477  if((mol->go_restricted(chain1,chain2,residDiff))) {
1478  return 0.0;
1479  }
1480 
1481  int LJIndex = -1;
1482  int LJbegin = pointerToGoBeg[atom1];
1483  int LJend = pointerToGoEnd[atom1];
1484  for(int i = LJbegin; i <= LJend; i++) {
1485  if(goIndxLJB[i] == atom2) {
1486  LJIndex = i;
1487  }
1488  }
1489 
1490  BigReal r2 = x*x + y*y + z*z;
1491  BigReal r = sqrt(r2);
1492 
1493  if (LJIndex == -1) {
1494  int exp_rep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
1495  BigReal epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1, chain2);
1496  BigReal sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1, chain2);
1497  double sigmaRepPow = pow(sigmaRep,exp_rep);
1498  BigReal LJ = (4*epsilonRep*exp_rep*sigmaRepPow*pow(r,-(exp_rep+1)));
1499  *goNative = 0;
1500  *goNonnative = (4*epsilonRep*sigmaRepPow*pow(r,-exp_rep));
1501  //*goNonnative = (4*epsilonRep * pow(sigmaRep/r,exp_rep));
1502  return (LJ/r);
1503  } else {
1504  // Code to calculate distances because the pair was found in one of the lists
1505  int exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
1506  int exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
1507  // We want the force, so we have to take the n+1 derivative
1508  BigReal powA = pow(r,-(exp_a + 1));
1509  BigReal powB = pow(r,-(exp_b + 1));
1510  BigReal powaa = pow(r,-exp_a);
1511  BigReal powbb = pow(r,-exp_b);
1512  BigReal epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
1513  BigReal LJ = 4 * epsilon * (exp_a*goSigmaPairA[LJIndex]*powA - exp_b*goSigmaPairB[LJIndex]*powB);
1514  *goNative = 4 * epsilon * (goSigmaPairA[LJIndex]*powaa - goSigmaPairB[LJIndex]*powbb);
1515  *goNonnative = 0;
1516  return (LJ/r);
1517  }
1518  return 0;
1519 }
int get_go_exp_a(int chain1, int chain2)
Definition: Molecule.h:1613
Real get_go_cutoff(int chain1, int chain2)
Definition: Molecule.h:1601
int * pointerToGoBeg
Definition: Molecule.h:657
short int32
Definition: dumpdcd.c:24
Real get_go_epsilonRep(int chain1, int chain2)
Definition: Molecule.h:1604
float Real
Definition: common.h:107
double * goSigmaPairA
Definition: Molecule.h:655
Real get_go_epsilon(int chain1, int chain2)
Definition: Molecule.h:1610
int32 * atomChainTypes
Definition: Molecule.h:643
int * pointerToGoEnd
Definition: Molecule.h:658
Real get_go_sigmaRep(int chain1, int chain2)
Definition: Molecule.h:1607
int get_go_exp_rep(int chain1, int chain2)
Definition: Molecule.h:1619
int32 * goResidIndices
Definition: Molecule.h:645
gridSize z
Bool go_restricted(int, int, int)
Definition: GoMolecule.C:525
BigReal goNonnative
double * goSigmaPairB
Definition: Molecule.h:656
int get_go_exp_b(int chain1, int chain2)
Definition: Molecule.h:1616
gridSize y
gridSize x
int * goIndxLJB
Definition: Molecule.h:654
double BigReal
Definition: common.h:112
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.

1339 {
1340  int resid1;
1341  int resid2;
1342  int residDiff;
1343  Real xcoorDiff;
1344  Real ycoorDiff;
1345  Real zcoorDiff;
1346  Real atomAtomDist;
1347  Real exp_a;
1348  Real exp_b;
1349  Real sigma_ij;
1350  Real epsilon;
1351  Real epsilonRep;
1352  Real sigmaRep;
1353  Real expRep;
1354  Real pow1;
1355  Real pow2;
1356 
1357  BigReal goForce = 0.0;
1358  *goNative = 0;
1359  *goNonnative = 0;
1360 
1361  // determine which Go chain pair we are working with
1362  DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ")" << std::endl);
1363  int goIndex1 = goSigmaIndices[atom1];
1364  int goIndex2 = goSigmaIndices[atom2];
1365 
1366  int32 chain1 = atomChainTypes[goIndex1];
1367  int32 chain2 = atomChainTypes[goIndex2];
1368 
1369  DebugM(3," chain1:" << chain1 << ", chain2:" << chain2 << std::endl);
1370  if (chain1 == 0 || chain2 == 0) return 0.0;
1371 
1372  // retrieve Go cutoff for this chain pair
1373  Real goCutoff = const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2);
1374  DebugM(3," goCutoff:" << goCutoff << std::endl);
1375  if (goCutoff == 0) return 0.0;
1376 
1377  // sigmas are initially set to -1.0 if the atom pair fails go_restricted
1378  // no goSigmas array anymore
1379  //Real sigma_ij = goSigmas[goSigmaIndices[atom1]*numGoAtoms + goSigmaIndices[atom2]];
1380 
1381  // XXX - used to be a condition for the following if
1382  //if the atoms are within 4 of each other
1383  //if ( !atoms_1to4(atom1,atom2) ) {
1384 
1385  // if goSigmaIndices aren't defined, don't calculate forces
1386  if ( goIndex1 != -1 && goIndex2 != -1 ) {
1387  resid1 = goResids[goIndex1];
1388  resid2 = goResids[goIndex2];
1389  residDiff = resid2 - resid1;
1390  if (residDiff < 0) residDiff = -residDiff;
1391  // if this is a Go atom pair that is not restricted
1392  if ( !(const_cast<Molecule*>(this)->go_restricted(chain1,chain2,residDiff)) ) {
1393  xcoorDiff = goCoordinates[goIndex1*3] - goCoordinates[goIndex2*3];
1394  ycoorDiff = goCoordinates[goIndex1*3 + 1] - goCoordinates[goIndex2*3 + 1];
1395  zcoorDiff = goCoordinates[goIndex1*3 + 2] - goCoordinates[goIndex2*3 + 2];
1396  atomAtomDist = sqrt(xcoorDiff*xcoorDiff + ycoorDiff*ycoorDiff + zcoorDiff*zcoorDiff);
1397 
1398  // if attractive then calculate attractive
1399  if ( atomAtomDist <= const_cast<Molecule*>(this)->get_go_cutoff(chain1,chain2) ) {
1400  exp_a = const_cast<Molecule*>(this)->get_go_exp_a(chain1,chain2);
1401  exp_b = const_cast<Molecule*>(this)->get_go_exp_b(chain1,chain2);
1402  sigma_ij = pow(static_cast<double>(exp_b/exp_a),(1.0/(exp_a-exp_b))) * atomAtomDist;
1403 
1404  // [JLai] print out atoms involved in native contacts
1405  // printf("ATOM1: %d C1: %d ATOM2: %d C2: %d\n", atom1,chain1,atom2,chain2);
1406 
1407  epsilon = const_cast<Molecule*>(this)->get_go_epsilon(chain1,chain2);
1408  pow1 = pow(sigma_ij/r,static_cast<double>(exp_a));
1409  pow2 = pow(sigma_ij/r,static_cast<double>(exp_b));
1410  //goForce = ((4/r) * epsilon * (exp_a * pow(sigma_ij/r,exp_a) - exp_b * pow(sigma_ij/r,exp_b)));
1411  goForce = ((4/(r*r)) * epsilon * (exp_a * pow1 - exp_b * pow2));
1412  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);
1413  //goEnergy = (4 * epsilon * ( pow(sigma_ij/r,exp_a) - pow(sigma_ij/r,exp_b) ) ); // JLai I changed some of the expressions
1414  *goNative = (4 * epsilon * ( pow1 - pow2 ) );
1415  *goNonnative = 0;
1416  }
1417 
1418  // if repulsive then calculate repulsive
1419  else {
1420  epsilonRep = const_cast<Molecule*>(this)->get_go_epsilonRep(chain1,chain2);
1421  sigmaRep = const_cast<Molecule*>(this)->get_go_sigmaRep(chain1,chain2);
1422  expRep = const_cast<Molecule*>(this)->get_go_exp_rep(chain1,chain2);
1423  pow1 = pow(sigmaRep/r,(BigReal)expRep);
1424  //goForce = ((12.0/r) * epsilonRep * pow(sigmaRep/r,12.0));
1425  goForce = ((4/(r*r)) * expRep * epsilonRep * pow1);
1426  DebugM(3,"get_go_force - (" << atom1 << "," << atom2 << ") chain1:" << chain1 << ", chain2:" << chain2 << ", epsilonRep:" << epsilonRep << ", sigmaRep:" << sigmaRep << ", r:" << r << ", goForce:" << goForce << std::endl);
1427  //goEnergy = (4 * epsilonRep * pow(sigmaRep/r,12.0)); // JLai I changed some of the expressions
1428  *goNonnative = (4 * epsilonRep * pow1);
1429  *goNative = 0;
1430  }
1431  }
1432  }
1433 
1434  //DebugM(3,"goForce:" << goForce << std::endl);
1435  return goForce;
1436 }
int get_go_exp_a(int chain1, int chain2)
Definition: Molecule.h:1613
Real get_go_cutoff(int chain1, int chain2)
Definition: Molecule.h:1601
short int32
Definition: dumpdcd.c:24
Real * goCoordinates
Definition: Molecule.h:648
Real get_go_epsilonRep(int chain1, int chain2)
Definition: Molecule.h:1604
float Real
Definition: common.h:107
#define DebugM(x, y)
Definition: Debug.h:59
Real get_go_epsilon(int chain1, int chain2)
Definition: Molecule.h:1610
int32 * atomChainTypes
Definition: Molecule.h:643
Real get_go_sigmaRep(int chain1, int chain2)
Definition: Molecule.h:1607
int get_go_exp_rep(int chain1, int chain2)
Definition: Molecule.h:1619
int * goResids
Definition: Molecule.h:649
Bool go_restricted(int, int, int)
Definition: GoMolecule.C:525
BigReal goNonnative
int get_go_exp_b(int chain1, int chain2)
Definition: Molecule.h:1616
BigReal goForce
int32 * goSigmaIndices
Definition: Molecule.h:644
double BigReal
Definition: common.h:112
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().

1607 { return go_array[MAX_GO_CHAINS*chain1 + chain2].sigmaRep; };
#define MAX_GO_CHAINS
Definition: Molecule.h:36
GoValue go_array[MAX_GO_CHAINS *MAX_GO_CHAINS]
Definition: Molecule.h:1584
Real sigmaRep
Definition: Molecule.h:110
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().

1277  {
1278  GridforceGrid *result = NULL;
1279  if (gridnum >= 0 && gridnum < numGridforceGrids) {
1280  result = gridfrcGrid[gridnum];
1281  }
1282  return result;
1283  }
int numGridforceGrids
Definition: Molecule.h:590
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().

1271  {
1272  k = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].k;
1273  q = gridfrcParams[gridnum][gridfrcIndexes[gridnum][atomnum]].q;
1274  }
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.

1180 {
1181  //Initialize return energies to zero
1182  *pairLJEnergy = 0.0;
1183  *pairGaussEnergy = 0.0;
1184 
1185  // Linear search for LJ data
1186  int LJIndex = -1;
1187  int LJbegin = pointerToLJBeg[atom1];
1188  int LJend = pointerToLJEnd[atom1];
1189  for(int i = LJbegin; i <= LJend; i++) {
1190  if(indxLJB[i] == atom2) {
1191  LJIndex = i;
1192  break;
1193  }
1194  }
1195 
1196  // Linear search for Gaussian data
1197  int GaussIndex = -1;
1198  int Gaussbegin = pointerToGaussBeg[atom1];
1199  int Gaussend = pointerToGaussEnd[atom1];
1200  for(int i = Gaussbegin; i <= Gaussend; i++) {
1201  if(indxGaussB[i] == atom2) {
1202  GaussIndex = i;
1203  break;
1204  }
1205  }
1206 
1207  if( LJIndex == -1 && GaussIndex == -1) {
1208  return 0;
1209  } else {
1210  // Code to calculate distances because the pair was found in one of the lists
1211  BigReal r2 = x*x + y*y + z*z;
1212  BigReal r = sqrt(r2);
1213  BigReal ri = 1/r;
1214  BigReal ri6 = (ri*ri*ri*ri*ri*ri);
1215  BigReal ri12 = ri6*ri6;
1216  BigReal ri13 = ri12*ri;
1217  BigReal LJ = 0;
1218  BigReal Gauss = 0;
1219  // Code to calculate LJ
1220  if (LJIndex != -1) {
1221  BigReal ri7 = ri6*ri;
1222  LJ = (12*(pairC12[LJIndex]*ri13) - 6*(pairC6[LJIndex]*ri7));
1223  *pairLJEnergy = (pairC12[LJIndex]*ri12 - pairC6[LJIndex]*ri6);
1224  //std::cout << pairC12[LJIndex] << " " << pairC6[LJIndex] << " " << ri13 << " " << ri7 << " " << LJ << " " << r << "\n";
1225  }
1226  // Code to calculate Gaussian
1227  if (GaussIndex != -1) {
1228  BigReal gr = 12*gRepulsive[GaussIndex]*ri13;
1229  BigReal r1prime = r - gMu1[GaussIndex];
1230  BigReal tmp1 = r1prime * r1prime;
1231  BigReal r2prime = r - gMu2[GaussIndex];
1232  BigReal tmp2 = r2prime * r2prime;
1233  BigReal tmp_gauss1 = 0;
1234  BigReal one_gauss1 = 1;
1235  BigReal tmp_gauss2 = 0;
1236  BigReal one_gauss2 = 1;
1237  if (giSigma1[GaussIndex] != 0) {
1238  tmp_gauss1 = exp(-tmp1*giSigma1[GaussIndex]);
1239  one_gauss1 = 1 - tmp_gauss1;
1240  }
1241  if (giSigma2[GaussIndex] != 0) {
1242  tmp_gauss2 = exp(-tmp2*giSigma2[GaussIndex]);
1243  one_gauss2 = 1 - tmp_gauss2;
1244  }
1245  BigReal A = gA[GaussIndex];
1246  Gauss = gr*one_gauss1*one_gauss2 - 2*A*tmp_gauss1*one_gauss2*r1prime*giSigma1[GaussIndex] \
1247  - 2*tmp_gauss1*one_gauss2*r1prime*giSigma1[GaussIndex]*gRepulsive[GaussIndex]*ri12 - 2*A*tmp_gauss2*one_gauss1*r2prime*giSigma2[GaussIndex] \
1248  - 2*tmp_gauss2*one_gauss1*r2prime*giSigma2[GaussIndex]*gRepulsive[GaussIndex]*ri12;
1249  *pairGaussEnergy = A*(-1+(one_gauss1)*(one_gauss2)*(1+gRepulsive[GaussIndex]*ri12/A));
1250  }
1251  //std::cout << "Net force: " << (LJ + Gauss) << " with ri " << (LJ + Gauss)*ri << "\n";
1252  return (LJ + Gauss)*ri;
1253  }
1254  return 0;
1255 }
Real * gMu2
Definition: Molecule.h:679
Real * giSigma2
Definition: Molecule.h:680
const BigReal A
Real * pairC6
Definition: Molecule.h:667
int * pointerToLJEnd
Definition: Molecule.h:664
Real * gA
Definition: Molecule.h:676
int * indxLJB
Definition: Molecule.h:666
int * pointerToGaussBeg
Definition: Molecule.h:671
Real * gMu1
Definition: Molecule.h:677
int * pointerToLJBeg
Definition: Molecule.h:663
Real * gRepulsive
Definition: Molecule.h:681
int * indxGaussB
Definition: Molecule.h:675
gridSize z
int * pointerToGaussEnd
Definition: Molecule.h:672
Real * pairC12
Definition: Molecule.h:668
Real * giSigma1
Definition: Molecule.h:678
gridSize y
gridSize x
double BigReal
Definition: common.h:112
int Molecule::get_groupSize ( int  )
Improper* Molecule::get_improper ( int  inum) const
inline

Definition at line 1071 of file Molecule.h.

Referenced by dumpbench().

1071 {return (&(impropers[inum]));}
int32* Molecule::get_impropers_for_atom ( int  anum)
inline

Definition at line 1151 of file Molecule.h.

Referenced by dumpbench().

1152  { return impropersByAtom[anum]; }
Lphost* Molecule::get_lphost ( int  atomid) const
inline

Definition at line 1081 of file Molecule.h.

1081  {
1082  // don't call unless simParams->drudeOn == TRUE
1083  // otherwise lphostIndexes array doesn't exist!
1084  int index = lphostIndexes[atomid];
1085  return (index != -1 ? &(lphosts[index]) : NULL);
1086  }
const int32* Molecule::get_mod_exclusions_for_atom ( int  anum) const
inline

Definition at line 1159 of file Molecule.h.

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

1322  {
1323  v = movDragParams[movDragIndexes[atomnum]].v;
1324  }
Bool Molecule::get_noPC ( )
inline

Definition at line 856 of file Molecule.h.

Referenced by ComputeQM::initialize().

856 { return qmNoPC; } ;
int Molecule::get_numQMAtoms ( )
inline
Real* Molecule::get_qmAtmChrg ( )
inline
const int* Molecule::get_qmAtmIndx ( )
inline
const Real* Molecule::get_qmAtomGroup ( ) const
inline
Real Molecule::get_qmAtomGroup ( int  indx) const
inline

Definition at line 820 of file Molecule.h.

820 {return qmAtomGroup[indx]; } ;
Bool Molecule::get_qmcSMD ( )
inline

Definition at line 867 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

867 { return qmcSMD;} ;
int* Molecule::get_qmCustomPCIdxs ( )
inline

Definition at line 865 of file Molecule.h.

Referenced by ComputeQM::initialize().

865 { return qmCustomPCIdxs; } ;
int* Molecule::get_qmCustPCSizes ( )
inline

Definition at line 864 of file Molecule.h.

Referenced by ComputeQM::initialize().

864 { return qmCustPCSizes; } ;
const BigReal* Molecule::get_qmDummyBondVal ( )
inline

Definition at line 834 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

834 { return qmDummyBondVal; } ;
const char* const* Molecule::get_qmDummyElement ( )
inline

Definition at line 839 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

839 { return qmDummyElement; } ;
const char* const* Molecule::get_qmElements ( )
inline

Definition at line 826 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

826 {return qmElementArray;} ;
const int* const* Molecule::get_qmGrpBonds ( )
inline

Definition at line 837 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

837 { return qmGrpBonds; } ;
Real* Molecule::get_qmGrpChrg ( )
inline

Definition at line 830 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

830 {return qmGrpChrg; } ;
Real* Molecule::get_qmGrpID ( )
inline

Definition at line 829 of file Molecule.h.

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

829 { return qmGrpID; } ;
Real* Molecule::get_qmGrpMult ( )
inline

Definition at line 831 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

831 {return qmGrpMult; } ;
const int* Molecule::get_qmGrpNumBonds ( )
inline

Definition at line 835 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

835 { return qmGrpNumBonds; } ;
const int* Molecule::get_qmGrpSizes ( )
inline

Definition at line 828 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

828 {return qmGrpSizes; } ;
int Molecule::get_qmLSSFreq ( )
inline

Definition at line 847 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

847 { return qmLSSFreq; } ;
int* Molecule::get_qmLSSIdxs ( )
inline

Definition at line 845 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

845 { return qmLSSIdxs;} ;
Mass* Molecule::get_qmLSSMass ( )
inline

Definition at line 846 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

846 { return qmLSSMass; } ;
int* Molecule::get_qmLSSRefIDs ( )
inline

Definition at line 849 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

849 { return qmLSSRefIDs ; } ;
Mass* Molecule::get_qmLSSRefMass ( )
inline

Definition at line 851 of file Molecule.h.

851 {return qmLSSRefMass; } ;
int* Molecule::get_qmLSSRefSize ( )
inline

Definition at line 850 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

850 { return qmLSSRefSize ; } ;
int Molecule::get_qmLSSResSize ( )
inline

Definition at line 848 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

848 { return qmLSSResidueSize; } ;
int* Molecule::get_qmLSSSize ( )
inline

Definition at line 844 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

844 { return qmLSSSize;} ;
int* Molecule::get_qmMeMMindx ( )
inline

Definition at line 858 of file Molecule.h.

Referenced by ComputeQM::initialize().

858 { return qmMeMMindx; } ;
int Molecule::get_qmMeNumBonds ( )
inline

Definition at line 857 of file Molecule.h.

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

857 { return qmMeNumBonds; } ;
Real* Molecule::get_qmMeQMGrp ( )
inline

Definition at line 859 of file Molecule.h.

Referenced by ComputeQM::initialize().

859 { return qmMeQMGrp; } ;
const int* const* Molecule::get_qmMMBond ( )
inline

Definition at line 836 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

836 { return qmMMBond; } ;
const int* const* Molecule::get_qmMMBondedIndx ( )
inline

Definition at line 838 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

838 { return qmMMBondedIndx; } ;
const int* const* Molecule::get_qmMMChargeTarget ( )
inline

Definition at line 841 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

841 { return qmMMChargeTarget; } ;
const int* Molecule::get_qmMMNumTargs ( )
inline

Definition at line 842 of file Molecule.h.

Referenced by ComputeQMMgr::recvPntChrg().

842 { return qmMMNumTargs; } ;
std::map<int,int>& Molecule::get_qmMMSolv ( )
inline

Definition at line 852 of file Molecule.h.

852 { return qmClassicSolv;} ;
int Molecule::get_qmNumBonds ( )
inline

Definition at line 833 of file Molecule.h.

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

827 {return qmNumGrps; } ;
int Molecule::get_qmPCFreq ( )
inline

Definition at line 861 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

861 { return qmPCFreq; } ;
Bool Molecule::get_qmReplaceAll ( )
inline

Definition at line 854 of file Molecule.h.

Referenced by ComputeQMMgr::recvPartQM().

854 {return qmReplaceAll; } ;
int Molecule::get_qmTotCustPCs ( )
inline

Definition at line 863 of file Molecule.h.

Referenced by ComputeQM::initialize().

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

145  {
146 
147  if (atomNames == NULL || resLookup == NULL)
148  {
149  NAMD_die("Tried to find atom from name on node other than node 0");
150  }
151  int i = 0;
152  int end = 0;
153  if ( resLookup->lookup(segid,resid,&i,&end) ) return 0;
154  return ( end - i );
155 }
int lookup(const char *segid, int resid, int *begin, int *end) const
Definition: Molecule.C:76
void NAMD_die(const char *err_msg)
Definition: common.C:83
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().

1329  {
1330  v = rotDragParams[rotDragIndexes[atomnum]].v;
1331  a = rotDragParams[rotDragIndexes[atomnum]].a;
1332  p = rotDragParams[rotDragIndexes[atomnum]].p;
1333  }
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().

1352  {
1353  return(ssAtomFlags[anum]);
1354  }
void Molecule::get_stir_refPos ( Vector refPos,
int  atomnum 
) const
inline

Definition at line 1302 of file Molecule.h.

1303  {
1304  refPos = stirParams[stirIndexes[atomnum]].refPos;
1305  }
Real Molecule::get_stir_startTheta ( int  atomnum) const
inline

Definition at line 1314 of file Molecule.h.

Referenced by ComputeStir::doForce().

1315  {
1316  return stirParams[stirIndexes[atomnum]].startTheta;
1317  }
Bond* Molecule::getAllAcceptors ( ) const
inline

Definition at line 1108 of file Molecule.h.

1108 {return acceptors;}
Angle* Molecule::getAllAngles ( ) const
inline

Definition at line 1091 of file Molecule.h.

Referenced by buildAngleData().

1091 {return angles;}
Bond* Molecule::getAllBonds ( ) const
inline

Definition at line 1090 of file Molecule.h.

Referenced by buildBondData().

1090 {return bonds;}
Crossterm* Molecule::getAllCrossterms ( ) const
inline

Definition at line 1094 of file Molecule.h.

Referenced by buildCrosstermData().

1094 {return crossterms;}
Dihedral* Molecule::getAllDihedrals ( ) const
inline

Definition at line 1093 of file Molecule.h.

Referenced by buildDihedralData().

1093 {return dihedrals;}
Bond* Molecule::getAllDonors ( ) const
inline

Definition at line 1107 of file Molecule.h.

1107 {return donors;}
Improper* Molecule::getAllImpropers ( ) const
inline

Definition at line 1092 of file Molecule.h.

Referenced by buildImproperData().

1092 {return impropers;}
Lphost* Molecule::getAllLphosts ( ) const
inline

Definition at line 1098 of file Molecule.h.

1098 { return lphosts; }
BigReal Molecule::GetAtomAlpha ( int  i) const
inline

Definition at line 482 of file Molecule.h.

References drude_constants::alpha.

482  {
483  return drudeConsts[i].alpha;
484  }
AtomNameInfo* Molecule::getAtomNames ( ) const
inline

Definition at line 491 of file Molecule.h.

Referenced by buildAtomData().

491 { return atomNames; }
Atom* Molecule::getAtoms ( ) const
inline

Definition at line 490 of file Molecule.h.

References atoms.

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

490 { return atoms; }
AtomSegResInfo* Molecule::getAtomSegResInfo ( ) const
inline

Definition at line 494 of file Molecule.h.

Referenced by buildAtomData().

494 { return atomSegResids; }
const float* Molecule::getBFactorData ( )
inline

Definition at line 1032 of file Molecule.h.

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

1032 { return (const float *)bfactor; }
BigReal Molecule::getEnergyTailCorr ( const BigReal  ,
const int   
)
int const* Molecule::getLcpoParamType ( )
inline

Definition at line 478 of file Molecule.h.

Referenced by HomePatch::setLcpoType().

478  {
479  return lcpoParamType;
480  }
const float* Molecule::getOccupancyData ( )
inline

Definition at line 1028 of file Molecule.h.

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

1028 { return (const float *)occupancy; }
BigReal Molecule::getVirialTailCorr ( const BigReal  )
Bool Molecule::go_restricted ( int  chain1,
int  chain2,
int  rDiff 
)

Definition at line 525 of file GoMolecule.C.

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

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

526 {
527  int i; // Loop counter
528  for(i=0; i<MAX_RESTRICTIONS; i++) {
529  if (go_array[(MAX_GO_CHAINS*chain1) + chain2].restrictions[i] == rDiff) {
530  return TRUE;
531  } else if (go_array[(MAX_GO_CHAINS*chain1) + chain2].restrictions[i] == -1) {
532  return FALSE;
533  }
534  }
535  return FALSE;
536 }
#define MAX_RESTRICTIONS
Definition: Molecule.h:37
#define FALSE
Definition: common.h:116
#define MAX_GO_CHAINS
Definition: Molecule.h:36
GoValue go_array[MAX_GO_CHAINS *MAX_GO_CHAINS]
Definition: Molecule.h:1584
#define TRUE
Definition: common.h:117
void Molecule::goInit ( )

Definition at line 54 of file GoMolecule.C.

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

54  {
55  numGoAtoms=0;
56  energyNative=0;
58  atomChainTypes=NULL;
59  goSigmaIndices=NULL;
60  goSigmas=NULL;
61  goWithinCutoff=NULL;
62  goCoordinates=NULL;
63  goResids=NULL;
64  goPDB=NULL;
65 }
PDB * goPDB
Definition: Molecule.h:650
Real * goCoordinates
Definition: Molecule.h:648
bool * goWithinCutoff
Definition: Molecule.h:647
int32 * atomChainTypes
Definition: Molecule.h:643
int * goResids
Definition: Molecule.h:649
Real * goSigmas
Definition: Molecule.h:646
BigReal energyNonnative
Definition: Molecule.h:685
BigReal energyNative
Definition: Molecule.h:684
int32 * goSigmaIndices
Definition: Molecule.h:644
int numGoAtoms
Definition: Molecule.h:642
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.

1247  {
1248  if (numConsTorque)
1249  {
1250  // Check the index to see if it is constrained
1251  return(consTorqueIndexes[atomnum] != -1);
1252  }
1253  else
1254  {
1255  // No constraints at all, so just return FALSE
1256  return(FALSE);
1257  }
1258  }
#define FALSE
Definition: common.h:116
int32 * consTorqueIndexes
Definition: Molecule.h:615
int numConsTorque
Definition: Molecule.h:595
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().

1198  {
1199  if (numConstraints)
1200  {
1201  // Check the index to see if it is constrained
1202  return(consIndexes[atomnum] != -1);
1203  }
1204  else
1205  {
1206  // No constraints at all, so just return FALSE
1207  return(FALSE);
1208  }
1209  }
#define FALSE
Definition: common.h:116
int numConstraints
Definition: Molecule.h:588
Bool Molecule::is_atom_exPressure ( int  atomnum) const
inline

Definition at line 1447 of file Molecule.h.

References numExPressureAtoms.

Referenced by Sequencer::langevinPiston().

1448  {
1449  return (numExPressureAtoms && exPressureAtomFlags[atomnum]);
1450  }
int numExPressureAtoms
Definition: Molecule.h:598
Bool Molecule::is_atom_fixed ( int  atomnum) const
inline

Definition at line 1407 of file Molecule.h.

References numFixedAtoms.

Referenced by WorkDistrib::createAtomLists().

1408  {
1409  return (numFixedAtoms && fixedAtomFlags[atomnum]);
1410  }
int numFixedAtoms
Definition: Molecule.h:596
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().

1182  {
1183  if (numGridforceGrids)
1184  {
1185  return(gridfrcIndexes[gridnum][atomnum] != -1);
1186  }
1187  else
1188  {
1189  return(FALSE);
1190  }
1191  }
#define FALSE
Definition: common.h:116
int numGridforceGrids
Definition: Molecule.h:590
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().

1215  {
1216  if (numMovDrag)
1217  {
1218  // Check the index to see if it is constrained
1219  return(movDragIndexes[atomnum] != -1);
1220  }
1221  else
1222  {
1223  // No constraints at all, so just return FALSE
1224  return(FALSE);
1225  }
1226  }
#define FALSE
Definition: common.h:116
int numMovDrag
Definition: Molecule.h:593
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().

1231  {
1232  if (numRotDrag)
1233  {
1234  // Check the index to see if it is constrained
1235  return(rotDragIndexes[atomnum] != -1);
1236  }
1237  else
1238  {
1239  // No constraints at all, so just return FALSE
1240  return(FALSE);
1241  }
1242  }
#define FALSE
Definition: common.h:116
int numRotDrag
Definition: Molecule.h:594
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().

1429  {
1430  if (numStirredAtoms)
1431  {
1432  // Check the index to see if it is constrained
1433  return(stirIndexes[atomnum] != -1);
1434  }
1435  else
1436  {
1437  // No constraints at all, so just return FALSE
1438  return(FALSE);
1439  }
1440  }
#define FALSE
Definition: common.h:116
int numStirredAtoms
Definition: Molecule.h:597
Bool Molecule::is_drude ( int  )
Bool Molecule::is_group_fixed ( int  atomnum) const
inline

Definition at line 1443 of file Molecule.h.

References numFixedAtoms.

1444  {
1445  return (numFixedAtoms && (fixedAtomFlags[atomnum] == -1));
1446  }
int numFixedAtoms
Definition: Molecule.h:596
Bool Molecule::is_hydrogen ( int  )
Bool Molecule::is_hydrogenGroupParent ( int  )
Bool Molecule::is_lp ( int  )
Bool Molecule::is_oxygen ( int  )
Bool Molecule::is_water ( int  )
Real Molecule::langevin_param ( int  atomnum) const
inline

Definition at line 1296 of file Molecule.h.

Referenced by WorkDistrib::createAtomLists().

1297  {
1298  return(langevinParams ? langevinParams[atomnum] : 0.);
1299  }
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().

523  {
524  // local variables prefixed by s_
525  int64_t s_NumDegFreedom = 3 * (int64_t) numAtoms;
526  int64_t s_NumFixedAtoms = num_fixed_atoms();
527  if (s_NumFixedAtoms) s_NumDegFreedom -= 3 * s_NumFixedAtoms;
528  if (numLonepairs) s_NumDegFreedom -= 3 * numLonepairs;
529  if ( ! (s_NumFixedAtoms || numConstraints
530  || simParams->comMove || simParams->langevinOn) ) {
531  s_NumDegFreedom -= 3;
532  }
533  if ( ! isInitialReport && simParams->pairInteractionOn) {
534  //
535  // DJH: a kludge? We want to initially report system degrees of freedom
536  //
537  // this doesn't attempt to deal with fixed atoms or constraints
538  s_NumDegFreedom = 3 * numFepInitial;
539  }
540  int s_NumFixedRigidBonds =
541  (simParams->fixedAtomsOn ? numFixedRigidBonds : 0);
542  if (simParams->watmodel == WAT_TIP4) {
543  // numLonepairs is subtracted here because all lonepairs have a rigid bond
544  // to oxygen, but all of the LP degrees of freedom are dealt with above
545  s_NumDegFreedom -= (numRigidBonds - s_NumFixedRigidBonds - numLonepairs);
546  }
547  else {
548  // Non-polarized systems don't have LPs.
549  // For Drude model, bonds that attach LPs are not counted as rigid;
550  // LPs have already been subtracted from degrees of freedom.
551  s_NumDegFreedom -= (numRigidBonds - s_NumFixedRigidBonds);
552  }
553  return s_NumDegFreedom;
554  }
int num_fixed_atoms() const
Definition: Molecule.h:497
Bool pairInteractionOn
int numLonepairs
Number of lone pairs.
Definition: Molecule.h:580
int numFixedRigidBonds
Definition: Molecule.h:605
int numFepInitial
Definition: Molecule.h:607
int numAtoms
Definition: Molecule.h:556
int numConstraints
Definition: Molecule.h:588
#define WAT_TIP4
Definition: common.h:189
int numRigidBonds
Definition: Molecule.h:604
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().

497  {
498  // local variables prefixed by s_
499  int s_NumFixedAtoms = (simParams->fixedAtomsOn ? numFixedAtoms : 0);
500  return s_NumFixedAtoms; // value is "turned on" SimParameters
501  }
int numFixedAtoms
Definition: Molecule.h:596
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().

503  {
504  // local variables prefixed by s_
505  int s_NumFixedAtoms = num_fixed_atoms();
506  int s_NumFixedGroups = (s_NumFixedAtoms ? numFixedGroups : 0);
507  return s_NumFixedGroups; // value is "turned on" SimParameters
508  }
int numFixedGroups
Definition: Molecule.h:603
int num_fixed_atoms() const
Definition: Molecule.h:497
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().

510  {
511  // local variables prefixed by s_
512  int64_t s_NumGroupDegFreedom = 3 * (int64_t) numHydrogenGroups;
513  int64_t s_NumFixedAtoms = num_fixed_atoms();
514  int s_NumFixedGroups = num_fixed_groups();
515  if (s_NumFixedGroups) s_NumGroupDegFreedom -= 3 * s_NumFixedGroups;
516  if ( ! (s_NumFixedAtoms || numConstraints
517  || simParams->comMove || simParams->langevinOn) ) {
518  s_NumGroupDegFreedom -= 3;
519  }
520  return s_NumGroupDegFreedom;
521  }
int numHydrogenGroups
Definition: Molecule.h:599
int num_fixed_atoms() const
Definition: Molecule.h:497
int numConstraints
Definition: Molecule.h:588
int num_fixed_groups() const
Definition: Molecule.h:503
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().

110  {
111 #ifdef MEM_OPT_VERSION
112  NAMD_die("QMMM interface is not supported in memory optimized builds.");
113 #else
114 
115  PDB *pdbP; // Pointer to PDB object to use
116 
117  qmNumQMAtoms = 0;
118 
119  qmNumGrps = 0 ;
120 
121  iout << iINFO << "Using the following PDB file for QM parameters: " <<
122  pdbFileName << "\n" << endi;
123  if (pdbFileName)
124  pdbP = new PDB(pdbFileName);
125  else
126  iout << "pdbFile name not defined!\n\n" << endi ;
127 
128  if (pdbP->num_atoms() != numAtoms)
129  {
130  NAMD_die("Number of atoms in QM parameter PDB file doesn't match coordinate PDB");
131  }
132 
133  qmElementArray = new char*[numAtoms] ;
134 
135  qmAtomGroup = new Real[numAtoms] ;
136 
137  BigReal bondVal;
138  Real atmGrp;
139 
140  std::set<Real> atmGrpSet;
141 
142  std::vector<Real> grpChrgVec;
143 
144  // Value in the bond column fro QM-MM bonds.
145  std::vector<BigReal> qmBondValVec;
146  // Identifiers of different QM regions
147  std::vector<Real> qmGroupIDsVec;
148  // This maps the qm group ID with the serail index of this group in
149  // various arrays.
150  std::map<Real,int> qmGrpIDMap ;
151 
152  std::vector<int> qmAtmIndxVec, qmGrpSizeVec;
153 
154  // Used to parse user selection in different places.
155  std::vector<std::string> strVec;
156 
157  qmNoPC = simParams->qmNoPC;
158 
159  // We set to zero by default, So there is no need for extra processing.
160  qmPCFreq = 0;
161  if (simParams->qmPCSelFreq > 1)
162  qmPCFreq = simParams->qmPCSelFreq;
163 
164 
167 
168 
169  qmLSSTotalNumAtms = 0;
170  SortedArray< qmSolvData> lssGrpRes;
171  std::map<Real,std::vector<int> > lssGrpRefIDs;
172  refSelStrMap lssRefUsrSel;
173  int totalNumRefAtms = 0;
174  int lssClassicResIndx = 0 ;
175  int lssCurrClassResID = -1 ;
176  char lssCurrClassResSegID[5];
177  if (simParams->qmLSSOn) {
178  DebugM(4, "LSS is ON! Processing QM solvent.\n") ;
179 
180  if (resLookup == NULL)
181  NAMD_die("We need residue data to conduct LSS.");
182 
183  if (simParams->qmLSSMode == QMLSSMODECOM) {
184 
185  StringList *current = cfgList->find("QMLSSRef");
186  for ( ; current; current = current->next ) {
187 
188  strVec = split( std::string(current->data) , " ");
189 
190  if (strVec.size() != 3 ) {
191  iout << iERROR << "Format error in QM LSS size: "
192  << current->data
193  << "\n" << endi;
194  NAMD_die("Error processing QM information.");
195  }
196 
197  std::stringstream storConv ;
198 
199  storConv << strVec[0] ;
200  Real grpID ;
201  storConv >> grpID;
202  if (storConv.fail()) {
203  iout << iERROR << "Error parsing QMLSSRef selection: "
204  << current->data
205  << "\n" << endi;
206  NAMD_die("Error processing QM information.");
207  }
208 
209  std::string segName = strVec[1].substr(0,4);
210 
211  storConv.clear() ;
212  storConv << strVec[2];
213  int resID ;
214  storConv >> resID;
215  if (storConv.fail()) {
216  iout << iERROR << "Error parsing QMLSSRef selection: "
217  << current->data
218  << "\n" << endi;
219  NAMD_die("Error processing QM information.");
220  }
221 
222  auto it = lssRefUsrSel.find(grpID) ;
223  if (it == lssRefUsrSel.end())
224  lssRefUsrSel.insert(refSelStrPair(grpID,refSelStrVec()));
225 
226  lssRefUsrSel[grpID].push_back(refSelStr(segName, resID));
227  }
228 
229  for (auto it=lssRefUsrSel.begin(); it!=lssRefUsrSel.end(); it++) {
230  iout << iINFO << "LSS user selection for COM composition of group "
231  << it->first << ":\n" << endi ;
232  for (int i=0; i<it->second.size();i++) {
233  iout << iINFO << "Segment " << it->second[i].segid
234  << " ; residue " << it->second[i].resid << "\n" << endi ;
235  }
236  }
237  }
238  }
239 
240 
241 
244 
245 
246  for (int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
247 
248  // If we are looking for QM-MM bonds, then we need to store extra info.
249  if (simParams->qmBondOn) {
250 
251  // We store both the qm group and the bond value
252  if ( strcmp(simParams->qmColumn,"beta") == 0 ){
253  atmGrp = pdbP->atom(atmInd)->temperaturefactor() ;
254 
255  bondVal = pdbP->atom(atmInd)->occupancy() ;
256  }
257  else {
258  atmGrp = pdbP->atom(atmInd)->occupancy() ;
259 
260  bondVal = pdbP->atom(atmInd)->temperaturefactor() ;
261  }
262 
263  qmBondValVec.push_back(bondVal);
264  }
265  else {
266 
267  if ( strcmp(simParams->qmColumn,"beta") == 0 ){
268  atmGrp = pdbP->atom(atmInd)->temperaturefactor() ;
269  }
270  else {
271  atmGrp = pdbP->atom(atmInd)->occupancy() ;
272  }
273  }
274 
275  // We store all atom QM Group IDs in an array for
276  // future transport to all PEs.
277  qmAtomGroup[atmInd] = atmGrp;
278 
279  // We store all atom's elements for quick access in the QM code.
280  // Elements are used to tell the QM code the atom type.
281  qmElementArray[atmInd] = new char[3];
282  strncpy(qmElementArray[atmInd],pdbP->atom(atmInd)->element(),3);
283 
284  // For QM atoms, we keep global atom indices and parse different QM Group
285  // IDs, keeping track of each groups size
286  if (atmGrp > 0){
287 
288  if (simParams->fixedAtomsOn) {
289  if ( fixedAtomFlags[atmInd] == 1 ) {
290  iout << iERROR << "QM atom cannot be fixed in space!\n"
291  << endi;
292  NAMD_die("Error processing QM information.");
293  }
294  }
295 
296  // Changes the VdW type of QM atoms.
297  // This is may be used to counter balance the overpolarization
298  // that QM atoms sufer.
299  if (simParams->qmVDW) {
300  // Creating a new type string
301  std::string qmType(qmElementArray[atmInd]) ;
302  // Erases empty spaces
303  qmType.erase(
304  std::remove_if(qmType.begin(), qmType.end(), isspace ),
305  qmType.end());
306  // pre-pends a "q" to the element name.
307  qmType = std::string("q") + qmType;
308 
309 // iout << "QM VdW type: " << atoms[atmInd].vdw_type
310 // << " atom type: " << atomNames[atmInd].atomtype
311 // << " new type "<< qmType << "\n" << endi;
312 
313  /* Look up the vdw constants for this atom */
314  // I am casting a non-const here because the function doesn't actually
315  // change the string value, but it doesn't ask for a const char* as
316  // an argument.
317  params->assign_vdw_index(const_cast<char*>( qmType.c_str()),
318  &(atoms[atmInd]));
319 
320 // iout << "----> new VdW type: " << atoms[atmInd].vdw_type << "\n" << endi;
321  }
322 
323  // Indexes all global indices of QM atoms.
324  qmAtmIndxVec.push_back(atmInd);
325 
326  auto retVal = atmGrpSet.insert(atmGrp);
327 
328  // If a new QM group is found, store the reverse reference in a map
329  // and increase the total count.
330  if (retVal.second) {
331  // This mak makes the reverse identification from the group ID
332  // to the sequential order in which each group was first found.
333  qmGrpIDMap.insert(std::pair<BigReal,int>(atmGrp, qmNumGrps));
334 
335  // This vector keeps track of the group ID for each group
336  qmGroupIDsVec.push_back(atmGrp);
337 
338  // This counter is used to keep track of the sequential order in
339  // which QM groups were first seen in the reference PDB file.
340  qmNumGrps++ ;
341 
342  // If this is a new QM group, initialize a new variable in the
343  // vector to keep track of the number of atoms in this group.
344  qmGrpSizeVec.push_back(1);
345 
346  // For each new QM group, a new entry in the total charge
347  // vector is created
348  grpChrgVec.push_back( atoms[atmInd].charge );
349  }
350  else {
351  // If the QM group has already been seen, increase the count
352  // of atoms in that group.
353  qmGrpSizeVec[ qmGrpIDMap[atmGrp] ]++ ;
354 
355  // If the QM group exists, adds current atom charge to total charge.
356  grpChrgVec[ qmGrpIDMap[atmGrp] ] += atoms[atmInd].charge;
357  }
358 
359  // In case we are using live solvent selection
360  if(simParams->qmLSSOn) {
361 
362  int resID = pdbP->atom(atmInd)->residueseq();
363  char segName[5];
364  strncpy(segName, pdbP->atom(atmInd)->segmentname(),5);
365 
366  int resSize = get_residue_size(segName,resID);
367 
368  int i =-1, end =-1;
369 
370  resLookup->lookup(segName,resID,&i, &end);
371 
372  if ( strncasecmp(pdbP->atom(atmInd)->residuename(),
373  simParams->qmLSSResname, 4) == 0) {
374 
375  // We try to insert every residue from every atom
376  qmSolvData *resP = lssGrpRes.find(qmSolvData(resID, i, resSize));
377 
378  if (resP != NULL) {
379  resP->atmIDs.push_back(atmInd) ;
380 // DebugM(3, "Existing residue " << resP->resID
381 // << " from segID " << resP->segName
382 // << " received atom "
383 // << atmInd << "\n" );
384  }
385  else {
386  lssGrpRes.insert(qmSolvData(resID, i, resSize, segName, atmGrp));
387 // DebugM(3, lssGrpRes.size() << ") Inserted new residue: "
388 // << resID << " from segID " << segName
389 // << " with atom "
390 // << i << "\n" ) ;
391  }
392  }
393  else {
394  // We store the QM atoms, per group, which are NOT part of
395  // solvent molecules.
396 
397  // Checks if we have a vector for this QM group.
398  auto itNewGrp = lssGrpRefIDs.find(atmGrp) ;
399  if (itNewGrp == lssGrpRefIDs.end()) {
400  lssGrpRefIDs.insert(std::pair<Real, std::vector<int> >(
401  atmGrp, std::vector<int>() ) );
402  }
403 
404  switch (simParams->qmLSSMode)
405  {
406 
407  case QMLSSMODECOM:
408  {
409  // If we are using COM selection, checks if the atom
410  // is part of the user selected residues
411  for(int i=0; i<lssRefUsrSel[atmGrp].size(); i++) {
412  if (lssRefUsrSel[atmGrp][i].resid == resID &&
413  strncmp(lssRefUsrSel[atmGrp][i].segid.c_str(),
414  segName,5) == 0 ) {
415 
416  lssGrpRefIDs[atmGrp].push_back(atmInd);
417  totalNumRefAtms++;
418  }
419  }
420 
421  } break;
422 
423  case QMLSSMODEDIST:
424  {
425 
426  lssGrpRefIDs[atmGrp].push_back(atmInd);
427  totalNumRefAtms++;
428 
429  } break;
430 
431  }
432  }
433 
434  }
435 
436  }
437  else if (atmGrp == 0) {
438 
439  if(simParams->qmLSSOn) {
440 
441  if ( strncasecmp(pdbP->atom(atmInd)->residuename(),
442  simParams->qmLSSResname, 4) == 0) {
443 
444  int resID = pdbP->atom(atmInd)->residueseq();
445  char segName[5];
446  strncpy(segName, pdbP->atom(atmInd)->segmentname(),5);
447 
448  if (lssCurrClassResID < 0) {
449  lssCurrClassResID = resID ;
450  strncpy(lssCurrClassResSegID, segName,5);
451  lssClassicResIndx = 0;
452  }
453  else if (lssCurrClassResID != resID ||
454  strcmp(lssCurrClassResSegID, segName) != 0 ) {
455  lssCurrClassResID = resID ;
456  strncpy(lssCurrClassResSegID, segName,5);
457  lssClassicResIndx++;
458  }
459 
460  qmClassicSolv.insert(std::pair<int,int>(atmInd,lssClassicResIndx));
461 
462  }
463  }
464  }
465  else if(atmGrp < 0) {
466  iout << iERROR << "QM group ID cannot be less than zero!\n" << endi;
467  NAMD_die("Error processing QM information.");
468  }
469  }
470 
471  // Sanity check
472  if (simParams->qmLSSOn) {
473  if (lssGrpRes.size() == 0)
474  NAMD_die("LSS was activated but there are no QM solvent molecules!\n");
475 
476  for (auto it=lssRefUsrSel.begin(); it!=lssRefUsrSel.end(); it++) {
477  auto itTarget = qmGrpIDMap.find(it->first);
478  if (itTarget == qmGrpIDMap.end()) {
479  iout << iERROR << "QM group ID for LSS could not be found in input!"
480  << " QM group ID: " << it->first << "\n" << endi;
481  NAMD_die("Error processing QM information.");
482  }
483  }
484 
485  DebugM(3,"We found " << lssClassicResIndx << " classical solvent molecules totalling "
486  << qmClassicSolv.size() << " atoms.\n" );
487 
488  }
489 
490  qmNumQMAtoms = qmAtmIndxVec.size();
491 
492  if (qmNumQMAtoms == 0)
493  NAMD_die("No QM atoms were found in this QM simulation!") ;
494 
495  iout << iINFO << "Number of QM atoms (excluding Dummy atoms): " <<
496  qmNumQMAtoms << "\n" << endi;
497 
498  qmAtmChrg = new Real[qmNumQMAtoms] ;
499  qmAtmIndx = new int[qmNumQMAtoms] ;
500  for (int i=0; i<qmNumQMAtoms; i++) {
501  // qmAtmChrg gets initialized with the PSF charges at the end of this
502  // function, but values may change as QM steps are taken.
503  qmAtmChrg[i] = 0;
504  qmAtmIndx[i] = qmAtmIndxVec[i] ;
505  }
506 
507  // This map relates the QM group index with a vector of pairs
508  // of bonded MM-QM atoms (with the bonded atoms ins this order:
509  // MM first, QM second).
510  std::map<int, std::vector<std::pair<int,int> > > qmGrpIDBonds ;
511  int bondCounter = 0;
512  if (simParams->qmBondOn) {
513 
514  // Checks all bonds for QM-MM pairs.
515  // Makes sanity checks against QM-QM pairs and MM-MM pairs that
516  // are flagged by the user to be bonds between QM and MM regions.
517  // QM-QM bonds will be removed in another function.
518  for (int bndIt = 0; bndIt < numBonds; bndIt++) {
519 
520  bond curr = bonds[bndIt] ;
521 
522  // In case either atom has a non-zero
523  if ( qmBondValVec[curr.atom1] != 0 &&
524  qmBondValVec[curr.atom2] != 0 ) {
525 
526  // Sanity checks (1 of 2)
527  if (qmAtomGroup[curr.atom1] != 0 &&
528  qmAtomGroup[curr.atom2] != 0) {
529  iout << iERROR << "Atoms " << curr.atom1 << " and " <<
530  curr.atom2 << " are assigned as QM atoms.\n" << endi;
531  NAMD_die("Error in QM-MM bond assignment.") ;
532  }
533 
534  // Sanity checks (2 of 2)
535  if (qmAtomGroup[curr.atom1] == 0 &&
536  qmAtomGroup[curr.atom2] == 0) {
537  iout << iERROR << "Atoms " << curr.atom1 << " and " <<
538  curr.atom2 << " are assigned as MM atoms.\n" << endi;
539  NAMD_die("Error in QM-MM bond assignment.") ;
540  }
541 
542  int currGrpID = 0;
543  std::pair<int,int> newPair(0,0);
544 
545  // We create a QM-MM pair with the MM atom first
546  if (qmAtomGroup[curr.atom1] != 0) {
547  newPair = std::pair<int,int>(curr.atom2, curr.atom1) ;
548  currGrpID = qmAtomGroup[curr.atom1];
549  } else {
550  newPair = std::pair<int,int>(curr.atom1, curr.atom2) ;
551  currGrpID = qmAtomGroup[curr.atom2];
552  }
553 
554  int grpIndx = qmGrpIDMap[currGrpID] ;
555 
556  // Stores bonds in vectors key-ed by the QM group ID of the QM atom.
557  auto retIter = qmGrpIDBonds.find(grpIndx) ;
558  // In case thi QM-MM bonds belong to a QM group we have not seen
559  // yet, stores a new vector in the map.
560  if (retIter == qmGrpIDBonds.end()) {
561  qmGrpIDBonds.insert(std::pair<BigReal,std::vector<std::pair<int,int> > >(
562  grpIndx, std::vector<std::pair<int,int> >() ) );
563  }
564 
565  qmGrpIDBonds[grpIndx].push_back( newPair );
566 
567  bondCounter++ ;
568  }
569 
570 
571  }
572 
573 // iout << "Finished processing QM-MM bonds.\n" << endi ;
574 
575  if (bondCounter == 0)
576  iout << iWARN << "We found " << bondCounter << " QM-MM bonds.\n" << endi ;
577  else
578  iout << iINFO << "We found " << bondCounter << " QM-MM bonds.\n" << endi ;
579 
580  }
581 
582  // Initializes several arrays used to setup the QM simulation.
583  qmNumBonds = bondCounter ;
584 
585  qmMMBond = new int*[qmNumBonds];
586 
587  qmDummyBondVal = new BigReal[qmNumBonds];
588 
589  qmMMChargeTarget = new int*[qmNumBonds] ;
590  qmMMNumTargs = new int[qmNumBonds] ;
591 
592  qmDummyElement = new char*[qmNumBonds] ;
593 
594 
595  qmNumGrps = atmGrpSet.size();
596 
597  qmGrpSizes = new int[qmNumGrps];
598 
599  qmGrpID = new Real[qmNumGrps];
600 
601  qmGrpChrg = new Real[qmNumGrps];
602 
603  qmGrpMult = new Real[qmNumGrps] ;
604 
605  qmGrpNumBonds = new int[qmNumGrps];
606 
607  qmGrpBonds = new int*[qmNumGrps];
608  qmMMBondedIndx = new int*[qmNumGrps];
609 
610 
613 
614 
615  // We first initialize the multiplicity vector with
616  // default values, then populate it with user defined values.
617  for (int grpIter = 0; grpIter < qmNumGrps; grpIter++) {
618  qmGrpMult[grpIter] = 0;
619  }
620 
621  int multCount = 0;
622  StringList *current = cfgList->find("QMMult");
623  for ( ; current; current = current->next ) {
624 
625  auto strVec = split( std::string(current->data) , " ");
626 
627  if (strVec.size() != 2 ) {
628  iout << iERROR << "Format error in QM Multiplicity string: "
629  << current->data
630  << "\n" << endi;
631  NAMD_die("Error processing QM information.");
632  }
633 
634  std::stringstream storConv ;
635 
636  storConv << strVec[0] ;
637  Real grpID ;
638  storConv >> grpID;
639  if (storConv.fail()) {
640  iout << iERROR << "Error parsing QMMult selection: "
641  << current->data
642  << "\n" << endi;
643  NAMD_die("Error processing QM information.");
644  }
645 
646  storConv.clear() ;
647  storConv << strVec[1];
648  Real multiplicity ;
649  storConv >> multiplicity;
650  if (storConv.fail()) {
651  iout << iERROR << "Error parsing QMMult selection: "
652  << current->data
653  << "\n" << endi;
654  NAMD_die("Error processing QM information.");
655  }
656 
657  auto it = qmGrpIDMap.find(grpID);
658 
659  if (it == qmGrpIDMap.end()) {
660  iout << iERROR << "Error parsing QMMult selection. Could not find QM group ID: "
661  << grpID
662  << "\n" << endi;
663  NAMD_die("Error processing QM information.");
664  }
665  else {
666  iout << iINFO << "Applying user defined multiplicity "
667  << multiplicity << " to QM group ID " << grpID << "\n" << endi;
668  qmGrpMult[it->second] = multiplicity;
669  }
670 
671  multCount++;
672  }
673 
674  if (multCount != qmNumGrps && simParams->qmFormat == QMFormatORCA) {
675  NAMD_die("ORCA neds multiplicity values for all QM regions.");
676  }
677 
678 
681 
682 
683  for (size_t grpIndx = 0; grpIndx < qmGrpSizeVec.size(); grpIndx++) {
684 
685  bool nonInteger = true;
686  if ((fabsf(roundf(grpChrgVec[grpIndx]) - grpChrgVec[grpIndx]) <= 0.001f) ) {
687  grpChrgVec[grpIndx] = roundf(grpChrgVec[grpIndx]) ;
688  nonInteger = false;
689  }
690 
691  iout << iINFO << grpIndx + 1 << ") Group ID: " << qmGroupIDsVec[grpIndx]
692  << " ; Group size: " << qmGrpSizeVec[grpIndx] << " atoms"
693  << " ; Total PSF charge: " << grpChrgVec[grpIndx] << "\n" << endi ;
694 
695  if (nonInteger && simParams->PMEOn)
696  NAMD_die("QM atoms do not add up to a whole charge, which is needed for PME.") ;
697  }
698 
699  int chrgCount = 0;
700  current = cfgList->find("QMCharge");
701  for ( ; current; current = current->next ) {
702 
703  auto strVec = split( std::string(current->data) , " ");
704 
705  if (strVec.size() != 2 ) {
706  iout << iERROR << "Format error in QM Charge string: "
707  << current->data
708  << "\n" << endi;
709  NAMD_die("Error processing QM information.");
710  }
711 
712  std::stringstream storConv ;
713 
714  storConv << strVec[0] ;
715  Real grpID ;
716  storConv >> grpID;
717  if (storConv.fail()) {
718  iout << iERROR << "Error parsing QMCharge selection: "
719  << current->data
720  << "\n" << endi;
721  NAMD_die("Error processing QM information.");
722  }
723 
724  storConv.clear() ;
725  storConv << strVec[1];
726  Real charge ;
727  storConv >> charge;
728  if (storConv.fail()) {
729  iout << iERROR << "Error parsing QMCharge selection: "
730  << current->data
731  << "\n" << endi;
732  NAMD_die("Error processing QM information.");
733  }
734 
735  auto it = qmGrpIDMap.find(grpID);
736 
737  if (it == qmGrpIDMap.end()) {
738  iout << iERROR << "Error parsing QMCharge selection. Could not find QM group ID: "
739  << grpID
740  << "\n" << endi;
741  NAMD_die("Error processing QM information.");
742  }
743  else {
744  iout << iINFO << "Found user defined charge "
745  << charge << " for QM group ID " << grpID << ". Will ignore PSF charge.\n" << endi;
746  grpChrgVec[it->second] = charge;
747  }
748 
749  chrgCount++;
750  }
751 
752  simParams->qmMOPACAddConfigChrg = false;
753  // Checks if QM group charges can be defined for MOPAC.
754  // Since charge can be defined in QMConfigLine, we need extra logic.
755  if (simParams->qmFormat == QMFormatMOPAC) {
756 
757  // Checks if group charge was supplied in config line for MOPAC.
758  std::string::size_type chrgLoc = std::string::npos ;
759 
760  // This will hold a sting with the first user supplied configuration line for MOPAC.
761  std::string configLine(cfgList->find("QMConfigLine")->data) ;
762  chrgLoc = configLine.find("CHARGE") ;
763 
764  if ( chrgLoc != std::string::npos ) {
765  iout << iINFO << "Found user defined charge in command line. This \
766 will be used for all QM regions and will take precedence over all other charge \
767 definitions.\n" << endi;
768  }
769  else if ( chrgLoc == std::string::npos && (simParams->qmChrgFromPSF || chrgCount == qmNumGrps) ) {
770  // If no charge was defined in the configuration line, gets from PSF and/or
771  // from user defined charges (through the QMCharge keyword).
772  simParams->qmMOPACAddConfigChrg = true;
773  }
774  else
775  {
776  // If we could nont find a charge definition in the config line AND
777  // no specific charge was selected for each QM region through QMCharge AND
778  // the QMChargeFromPSF was not turned ON, the we stop NAMD and scream at the user.
779 // if ( chrgLoc == std::string::npos && (chrgCount != qmNumGrps ) && !simParams->qmChrgFromPSF)
780  iout << iERROR << "Could not find charge for all QM groups. For MOPAC, \
781 charges can be defined through QMConfigLine, QMCharge or QMChargeFromPSF keywords.\n" << endi;
782  NAMD_die("Error processing QM information.");
783  }
784 
785  }
786 
787  if (simParams->qmFormat == QMFormatORCA ) {
788  if ((chrgCount != qmNumGrps ) && !simParams->qmChrgFromPSF) {
789  // If we are not supposed to get charges from the PSF, and not
790  // enough charges were set to cover all QM regions,
791  // we stop NAMD and scream at the user.
792  iout << iERROR << "Could not find charge for all QM groups. For ORCA, \
793 charges can be defined through QMCharge or QMChargeFromPSF keywords.\n" << endi;
794  NAMD_die("Error processing QM information.");
795  }
796  }
797 
798  // If mechanichal embedding was requested but we have QM-MM bonds, we need
799  // to send extra info to ComputeQM to preserve calculation speed.
800  if (qmNumBonds > 0 && qmNoPC) {
801  qmMeNumBonds = qmNumBonds;
802  qmMeMMindx = new int[qmMeNumBonds] ;
803  qmMeQMGrp = new Real[qmMeNumBonds] ;
804  }
805  else {
806  qmMeNumBonds = 0 ;
807  qmMeMMindx = NULL;
808  qmMeQMGrp = NULL;
809  }
810 
811 
814 
815 
816  bondCounter = 0;
817  for (int grpIter = 0; grpIter < qmNumGrps; grpIter++) {
818 
819  qmGrpID[grpIter] = qmGroupIDsVec[grpIter] ;
820 
821  qmGrpSizes[grpIter] = qmGrpSizeVec[grpIter] ;
822 
823  qmGrpChrg[grpIter] = grpChrgVec[grpIter];
824 
825 // iout << "Loaded " << qmGrpSizes[grpIter] << " atoms to this group.\n" << endi;
826 
827  int currNumbBonds = qmGrpIDBonds[grpIter].size() ;
828 
829  // Assigns the number of bonds that the current QM group has.
830  qmGrpNumBonds[grpIter] = currNumbBonds;
831 
832  if (currNumbBonds > 0) {
833 
834  qmGrpBonds[grpIter] = new int[currNumbBonds];
835  qmMMBondedIndx[grpIter] = new int[currNumbBonds];
836 
837  for (int bndIter = 0; bndIter<currNumbBonds; bndIter++) {
838 
839  // Adds the bonds to the overall sequential list.
840  qmMMBond[bondCounter] = new int[2] ;
841  qmMMBond[bondCounter][0] = qmGrpIDBonds[grpIter][bndIter].first ;
842  qmMMBond[bondCounter][1] = qmGrpIDBonds[grpIter][bndIter].second ;
843 
844  // For the current QM group, and the current bond, gets the bond index.
845  qmGrpBonds[grpIter][bndIter] = bondCounter;
846 
847  // For the current QM group, and the current bond, gets the MM atom.
848  qmMMBondedIndx[grpIter][bndIter] = qmGrpIDBonds[grpIter][bndIter].first ;
849 
850  // Assign the default value of dummy element
851  qmDummyElement[bondCounter] = new char[3];
852  strcpy(qmDummyElement[bondCounter],"H\0");
853 
854  // Determines the distance that will separate the new Dummy atom
855  // and the Qm atom to which it will be bound.
856  bondVal = 0;
857  if (simParams->qmBondDist) {
858  if (qmBondValVec[qmMMBond[bondCounter][0]] !=
859  qmBondValVec[qmMMBond[bondCounter][1]]
860  ) {
861  iout << iERROR << "qmBondDist is ON but the values in the bond column are different for atom "
862  << qmMMBond[bondCounter][0] << " and " << qmMMBond[bondCounter][1]
863  << ".\n" << endi ;
864  NAMD_die("Error in QM-MM bond processing.");
865  }
866 
867  bondVal = qmBondValVec[qmMMBond[bondCounter][0]] ;
868  }
869  else {
870 
871  if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"C") == 0 ) {
872  bondVal = 1.09 ;
873  }
874  else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"O") == 0 ) {
875  bondVal = 0.98 ;
876  }
877  else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"N") == 0 ) {
878  bondVal = 0.99 ;
879  }
880  else {
881  iout << iERROR << "qmBondDist is OFF but the QM atom has no default distance value. Atom "
882  << qmMMBond[bondCounter][1] << ", with element: "
883  << qmElementArray[qmMMBond[bondCounter][1]]
884  << ".\n" << endi ;
885  NAMD_die("Error in QM-MM bond processing.");
886  }
887 
888  }
889 
890  iout << iINFO << "MM-QM pair: " << qmMMBond[bondCounter][0] << ":"
891  << qmMMBond[bondCounter][1]
892  << " -> Value (distance or ratio): " << bondVal
893  << " (QM Group " << grpIter << " ID " << qmGrpID[grpIter] << ")"
894  << "\n" << endi ;
895 
896  qmDummyBondVal[bondCounter] = bondVal;
897 
898  // In case we are preparing for a mechanical embedding simulation
899  // with no point charges, populate the following vectors
900  if (qmMeNumBonds > 0) {
901  qmMeMMindx[bondCounter] = qmMMBond[bondCounter][0];
902  qmMeQMGrp[bondCounter] = qmGrpID[grpIter];
903  }
904 
905  bondCounter++ ;
906  }
907  }
908 
909  }
910 
911 
914 
915  current = NULL;
916  if (qmNumBonds > 0)
917  current = cfgList->find("QMLinkElement");
918 
919  int numParsedLinkElem = 0;
920  for ( ; current != NULL; current = current->next ) {
921 
922  DebugM(3,"Parsing link atom element: " << current->data << "\n" );
923 
924  strVec = split( std::string(current->data) , " ");
925 
926  // We need the two atoms that compose the QM-MM bonds and
927  // then the element.
928  if (strVec.size() != 3 && qmNumBonds > 1) {
929  iout << iERROR << "Format error in QM link atom element selection: "
930  << current->data
931  << "\n" << endi;
932  NAMD_die("Error processing QM information.");
933  }
934 
935  // If there is only one QM-MM bond, we can accept the element only.
936  if (strVec.size() != 1 && strVec.size() != 3 && qmNumBonds == 1) {
937  iout << iERROR << "Format error in QM link atom element selection: "
938  << current->data
939  << "\n" << endi;
940  NAMD_die("Error processing QM information.");
941  }
942 
943  std::stringstream storConv ;
944  int bondAtom1, bondAtom2;
945  std::string linkElement ;
946 
947  if (strVec.size() == 1) {
948  linkElement = strVec[0].substr(0,2);
949  }
950  else if (strVec.size() == 3) {
951 
952  storConv << strVec[0] ;
953  storConv >> bondAtom1;
954  if (storConv.fail()) {
955  iout << iERROR << "Error parsing link atom element selection: "
956  << current->data
957  << "\n" << endi;
958  NAMD_die("Error processing QM information.");
959  }
960 
961  storConv.clear() ;
962  storConv << strVec[1];
963  storConv >> bondAtom2;
964  if (storConv.fail()) {
965  iout << iERROR << "Error parsing link atom element selection: "
966  << current->data
967  << "\n" << endi;
968  NAMD_die("Error processing QM information.");
969  }
970 
971  linkElement = strVec[2].substr(0,2);
972  }
973 
974  numParsedLinkElem++;
975 
976  if (numParsedLinkElem>qmNumBonds) {
977  iout << iERROR << "More link atom elements were set than there"
978  " are QM-MM bonds.\n" << endi;
979  NAMD_die("Error processing QM information.");
980  }
981 
982  int bondIter;
983 
984  if (strVec.size() == 1) {
985  bondIter = 0;
986  }
987  else if (strVec.size() == 3) {
988 
989  Bool foundBond = false;
990 
991  for (bondIter=0; bondIter<qmNumBonds; bondIter++) {
992 
993  if ( (qmMMBond[bondIter][0] == bondAtom1 &&
994  qmMMBond[bondIter][1] == bondAtom2 ) ||
995  (qmMMBond[bondIter][0] == bondAtom2 &&
996  qmMMBond[bondIter][1] == bondAtom1 ) ) {
997 
998  foundBond = true;
999  break;
1000  }
1001  }
1002 
1003  if (! foundBond) {
1004  iout << iERROR << "Error parsing link atom element selection: "
1005  << current->data
1006  << "\n" << endi;
1007  iout << iERROR << "Atom pair was not found to be a QM-MM bond.\n"
1008  << endi;
1009  NAMD_die("Error processing QM information.");
1010  }
1011  }
1012 
1013  strcpy(qmDummyElement[bondIter],linkElement.c_str());
1014  qmDummyElement[bondIter][2] = '\0';
1015 
1016  iout << iINFO << "Applying user defined link atom element "
1017  << qmDummyElement[bondIter] << " to QM-MM bond "
1018  << bondIter << ": " << qmMMBond[bondIter][0]
1019  << " - " << qmMMBond[bondIter][1]
1020  << "\n" << endi;
1021  }
1022 
1023 
1024 
1027 
1028 
1029  int32 **bondsWithAtomLocal = NULL ;
1030  int32 *byAtomSizeLocal = NULL;
1031  ObjectArena <int32 >* tmpArenaLocal = NULL;
1032  if (qmNumBonds > 0) {
1033 
1034  bondsWithAtomLocal = new int32 *[numAtoms];
1035  byAtomSizeLocal = new int32[numAtoms];
1036  tmpArenaLocal = new ObjectArena<int32>;
1037 
1038  // Build the bond lists
1039  for (int i=0; i<numAtoms; i++)
1040  {
1041  byAtomSizeLocal[i] = 0;
1042  }
1043  for (int i=0; i<numRealBonds; i++)
1044  {
1045  byAtomSizeLocal[bonds[i].atom1]++;
1046  byAtomSizeLocal[bonds[i].atom2]++;
1047  }
1048  for (int i=0; i<numAtoms; i++)
1049  {
1050  bondsWithAtomLocal[i] = tmpArenaLocal->getNewArray(byAtomSizeLocal[i]+1);
1051  bondsWithAtomLocal[i][byAtomSizeLocal[i]] = -1;
1052  byAtomSizeLocal[i] = 0;
1053  }
1054  for (int i=0; i<numRealBonds; i++)
1055  {
1056  int a1 = bonds[i].atom1;
1057  int a2 = bonds[i].atom2;
1058  bondsWithAtomLocal[a1][byAtomSizeLocal[a1]++] = i;
1059  bondsWithAtomLocal[a2][byAtomSizeLocal[a2]++] = i;
1060  }
1061  }
1062 
1063  // In this loops we try to find other bonds in which the MM atoms from
1064  // QM-MM bonds may be involved. The other MM atoms (which we will call M2 and M3)
1065  // will be involved in charge manipulation. See ComputeQM.C for details.
1066  for (int qmBndIt = 0; qmBndIt < qmNumBonds; qmBndIt++) {
1067 
1068  // The charge targets are accumulated in a temporary vector and then
1069  // transfered to an array that will be transmited to the ComputeQMMgr object.
1070  std::vector<int> chrgTrgt ;
1071 
1072  int MM1 = qmMMBond[qmBndIt][0], MM2, MM2BondIndx, MM3, MM3BondIndx;
1073 
1074  switch (simParams->qmBondScheme) {
1075 
1076  case QMSCHEMERCD:
1077 
1078  case QMSCHEMECS:
1079  {
1080  // Selects ALL MM2 atoms.
1081  for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
1082  MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
1083 
1084  // Checks which of the atoms in the bond structure is the
1085  // MM2 atom.
1086  if (bonds[MM2BondIndx].atom1 == MM1)
1087  MM2 = bonds[MM2BondIndx].atom2;
1088  else
1089  MM2 = bonds[MM2BondIndx].atom1;
1090 
1091  // In case the bonded atom is a QM atom,
1092  // skips the index.
1093  if (qmAtomGroup[MM2] > 0)
1094  continue;
1095 
1096  chrgTrgt.push_back(MM2);
1097  }
1098 
1099  } break;
1100 
1101  case QMSCHEMEZ3:
1102  {
1103  // Selects all MM3 atoms
1104  for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
1105  MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
1106 
1107  // Checks which of the atoms in the bond structure is the
1108  // MM2 atom.
1109  if (bonds[MM2BondIndx].atom1 == MM1)
1110  MM2 = bonds[MM2BondIndx].atom2;
1111  else
1112  MM2 = bonds[MM2BondIndx].atom1;
1113 
1114  // In case the bonded atom is a QM atom,
1115  // skips the index.
1116  if (qmAtomGroup[MM2] > 0)
1117  continue;
1118 
1119  for (int i=0; i<byAtomSizeLocal[MM2]; i++) {
1120  MM3BondIndx = bondsWithAtomLocal[MM2][i] ;
1121 
1122  // Checks which of the atoms in the bond structure is the
1123  // MM3 atom.
1124  if (bonds[MM3BondIndx].atom1 == MM2)
1125  MM3 = bonds[MM3BondIndx].atom2;
1126  else
1127  MM3 = bonds[MM3BondIndx].atom1;
1128 
1129  // In case the bonded atom is a QM atom,
1130  // skips the index.
1131  // We also keep the search from going back to MM1.
1132  if (qmAtomGroup[MM3] > 0 || MM3 == MM1)
1133  continue;
1134 
1135  chrgTrgt.push_back(MM3);
1136  }
1137 
1138  }
1139 
1140  };
1141 
1142  case QMSCHEMEZ2:
1143  {
1144  // Selects all MM2 atoms
1145  for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
1146  MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
1147 
1148  // Checks which of the atoms in the bond structure is the
1149  // MM2 atom.
1150  if (bonds[MM2BondIndx].atom1 == MM1)
1151  MM2 = bonds[MM2BondIndx].atom2;
1152  else
1153  MM2 = bonds[MM2BondIndx].atom1;
1154 
1155  // In case the bonded atom is a QM atom,
1156  // skips the index.
1157  if (qmAtomGroup[MM2] > 0)
1158  continue;
1159 
1160  chrgTrgt.push_back(MM2);
1161  }
1162 
1163  };
1164 
1165  case QMSCHEMEZ1:
1166  {
1167  // Selects all MM1 atoms
1168  chrgTrgt.push_back(MM1);
1169  } break;
1170  }
1171 
1172 
1173  qmMMChargeTarget[qmBndIt] = new int[chrgTrgt.size()] ;
1174  qmMMNumTargs[qmBndIt] = chrgTrgt.size();
1175 
1176  DebugM(3, "MM-QM bond " << qmBndIt << "; MM atom "
1177  << qmMMBond[qmBndIt][0] << " conections: \n" );
1178 
1179  for (size_t i=0; i < chrgTrgt.size(); i++) {
1180  qmMMChargeTarget[qmBndIt][i] = chrgTrgt[i];
1181  DebugM(3,"MM Bonded to: " << chrgTrgt[i] << "\n" );
1182  }
1183 
1184  chrgTrgt.clear();
1185  }
1186 
1187  if (bondsWithAtomLocal != NULL)
1188  delete [] bondsWithAtomLocal; bondsWithAtomLocal = 0;
1189  if (byAtomSizeLocal != NULL)
1190  delete [] byAtomSizeLocal; byAtomSizeLocal = 0;
1191  if (tmpArenaLocal != NULL)
1192  delete tmpArenaLocal; tmpArenaLocal = 0;
1193 
1194 
1197 
1198 
1199  if(simParams->qmLSSOn) {
1200 
1201  std::map<Real,int> grpLSSSize ;
1202  std::map<Real,int>::iterator itGrpSize;
1203 
1204  qmLSSTotalNumAtms = 0;
1205  qmLSSResidueSize = 0;
1206 
1207  if (simParams->qmLSSFreq == 0)
1208  qmLSSFreq = simParams->stepsPerCycle ;
1209  else
1210  qmLSSFreq = simParams->qmLSSFreq;
1211 
1212  #ifdef DEBUG_QM
1213  int resSize = -1;
1214  #endif
1215 
1216  std::map<Real, int> grpNumLSSRes;
1217  std::map<Real, int>::iterator itGrpNumRes;
1218 
1219  for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
1220 
1221  if (it->atmIDs.size() != it->size) {
1222  iout << iERROR << "The number of atoms loaded for residue "
1223  << it->resID << " does not match the expected for this residue type.\n"
1224  << endi;
1225  NAMD_die("Error parsing data for LSS.");
1226  }
1227 
1228  qmLSSTotalNumAtms += it->size;
1229 
1230  #ifdef DEBUG_QM
1231  if (resSize < 0) resSize = it->size ;
1232  if (resSize > 0 and resSize != it->size) {
1233  iout << iERROR << "The number of atoms loaded for residue "
1234  << it->resID << " does not match previously loaded residues.\n"
1235  << endi;
1236  NAMD_die("Error parsing data for LSS.");
1237  }
1238 
1239 // DebugM(3,"Residue " << it->resID << ": " << it->segName
1240 // << " - from " << it->begAtmID << " with size "
1241 // << it->size << " (QM ID: " << it->qmGrpID
1242 // << ") has " << it->atmIDs.size() << " atoms: \n" ) ;
1243 // for (int i=0; i<it->atmIDs.size(); i++)
1244 // DebugM(3, it->atmIDs[i] << "\n" );
1245  #endif
1246 
1247  // Calculating total number of atoms per group
1248  itGrpSize = grpLSSSize.find(it->qmGrpID) ;
1249  if (itGrpSize != grpLSSSize.end())
1250  itGrpSize->second += it->size;