NAMD
NamdState.C
Go to the documentation of this file.
1 /*
2 *** Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
3 *** The Board of Trustees of the University of Illinois.
4 *** All rights reserved.
5 **/
6 
7 /*
8  Holds pointers to large molecule data structure, simulation parameters...
9 */
10 
11 #include "InfoStream.h"
12 #include "common.h"
13 #include "Molecule.h"
14 #include "Parameters.h"
15 #include "SimParameters.h"
16 #include "ConfigList.h"
17 #include "PDB.h"
18 #include "NamdState.h"
19 #include "Controller.h"
20 #include "ScriptTcl.h"
21 #ifndef WIN32
22 #include <unistd.h>
23 #endif
24 #include <sys/stat.h>
25 #include "parm.h"
26 
27 //#define DEBUGM
28 #include "Debug.h"
29 
30 #include "CompressPsf.h"
31 #include "PluginIOMgr.h"
32 #include "BackEnd.h"
33 
35 {
36  configList = NULL;
37  simParameters = NULL;
38  parameters = NULL;
39  molecule = NULL;
40  pdb = NULL;
41 }
42 
43 int
45 {
46  int ret=0;
47  if (configList != NULL) {
48  DebugM(1, "Config List exists\n");
49  } else ret++;
50 
51  if (simParameters != NULL) {
52  DebugM(1, "SimParameters exists\n");
53  }
54  else ret++;
55 
56  if (parameters != NULL) {
57  DebugM(1,"Parameters exists\n");
58  }
59  else ret++;
60 
61  if (molecule != NULL) {
62  DebugM(1, "Molecule exists\n");
63  }
64  else ret++;
65 
66  if (pdb != NULL) {
67  DebugM(1,"PDB exists \n");
68  }
69  else ret++;
70 
71  return(ret);
72 }
73 
75 {
76  controller=controllerPtr;
77 }
78 
80 {
81  controller->run();
82 }
83 
84 extern void read_binary_coors(char *fname, PDB *pdbobj);
85 
86 #ifdef MEM_OPT_VERSION
87  //Check features that are not supported in the memory optimized
88  //version. --Chao Mei
89 void NamdState::checkMemOptCompatibility(){
90  if(simParameters->genCompressedPsf) {
91  iout << "In the memory optimized version, the compression of molecule information is not supported! "
92  << "Please use the non-memory optimized version.\n" <<endi;
93  NAMD_die("MEMOPT: unsupported usage");
94  }
95 
96  if(simParameters->langevinOn && simParameters->langevinDamping == 0.0) {
97  iout << iWARN << "langevinDamping MUST NOT BE 0.0 IF LANGEVIN IS"
98  <<" TURNED ON IN MEMORY OPTIMIZED VERSION\n" <<endi;
99  NAMD_die("MEMOPT: unsupported feature");
100  }
101  if(simParameters->tCoupleOn)
102  NAMD_die("MEMOPT: tCouple is not supported in memory optimized version");
103  if(simParameters->pairInteractionOn)
104  NAMD_die("MEMOPT: pairInteractionOn could not be enabled in memory optimized version");
105  if(simParameters->alchFepOn || simParameters->alchOn){
106  iout << iWARN << "ALCH: AUTOMATIC DELETION OF BONDED INTERACTIONS "
107  << "BETWEEN INITIAL AND FINAL GROUPS IS NOT SUPPORTED IN MEMORY "
108  << "OPTIMISED VERSION - MANUAL PROCESSING IS NECESSARY\n" << endi;
109  NAMD_die("MEMOPT: unsupported feature");
110  }
111  if(simParameters->alchThermIntOn)
112  NAMD_die("MEMOPT: alchThermIntOn could not be enabled in memory optimized version");
113  if(simParameters->lesOn)
114  NAMD_die("MEMOPT: lesOn could not be enabled in memory optimized version");
115  if(simParameters->soluteScalingOn)
116  NAMD_die("MEMOPT: soluteScalingOn could not be enabled in memory optimized version");
117  if(simParameters->lonepairs) {
118  NAMD_die("MEMOPT: lonepairs could not be enabled in memory optimized version");
119  }
120 }
121 #endif
122 
124  configList = cfgList;
125  if (!configList->okay()) {
126  NAMD_die("Simulation config file is incomplete or contains errors.");
127  }
128  DebugM(1,"NamdState::configFileInit configList okay\n");
129 
130  char *currentdir = 0;
131  simParameters = new SimParameters(configList,currentdir);
132  fflush(stdout);
133  lattice = simParameters->lattice;
134 
135  //Check features that are not supported in the memory optimized
136  //version. --Chao Mei
137 #ifdef MEM_OPT_VERSION
138  checkMemOptCompatibility();
139 #endif
140 
141  //Check rigidBonds type when generating the compressed psf files.
142  if(simParameters->genCompressedPsf) {
143  if(simParameters->rigidBonds == RIGID_NONE){
144  //default to RIGID_ALL
145  simParameters->rigidBonds = RIGID_ALL;
146  }
147  }
148 
149  return loadStructure(0,0,0);
150 }
151 
152 int NamdState::loadStructure(const char *molFilename, const char *pdbFilename, int reload) {
153 
154  StringList *molInfoFilename;
155  // If it's AMBER force field, read the AMBER style files;
156  // if it's GROMACS, read the GROMACS files;
157  // Otherwise read the CHARMM style files
158 
159  if (simParameters->amberOn) {
160  if ( reload ) NAMD_die("Molecular structure reloading not supported for Amber input files.\n");
161  StringList *parmFilename = configList->find("parmfile");
162  molInfoFilename = parmFilename;
163  StringList *coorFilename = configList->find("ambercoor");
164  // "amber" is a temporary data structure, which records all
165  // the data from the parm file. After copying them into
166  // molecule, parameter and pdb structures, it will be deleted.
167  Ambertoppar *amber;
168  amber = new Ambertoppar;
169  if (amber->readparm(parmFilename->data))
170  { parameters = new Parameters(amber, simParameters->vdwscale14);
171  molecule = new Molecule(simParameters, parameters, amber);
172  if (coorFilename != NULL)
173  pdb = new PDB(coorFilename->data,amber);
174  delete amber;
175  }
176  else
177  NAMD_die("Failed to read AMBER parm file!");
178  parameters->print_param_summary();
179  }
180  else if (simParameters->gromacsOn) {
181  if ( reload ) NAMD_die("Molecular structure reloading not supported for Gromacs input files.\n");
182  StringList *topFilename = configList->find("grotopfile");
183  molInfoFilename = topFilename;
184  StringList *coorFilename = configList->find("grocoorfile");
185  // "gromacsFile" is a temporary data structure, which records all
186  // the data from the topology file. After copying it into the
187  // molecule and parameter and pdb, it will be deleted.
188  GromacsTopFile *gromacsFile;
189  gromacsFile = new GromacsTopFile(topFilename->data);
190  parameters = new Parameters(gromacsFile,simParameters->minimizeCGOn);
191  if (coorFilename != NULL)
192  pdb = new PDB(coorFilename->data,gromacsFile);
193 
194  molecule = new Molecule(simParameters, parameters, gromacsFile);
195  // XXX does Molecule(needAll,these,arguments)?
196 
197  delete gromacsFile; // XXX unimplemented
198 
199  // XXX add error handling when the file doesn't exist
200  // XXX make sure the right things happen when the parameters are
201  // not even specified.
202  // NAMD_die("Failed to read AMBER parm file!");
203  parameters->print_param_summary();
204  }
205  else if (simParameters->usePluginIO){
206 #ifdef MEM_OPT_VERSION
207  NAMD_die("Using plugin IO is not supported in memory optimized version!");
208 #else
209  if ( pdbFilename ) {
210  NAMD_bug("NamdState::loadStructure pdbFilename non-null with usePluginIO\n");
211  }
212 
213  PluginIOMgr *pIOMgr = new PluginIOMgr();
214 
215  iout << iWARN << "Plugin-based I/O is still in development and may still have bugs\n" << endi;
216 
217  molfile_plugin_t *pIOHandle = pIOMgr->getPlugin();
218  if (pIOHandle == NULL) {
219  NAMD_die("ERROR: Failed to match requested plugin type");
220  }
221  if ( pIOHandle->open_file_read == NULL )
222  NAMD_die("ERROR: Selected plugin type cannot open files");
223  if ( pIOHandle->read_structure == NULL )
224  NAMD_die("ERROR: Selected plugin type cannot read structures");
225  if ( pIOHandle->read_next_timestep == NULL )
226  NAMD_die("ERROR: Selected plugin type cannot read coordinates");
227 
228  StringList *moleculeFilename = configList->find("structure");
229  molInfoFilename = moleculeFilename;
230  if ( ! molFilename ) molFilename = moleculeFilename->data;
231  if ( ! reload ) {
232  StringList *parameterFilename = configList->find("parameters");
233  //****** BEGIN CHARMM/XPLOR type changes
234  // For AMBER use different constructor based on parm_struct!!! -JCP
235  parameters = new Parameters(simParameters, parameterFilename);
236  parameters->print_param_summary();
237  }
238 
239  int numAtoms = 0;
240  //TODO: not sure about the name field in the handler
241  void *plgFile = pIOHandle->open_file_read(molFilename,
242  pIOHandle->name, &numAtoms);
243  if(plgFile == NULL) {
244  NAMD_die("ERROR: Opening structure file failed!");
245  }
246 
247  double fileReadTime = CmiWallTimer();
248  molecule = new Molecule(simParameters, parameters, pIOHandle, plgFile, numAtoms);
249  iout << iINFO << "TIME FOR LOAD MOLECULE STRUCTURE INFORMATION: " << CmiWallTimer() - fileReadTime << "\n" << endi;
250 
251  /* If we are only generating compressed molecule information, the PDB object is not needed */
252  if(!simParameters->genCompressedPsf) {
253  fileReadTime = CmiWallTimer();
254  //get the occupancy data from the Molecule object and then free it
255  //as it is stored in the Molecule object.
256  pdb = new PDB(pIOHandle, plgFile, molecule->numAtoms, molecule->getOccupancyData(), molecule->getBFactorData());
257  molecule->freeOccupancyData();
258  molecule->freeBFactorData();
259  iout << iINFO << "TIME FOR LOADING ATOMS' COORDINATES INFORMATION: " << CmiWallTimer() - fileReadTime << "\n" << endi;
260  }
261 
262  pIOHandle->close_file_read(plgFile);
263  delete pIOMgr;
264 #endif
265  }
266  else
267  {
268  StringList *moleculeFilename = configList->find("structure");
269  molInfoFilename = moleculeFilename;
270  if ( ! molFilename ) molFilename = moleculeFilename->data;
271  if ( ! reload ) {
272  StringList *parameterFilename = configList->find("parameters");
273  //****** BEGIN CHARMM/XPLOR type changes
274  // For AMBER use different constructor based on parm_struct!!! -JCP
275  parameters = new Parameters(simParameters, parameterFilename);
276  //****** END CHARMM/XPLOR type changes
277 
278  parameters->print_param_summary();
279  }
280 
281  double fileReadTime = CmiWallTimer();
282  molecule = new Molecule(simParameters, parameters, (char*)molFilename, configList);
283  iout << iINFO << "TIME FOR READING PSF FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
284 }
285 
286  fflush(stdout);
287 
288 #ifdef MEM_OPT_VERSION
289  //upon knowing the number of atoms, it's good time to estimate the number of
290  //input/output processors if their value is not set
291  if(simParameters->numinputprocs==0){
292  int numatoms = molecule->numAtoms;
293  long estval = (sizeof(InputAtom)+2*sizeof(int)+1)*((long)(numatoms));
294  int numprocs = estval>>26; //considering every input proc consumes about 64M.
295  if(numprocs==0){
296  numprocs=1;
297  }else if(numprocs>CkNumPes()){
298  numprocs=CkNumPes();
299  }
300  simParameters->numinputprocs=numprocs;
301  }
302  if(simParameters->numoutputprocs==0){
303  int numatoms = molecule->numAtoms;
304  long estval = (sizeof(Vector)*2)*((long)(numatoms));
305  int numprocs = estval>>26; //considering every input proc consumes about 64M.
306  if(numprocs==0){
307  numprocs=1;
308  }else if(numprocs>CkNumPes()){
309  numprocs=CkNumPes();
310  }
311  simParameters->numoutputprocs=numprocs;
312  }
313  //check the number of output procs that simultaneously write to a file
314  if(simParameters->numoutputwrts > simParameters->numoutputprocs) {
315  simParameters->numoutputwrts = simParameters->numoutputprocs;
316  }
317 
318  if (simParameters->fixedAtomsOn){
319  double fileReadTime = CmiWallTimer();
320  molecule->load_fixed_atoms(configList->find("fixedAtomListFile"));
321  iout << iINFO << "TIME FOR READING FIXED ATOMS FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
322  }
323 
324  if (simParameters->constraintsOn){
325  double fileReadTime = CmiWallTimer();
326  molecule->load_constrained_atoms(configList->find("consAtomListFile"));
327  iout << iINFO << "TIME FOR READING CONSTRAINED ATOMS FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
328  }
329 #else
330  if (simParameters->extraBondsOn) {
331  //The extra bonds building will be executed in read_compressed_psf in
332  //the memory optimized version, so avoid calling this function in the
333  //memory optimized run.
334  if(!simParameters->useCompressedPsf)
335  molecule->build_extra_bonds(parameters, configList->find("extraBondsFile"));
336  }
337  if(simParameters->genCompressedPsf) {
338  double fileReadTime = CmiWallTimer();
339  compress_molecule_info(molecule, molInfoFilename->data, parameters, simParameters, configList);
340  iout << "Finished compressing molecule information, which takes " << CmiWallTimer()-fileReadTime <<"(s)\n"<<endi;
341  BackEnd::exit();
342  }
343 
344  //If using plugin-based IO, the PDB object is already created!
345  StringList *coordinateFilename = NULL;
346  if(!simParameters->usePluginIO) {
347  //In the memory opt version, the coordinates of atoms
348  //are read during startup in parallel with a bincoordinates input
349  //-Chao Mei
350  double fileReadTime = CmiWallTimer();
351  if ( pdbFilename ) {
352  iout << iINFO << "Reading pdb file " << pdbFilename << "\n" << endi;
353  pdb = new PDB(pdbFilename);
354  } else {
355  coordinateFilename = configList->find("coordinates");
356  if (coordinateFilename != NULL) {
357  iout << iINFO << "Reading pdb file " << coordinateFilename->data << "\n" << endi;
358  pdb = new PDB(coordinateFilename->data);
359  }
360  }
361  if (pdb->num_atoms() != molecule->numAtoms) {
362  NAMD_die("Number of pdb and psf atoms are not the same!");
363  }
364  iout << iINFO << "TIME FOR READING PDB FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
365  iout << iINFO << "\n" << endi;
366  }
367 
368  // If constraints are active, build the parameters necessary
369  if (simParameters->constraintsOn)
370  {
371  StringList *consRefFile = configList->find("consref");
372  StringList *consKFile = configList->find("conskfile");
373 
374  if (coordinateFilename != NULL) {
375  if(strcasecmp(coordinateFilename->data, consRefFile->data)==0)
376  consRefFile = NULL;
377  if(strcasecmp(coordinateFilename->data, consKFile->data)==0)
378  consKFile = NULL;
379  }
380 
381  molecule->build_constraint_params(consRefFile, consKFile,
382  configList->find("conskcol"),
383  pdb,
384  NULL);
385  }
386 #endif
387  //CkPrintf ("DEBUG--check if StirOn to build stir params..\n");
388 
389  if (simParameters->stirOn)
390  {
391  //CkPrintf ("DEBUG--now to build stir params..\n");
392 
393  molecule->build_stirred_atoms(configList->find("stirFilename"),
394  configList->find("stirredAtomsCol"),
395  pdb,
396  NULL);
397  }
398 
399 
400 #ifndef MEM_OPT_VERSION
401  if (simParameters->fixedAtomsOn)
402  {
403  molecule->build_fixed_atoms(configList->find("fixedatomsfile"),
404  configList->find("fixedatomscol"),
405  pdb,
406  NULL);
407  }
408 #endif
409 
410  /* BEGIN gf */
411  if (simParameters->mgridforceOn)
412  {
413  molecule->build_gridforce_params(configList->find("gridforcefile"),
414  configList->find("gridforcecol"),
415  configList->find("gridforcechargecol"),
416  configList->find("gridforcepotfile"),
417  pdb,
418  NULL);
419  }
420  /* END gf */
421 
422  // If constant forces are active, read the forces necessary
423  if (simParameters->consForceOn) {
424  char *filename = NULL;
425  if (configList->find("consforcefile"))
426  filename = configList->find("consforcefile")->data;
427  molecule->build_constant_forces(filename);
428  }
429 
430  if (simParameters->excludeFromPressure) {
431  molecule->build_exPressure_atoms(
432  configList->find("excludeFromPressureFile"),
433  configList->find("excludeFromPressureCol"),
434  pdb, NULL);
435  }
436 
437  // If moving drag is active, build the parameters necessary
438  if (simParameters->movDragOn) {
439  molecule->build_movdrag_params(configList->find("movDragFile"),
440  configList->find("movDragCol"),
441  configList->find("movDragVelFile"),
442  pdb,
443  NULL);
444  }
445 
446  // If rotating drag is active, build the parameters necessary
447  if (simParameters->rotDragOn) {
448  molecule->build_rotdrag_params(configList->find("rotDragFile"),
449  configList->find("rotDragCol"),
450  configList->find("rotDragAxisFile"),
451  configList->find("rotDragPivotFile"),
452  configList->find("rotDragVelFile"),
453  configList->find("rotDragVelCol"),
454  pdb,
455  NULL);
456  }
457 
458  // If "constant" torque is active, build the parameters necessary
459  if (simParameters->consTorqueOn) {
460  molecule->build_constorque_params(configList->find("consTorqueFile"),
461  configList->find("consTorqueCol"),
462  configList->find("consTorqueAxisFile"),
463  configList->find("consTorquePivotFile"),
464  configList->find("consTorqueValFile"),
465  configList->find("consTorqueValCol"),
466  pdb,
467  NULL);
468  }
469 
470 #ifndef MEM_OPT_VERSION
471  // If langevin dynamics or temperature coupling are active, build
472  // the parameters necessary
473  if (simParameters->langevinOn)
474  {
475  if (simParameters->langevinDamping == 0.0) {
476  molecule->build_langevin_params(configList->find("langevinfile"),
477  configList->find("langevincol"),
478  pdb,
479  NULL);
480  } else {
481  molecule->build_langevin_params(simParameters->langevinDamping,
482  simParameters->drudeDamping,
483  simParameters->langevinHydrogen);
484  }
485  }
486  else if (simParameters->tCoupleOn)
487  {
488  // Temperature coupling uses the same parameters, but with different
489  // names . . .
490  molecule->build_langevin_params(configList->find("tcouplefile"),
491  configList->find("tcouplecol"),
492  pdb,
493  NULL);
494  }
495 
496  // Modifications for alchemical fep
497  // identify the mutant atoms for fep simulation
498  if (simParameters->alchOn) {
499  molecule->build_fep_flags(configList->find("alchfile"),
500  configList->find("alchcol"), pdb, NULL, "alch" );
501  molecule->delete_alch_bonded();
502  if (simParameters->sdScaling) {
503  if (configList->find("unperturbedBondFile") == NULL) {
504  NAMD_die("Input file for Shobana's bond terms is required with sdScaling on");
505  }
506  molecule->build_alch_unpert_bond_lists(configList->find("unperturbedBondFile")->data);
507  }
508  }
509 //fepe
510 
511  if (simParameters->lesOn) {
512  if (simParameters->alchOn) NAMD_bug("FEP/TI and LES are incompatible!");
513  molecule->build_fep_flags(configList->find("lesfile"),
514  configList->find("lescol"), pdb, NULL, "les");
515  }
516  if (simParameters->soluteScalingOn) {
517  molecule->build_ss_flags(configList->find("soluteScalingFile"),
518  configList->find("soluteScalingCol"), pdb, NULL);
519  }
520  if (simParameters->pairInteractionOn) {
521  molecule->build_fep_flags(configList->find("pairInteractionFile"),
522  configList->find("pairInteractionCol"), pdb, NULL, "pairInteraction");
523  }
524  if (simParameters->pressureProfileAtomTypes > 1) {
525  molecule->build_fep_flags(configList->find("pressureProfileAtomTypesFile"),
526  configList->find("pressureProfileAtomTypesCol"), pdb, NULL,
527  "pressureProfileAtomTypes");
528  }
529 
530  #ifdef OPENATOM_VERSION
531  if (simParameters->openatomOn) {
532  molecules->build_qmmm_flags(configList->find("openatomPdbFile",
533  configList->find("openatomPdbCol"), pdb, NULL, "openatomPdb")
534  }
535  #endif // OPENATOM_VERSION
536 
537  if (simParameters->qmForcesOn){
538 
539 #ifdef MEM_OPT_VERSION
540  NAMD_die("QM forces are not supported in memory-optimized builds.");
541 #endif
542 
543 #ifdef NAMD_CUDA
544  NAMD_die("QM forces are not compatible with CUDA at this time");
545 #endif
546 
547  molecule->set_qm_replaceAll(simParameters->qmReplaceAll);
548 
549  if (simParameters->qmParamPDBDefined)
550  molecule->prepare_qm(simParameters->qmParamPDB,
551  parameters, configList);
552  else if (pdbFilename)
553  molecule->prepare_qm(pdbFilename,
554  parameters, configList);
555  else
556  molecule->prepare_qm(configList->find("coordinates")->data,
557  parameters, configList);
558 
559  }
560 
561 
562  if (simParameters->LJcorrection) {
563  molecule->compute_LJcorrection();
564  }
565 #endif
566 
567  // JLai checks to see if Go Forces are turned on
568  if (simParameters->goForcesOn) {
569 #ifdef MEM_OPT_VERSION
570  NAMD_die("Go forces are not supported in memory-optimized builds.");
571 #else
572  StringList *moleculeFilename = configList->find("structure");
573  StringList *parameterFilename = configList->find("parameters");
574  StringList *goFilename = configList->find("goParameters");
575  StringList *goStructureFilename = configList->find("goCoordinates");
576 
577  // Added by JLai -- 1.10.12 -- Code to build the Go parameters (from within the Go molecule instead of the parameters object)
578  molecule->build_go_params(goFilename);
579  // Added by JLai -- 6.3.11 -- flag to switch between the different goMethodologies
580  int goMethod = simParameters->goMethod;
581  if (goMethod == 1) { // should probably replace with switch statement
582  iout << iINFO << "Using Go method matrix\n" << endi;
583  molecule->build_go_sigmas(goStructureFilename, NULL);
584  } else if (goMethod == 2) {
585  iout << iINFO << "Using Go method jump table\n" << endi;
586  molecule->build_go_sigmas2(goStructureFilename, NULL);
587  } else if (goMethod == 3) {
588  iout << iINFO << "Using Go method lowmem\n" << endi;
589  molecule->build_go_arrays(goStructureFilename, NULL);
590  } else {
591  NAMD_die("Failed to read goMethod variable in NamdState.C");
592  }
593 #endif
594  }
595  // End of Go code -- JLai
596 
597 #ifndef MEM_OPT_VERSION
598  iout << iINFO << "****************************\n";
599  iout << iINFO << "STRUCTURE SUMMARY:\n";
600  iout << iINFO << molecule->numAtoms << " ATOMS\n";
601  iout << iINFO << molecule->numBonds << " BONDS\n";
602  iout << iINFO << molecule->numAngles << " ANGLES\n";
603  iout << iINFO << molecule->numDihedrals << " DIHEDRALS\n";
604  iout << iINFO << molecule->numImpropers << " IMPROPERS\n";
605  iout << iINFO << molecule->numCrossterms << " CROSSTERMS\n";
606  iout << iINFO << molecule->numExclusions << " EXCLUSIONS\n";
607 
608  //****** BEGIN CHARMM/XPLOR type changes
609  if ((molecule->numMultipleDihedrals) && (simParameters->paraTypeXplorOn))
610  {
611  iout << iINFO << molecule->numMultipleDihedrals
612  << " DIHEDRALS WITH MULTIPLE PERIODICITY (BASED ON PSF FILE)\n";
613  }
614  if ((molecule->numMultipleDihedrals) && (simParameters->paraTypeCharmmOn))
615  {
616  iout << iINFO << molecule->numMultipleDihedrals
617  << " DIHEDRALS WITH MULTIPLE PERIODICITY IGNORED (BASED ON PSF FILE) \n";
618  iout << iINFO
619  << " CHARMM MULTIPLICITIES BASED ON PARAMETER FILE INFO! \n";
620  }
621  //****** END CHARMM/XPLOR type changes
622 
623  if (molecule->numMultipleImpropers)
624  {
625  iout << iINFO << molecule->numMultipleImpropers
626  << " IMPROPERS WITH MULTIPLE PERIODICITY\n";
627  }
628 
629  if (simParameters->constraintsOn)
630  {
631  iout << iINFO << molecule->numConstraints << " CONSTRAINTS\n";
632  }
633 
634  if (simParameters->consForceOn)
635  iout << iINFO << molecule->numConsForce << " CONSTANT FORCES\n";
636 
637  if (simParameters->stirOn)
638  iout << iINFO << molecule->numStirredAtoms << " STIRRED ATOMS\n";
639 
640  if (simParameters->fixedAtomsOn)
641  {
642  iout << iINFO << molecule->numFixedAtoms << " FIXED ATOMS\n";
643  }
644 
645  if (simParameters->rigidBonds)
646  {
647  iout << iINFO << molecule->numRigidBonds << " RIGID BONDS\n";
648  }
649 
650  if (simParameters->fixedAtomsOn && simParameters->rigidBonds)
651  {
652  iout << iINFO << molecule->numFixedRigidBonds <<
653  " RIGID BONDS BETWEEN FIXED ATOMS\n";
654  }
655 
656  /* BEGIN gf */
657  if (simParameters->mgridforceOn)
658  {
659  int i;
660  iout << iINFO << molecule->numGridforceGrids
661  << " GRIDS ACTIVE\n";
662  }
663  /* END gf */
664 
665 //Modifications for alchemical fep
666  if (simParameters->alchOn) {
667  iout << iINFO << "ALCH: "
668  << molecule->numFepInitial <<
669  " ATOMS TO DISAPPEAR IN FINAL STATE\n";
670  iout << iINFO << "ALCH: "
671  << molecule->numFepFinal <<
672  " ATOMS TO APPEAR IN FINAL STATE\n";
673  if (molecule->suspiciousAlchBonds) {
674  iout << iWARN << "ALCH: SUSPICIOUS BONDS BETWEEN INITIAL AND " <<
675  "FINAL GROUPS WERE FOUND" << "\n" << endi;
676  }
677  if (molecule->alchDroppedAngles) {
678  iout << iINFO << "ALCH: "
679  << molecule->alchDroppedAngles <<
680  " ANGLES LINKING INITIAL AND FINAL ATOMS DELETED\n";
681  }
682  if (molecule->alchDroppedDihedrals) {
683  iout << iINFO << "ALCH: "
684  << molecule->alchDroppedDihedrals <<
685  " DIHEDRALS LINKING INITIAL AND FINAL ATOMS DELETED\n";
686  }
687  if (molecule->alchDroppedImpropers) {
688  iout << iINFO << "ALCH: "
689  << molecule->alchDroppedImpropers <<
690  " IMPROPERS LINKING INITIAL AND FINAL ATOMS DELETED\n";
691  }
692  }
693 //fepe
694 
695  if (simParameters->lesOn) {
696  iout << iINFO << molecule->numFepInitial <<
697  " LOCALLY ENHANCED ATOMS ENABLED\n";
698  }
699 
700  if (simParameters->soluteScalingOn) {
701  iout << iINFO << " SOLUTE SCALING ENABLED\n";
702  }
703 
704  if (simParameters->pairInteractionOn) {
705  iout << iINFO << "PAIR INTERACTION GROUP 1 CONTAINS "
706  << molecule->numFepInitial << " ATOMS\n";
707  if (!simParameters->pairInteractionSelf) {
708  iout << iINFO << "PAIR INTERACTION GROUP 2 CONTAINS "
709  << molecule->numFepFinal << " ATOMS\n";
710  }
711  }
712 
713 #if 1
714  if (molecule->numLonepairs != 0) {
715  iout << iINFO << molecule->numLonepairs << " LONE PAIRS\n";
716  }
717  if (molecule->numDrudeAtoms != 0) {
718  iout << iINFO << molecule->numDrudeAtoms << " DRUDE ATOMS\n";
719  }
720  iout << iINFO << molecule->num_deg_freedom(1)
721  << " DEGREES OF FREEDOM\n";
722  if (simParameters->drudeOn) {
723  int g_bond = 3 * molecule->numDrudeAtoms;
724  int g_com = molecule->num_deg_freedom(1) - g_bond;
725  iout << iINFO << g_com << " DRUDE COM DEGREES OF FREEDOM\n";
726  iout << iINFO << g_bond << " DRUDE BOND DEGREES OF FREEDOM\n";
727  }
728 #endif
729 #if 0
730  {
731  // Copied from Controller::printEnergies()
732  int64 numAtoms = molecule->numAtoms;
733  int64 numDegFreedom = 3 * numAtoms;
734  int numLonepairs = molecule->numLonepairs;
735  int numFixedAtoms = molecule->numFixedAtoms;
736  if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
737  if ( ! ( numFixedAtoms || molecule->numConstraints
738  || simParameters->comMove || simParameters->langevinOn ) ) {
739  numDegFreedom -= 3;
740  }
741  if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
742  int numRigidBonds = molecule->numRigidBonds;
743  int numFixedRigidBonds = molecule->numFixedRigidBonds;
744  // numLonepairs is subtracted here because all lonepairs have a rigid bond
745  // to oxygen, but all of the LP degrees of freedom are dealt with above
746  numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
747  iout << iINFO << numDegFreedom << " DEGREES OF FREEDOM\n";
748  }
749 #endif
750 
751  iout << iINFO << molecule->numHydrogenGroups << " HYDROGEN GROUPS\n";
752  iout << iINFO << molecule->maxHydrogenGroupSize
753  << " ATOMS IN LARGEST HYDROGEN GROUP\n";
754  iout << iINFO << molecule->numMigrationGroups << " MIGRATION GROUPS\n";
755  iout << iINFO << molecule->maxMigrationGroupSize
756  << " ATOMS IN LARGEST MIGRATION GROUP\n";
757  if (simParameters->fixedAtomsOn)
758  {
759  iout << iINFO << molecule->numFixedGroups <<
760  " HYDROGEN GROUPS WITH ALL ATOMS FIXED\n";
761  }
762 
763  {
764  BigReal totalMass = 0;
765  BigReal totalCharge = 0;
766  int i;
767  for ( i = 0; i < molecule->numAtoms; ++i ) {
768  totalMass += molecule->atommass(i);
769  totalCharge += molecule->atomcharge(i);
770  }
771  iout << iINFO << "TOTAL MASS = " << totalMass << " amu\n";
772  iout << iINFO << "TOTAL CHARGE = " << totalCharge << " e\n";
773 
774  BigReal volume = lattice.volume();
775  if ( volume ) {
776  iout << iINFO << "MASS DENSITY = "
777  << ((totalMass/volume) / 0.6022) << " g/cm^3\n";
778  iout << iINFO << "ATOM DENSITY = "
779  << (molecule->numAtoms/volume) << " atoms/A^3\n";
780  }
781  }
782 
783  iout << iINFO << "*****************************\n";
784  iout << endi;
785  fflush(stdout);
786 
787  StringList *binCoordinateFilename = configList->find("bincoordinates");
788  if ( binCoordinateFilename && ! reload ) {
789  read_binary_coors(binCoordinateFilename->data, pdb);
790  }
791 
792  DebugM(4, "::configFileInit() - printing Molecule Information\n");
793 
794  molecule->print_atoms(parameters);
795  molecule->print_bonds(parameters);
796  molecule->print_exclusions();
797  fflush(stdout);
798 #endif
799 
800  DebugM(4, "::configFileInit() - done printing Molecule Information\n");
801  DebugM(1, "::configFileInit() - done\n");
802 
803  return(0);
804 }
805 
int numFixedGroups
Definition: Molecule.h:603
void build_gridforce_params(StringList *, StringList *, StringList *, StringList *, PDB *, char *)
Definition: Molecule.C:6373
int suspiciousAlchBonds
Definition: Molecule.h:562
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
Definition: PDB.h:35
int alchDroppedDihedrals
Definition: Molecule.h:564
int numBonds
Definition: Molecule.h:559
void print_bonds(Parameters *)
Definition: Molecule.C:5380
int status()
Definition: NamdState.C:44
void build_extra_bonds(Parameters *parameters, StringList *file)
void print_exclusions()
Definition: Molecule.C:5417
int numHydrogenGroups
Definition: Molecule.h:599
molfile_plugin_t * getPlugin()
Definition: PluginIOMgr.h:22
static void exit(int status=0)
Definition: BackEnd.C:276
Definition: Vector.h:64
Bool excludeFromPressure
static int numatoms
Definition: ScriptTcl.C:64
int alchDroppedImpropers
Definition: Molecule.h:565
#define DebugM(x, y)
Definition: Debug.h:59
int numGridforceGrids
Definition: Molecule.h:590
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
void freeOccupancyData()
Definition: Molecule.h:1030
#define iout
Definition: InfoStream.h:87
int num_atoms(void)
Definition: PDB.C:323
int loadStructure(const char *, const char *, int)
Definition: NamdState.C:152
int readparm(char *)
Definition: parm.C:151
Bool pairInteractionOn
void build_go_arrays(StringList *, char *)
Definition: GoMolecule.C:950
int64_t num_deg_freedom(int isInitialReport=0) const
Definition: Molecule.h:523
int numFepFinal
Definition: Molecule.h:608
BigReal vdwscale14
void print_param_summary()
Definition: Parameters.C:4966
void runController(void)
Definition: NamdState.C:79
void build_langevin_params(BigReal coupling, BigReal drudeCoupling, Bool doHydrogen)
NamdState(void)
Definition: NamdState.C:34
int numMultipleImpropers
Definition: Molecule.h:637
void read_binary_coors(char *fname, PDB *pdbobj)
Definition: NamdOneTools.C:34
struct parm Ambertoppar
int numMultipleDihedrals
Definition: Molecule.h:635
int numLonepairs
Number of lone pairs.
Definition: Molecule.h:580
const float * getBFactorData()
Definition: Molecule.h:1032
void NAMD_bug(const char *err_msg)
Definition: common.C:123
int numFixedRigidBonds
Definition: Molecule.h:605
BigReal langevinDamping
void build_constraint_params(StringList *, StringList *, StringList *, PDB *, char *)
GoChoices goMethod
void build_constorque_params(StringList *, StringList *, StringList *, StringList *, StringList *, StringList *, PDB *, char *)
int numFepInitial
Definition: Molecule.h:607
void set_qm_replaceAll(Bool newReplaceAll)
Definition: Molecule.h:817
int numFixedAtoms
Definition: Molecule.h:596
void prepare_qm(const char *pdbFileName, Parameters *params, ConfigList *cfgList)
Definition: MoleculeQM.C:109
int numAngles
Definition: Molecule.h:560
void build_stirred_atoms(StringList *, StringList *, PDB *, char *)
const float * getOccupancyData()
Definition: Molecule.h:1028
Definition: parm.h:15
int numAtoms
Definition: Molecule.h:556
int numStirredAtoms
Definition: Molecule.h:597
int numCrossterms
Definition: Molecule.h:567
void NAMD_die(const char *err_msg)
Definition: common.C:83
void run(void)
Definition: Controller.C:268
int numConsForce
Definition: Molecule.h:611
void build_ss_flags(const StringList *ssfile, const StringList *sscol, PDB *initial_pdb, const char *cwd)
void build_rotdrag_params(StringList *, StringList *, StringList *, StringList *, StringList *, StringList *, PDB *, char *)
BigReal volume(void) const
Definition: Lattice.h:277
int maxMigrationGroupSize
Definition: Molecule.h:602
int numDihedrals
Definition: Molecule.h:561
void build_fep_flags(StringList *, StringList *, PDB *, char *, const char *)
int numImpropers
Definition: Molecule.h:566
int numConstraints
Definition: Molecule.h:588
void build_alch_unpert_bond_lists(char *)
void compress_molecule_info(Molecule *mol, char *psfFileName, Parameters *param, SimParameters *simParam, ConfigList *cfgList)
Definition: CompressPsf.C:436
void build_fixed_atoms(StringList *, StringList *, PDB *, char *)
void compute_LJcorrection()
long long int64
Definition: common.h:34
int pressureProfileAtomTypes
void print_atoms(Parameters *)
Definition: Molecule.C:5337
int maxHydrogenGroupSize
Definition: Molecule.h:600
char * data
Definition: ConfigList.h:48
int numMigrationGroups
Definition: Molecule.h:601
void useController(Controller *controllerPtr)
Definition: NamdState.C:74
void build_go_sigmas(StringList *, char *)
Definition: GoMolecule.C:577
Real atomcharge(int anum) const
Definition: Molecule.h:1048
void delete_alch_bonded(void)
StringList * find(const char *name) const
Definition: ConfigList.C:341
#define RIGID_ALL
Definition: SimParameters.h:78
void build_exPressure_atoms(StringList *, StringList *, PDB *, char *)
Bool pairInteractionSelf
char qmParamPDB[128]
void build_constant_forces(char *)
Bool qmParamPDBDefined
Real atommass(int anum) const
Definition: Molecule.h:1038
int numDrudeAtoms
Number of Drude particles.
Definition: Molecule.h:581
infostream & endi(infostream &s)
Definition: InfoStream.C:38
Bool okay(void)
Definition: ConfigList.h:120
int alchDroppedAngles
Definition: Molecule.h:563
void build_movdrag_params(StringList *, StringList *, StringList *, PDB *, char *)
#define RIGID_NONE
Definition: SimParameters.h:77
BigReal drudeDamping
int numRigidBonds
Definition: Molecule.h:604
int configListInit(ConfigList *)
Definition: NamdState.C:123
void freeBFactorData()
Definition: Molecule.h:1034
double BigReal
Definition: common.h:112
void build_go_params(StringList *)
Definition: GoMolecule.C:79
void build_go_sigmas2(StringList *, char *)
Definition: GoMolecule.C:747
int numExclusions
Definition: Molecule.h:570