Fix REST2 config parameter name lookup
[namd.git] / src / NamdState.C
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
34 NamdState::NamdState()
35 {
36     configList = NULL;
37     simParameters = NULL;
38     parameters = NULL;
39     molecule = NULL;
40     pdb = NULL;
41 }
42
43 int
44 NamdState::status()
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
74 void NamdState:: useController(Controller *controllerPtr)
75 {
76   controller=controllerPtr;
77 }
78
79 void NamdState:: runController(void)
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
123 int NamdState::configListInit(ConfigList *cfgList) {
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         }
503 //fepe
504         if (simParameters->lesOn) {
505            if (simParameters->alchOn) NAMD_bug("FEP/TI and LES are incompatible!");
506            molecule->build_fep_flags(configList->find("lesfile"),
507                 configList->find("lescol"), pdb, NULL, "les");
508         }
509         if (simParameters->soluteScalingOn) {
510            molecule->build_ss_flags(configList->find("soluteScalingFile"),
511                 configList->find("soluteScalingCol"), pdb, NULL);
512         }
513         if (simParameters->pairInteractionOn) {
514            molecule->build_fep_flags(configList->find("pairInteractionFile"),
515                 configList->find("pairInteractionCol"), pdb, NULL, "pairInteraction");
516         }      
517         if (simParameters->pressureProfileAtomTypes > 1) {
518           molecule->build_fep_flags(configList->find("pressureProfileAtomTypesFile"),
519                 configList->find("pressureProfileAtomTypesCol"), pdb, NULL, "pressureProfileAtomTypes");
520         }
521        
522         #ifdef OPENATOM_VERSION
523         if (simParameters->openatomOn) {
524           molecules->build_qmmm_flags(configList->find("openatomPdbFile",
525                 configList->find("openatomPdbCol"), pdb, NULL, "openatomPdb")
526         }
527         #endif // OPENATOM_VERSION
528
529         if (simParameters->qmForcesOn){
530             
531 #ifdef MEM_OPT_VERSION
532             NAMD_die("QM forces are not supported in memory-optimized builds.");
533 #endif
534             
535 #ifdef NAMD_CUDA
536            NAMD_die("QM forces are not compatible with CUDA at this time");
537 #endif
538             
539             molecule->set_qm_replaceAll(simParameters->qmReplaceAll);
540             
541             if (simParameters->qmParamPDBDefined)
542                 molecule->prepare_qm(simParameters->qmParamPDB,
543                                           parameters, configList);
544             else if (pdbFilename)
545                 molecule->prepare_qm(pdbFilename,
546                                           parameters, configList);
547             else
548                 molecule->prepare_qm(configList->find("coordinates")->data,
549                                           parameters, configList);
550             
551         }
552         
553         
554         if (simParameters->LJcorrection) {
555           molecule->compute_LJcorrection();
556         }
557 #endif
558
559         // JLai checks to see if Go Forces are turned on
560         if (simParameters->goForcesOn) {
561 #ifdef MEM_OPT_VERSION
562           NAMD_die("Go forces are not supported in memory-optimized builds.");
563 #else
564           StringList *moleculeFilename = configList->find("structure");
565           StringList *parameterFilename = configList->find("parameters");
566           StringList *goFilename = configList->find("goParameters");
567           StringList *goStructureFilename = configList->find("goCoordinates");
568           
569           // Added by JLai -- 1.10.12 -- Code to build the Go parameters (from within the Go molecule instead of the parameters object)
570           molecule->build_go_params(goFilename);
571           // Added by JLai -- 6.3.11 -- flag to switch between the different goMethodologies
572           int goMethod = simParameters->goMethod;
573           if (goMethod == 1) { // should probably replace with switch statement
574             iout << iINFO << "Using Go method matrix\n" << endi;
575             molecule->build_go_sigmas(goStructureFilename, NULL);
576           } else if (goMethod == 2) {
577             iout << iINFO << "Using Go method jump table\n" << endi;
578             molecule->build_go_sigmas2(goStructureFilename, NULL);
579           } else if (goMethod == 3) {
580             iout << iINFO << "Using Go method lowmem\n" << endi;
581             molecule->build_go_arrays(goStructureFilename, NULL);
582           } else {
583             NAMD_die("Failed to read goMethod variable in NamdState.C");
584           }
585 #endif
586         }
587         // End of Go code -- JLai
588
589 #ifndef MEM_OPT_VERSION
590         iout << iINFO << "****************************\n";
591         iout << iINFO << "STRUCTURE SUMMARY:\n";
592         iout << iINFO << molecule->numAtoms << " ATOMS\n";
593         iout << iINFO << molecule->numBonds << " BONDS\n";
594         iout << iINFO << molecule->numAngles << " ANGLES\n";
595         iout << iINFO << molecule->numDihedrals << " DIHEDRALS\n";
596         iout << iINFO << molecule->numImpropers << " IMPROPERS\n";
597         iout << iINFO << molecule->numCrossterms << " CROSSTERMS\n";
598         iout << iINFO << molecule->numExclusions << " EXCLUSIONS\n";
599
600         //****** BEGIN CHARMM/XPLOR type changes
601         if ((molecule->numMultipleDihedrals) && (simParameters->paraTypeXplorOn))
602         {
603                 iout << iINFO << molecule->numMultipleDihedrals 
604              << " DIHEDRALS WITH MULTIPLE PERIODICITY (BASED ON PSF FILE)\n";
605         }
606         if ((molecule->numMultipleDihedrals) && (simParameters->paraTypeCharmmOn))
607         {
608                 iout << iINFO << molecule->numMultipleDihedrals 
609          << " DIHEDRALS WITH MULTIPLE PERIODICITY IGNORED (BASED ON PSF FILE) \n";
610                 iout << iINFO  
611          << " CHARMM MULTIPLICITIES BASED ON PARAMETER FILE INFO! \n";
612         }
613         //****** END CHARMM/XPLOR type changes
614
615         if (molecule->numMultipleImpropers)
616         {
617                 iout << iINFO << molecule->numMultipleImpropers 
618                          << " IMPROPERS WITH MULTIPLE PERIODICITY\n";
619         }
620         
621         if (simParameters->constraintsOn)
622         {
623            iout << iINFO << molecule->numConstraints << " CONSTRAINTS\n";
624         }
625
626         if (simParameters->consForceOn)
627           iout << iINFO << molecule->numConsForce << " CONSTANT FORCES\n";
628
629         if (simParameters->stirOn)
630           iout << iINFO << molecule->numStirredAtoms << " STIRRED ATOMS\n";
631
632         if (simParameters->fixedAtomsOn)
633         {
634            iout << iINFO << molecule->numFixedAtoms << " FIXED ATOMS\n";
635         }
636
637         if (simParameters->rigidBonds)
638         {
639            iout << iINFO << molecule->numRigidBonds << " RIGID BONDS\n";
640         }
641
642         if (simParameters->fixedAtomsOn && simParameters->rigidBonds)
643         {
644            iout << iINFO << molecule->numFixedRigidBonds <<
645                         " RIGID BONDS BETWEEN FIXED ATOMS\n";
646         }
647         
648         /* BEGIN gf */
649         if (simParameters->mgridforceOn)
650         {
651             int i;
652             iout << iINFO << molecule->numGridforceGrids 
653                  << " GRIDS ACTIVE\n";
654         }
655         /* END gf */
656
657 //Modifications for alchemical fep
658         if (simParameters->alchOn) {
659           iout << iINFO << "ALCH: " 
660                << molecule->numFepInitial <<
661                " ATOMS TO DISAPPEAR IN FINAL STATE\n";
662            iout << iINFO << "ALCH: " 
663                <<  molecule->numFepFinal <<
664                " ATOMS TO APPEAR IN FINAL STATE\n";
665            if (molecule->suspiciousAlchBonds) {
666              iout << iWARN << "ALCH: SUSPICIOUS BONDS BETWEEN INITIAL AND " <<
667              "FINAL GROUPS WERE FOUND" << "\n" << endi;
668            }
669            if (molecule->alchDroppedAngles) {
670              iout << iINFO << "ALCH: " 
671                  << molecule->alchDroppedAngles <<
672                  " ANGLES LINKING INITIAL AND FINAL ATOMS DELETED\n";
673            }
674            if (molecule->alchDroppedDihedrals) {
675              iout << iINFO << "ALCH: "
676                  << molecule->alchDroppedDihedrals <<
677                  " DIHEDRALS LINKING INITIAL AND FINAL ATOMS DELETED\n";
678            }
679            if (molecule->alchDroppedImpropers) {
680              iout << iINFO << "ALCH: "
681                  << molecule->alchDroppedImpropers <<
682                  " IMPROPERS LINKING INITIAL AND FINAL ATOMS DELETED\n";
683            }
684         }
685 //fepe
686
687         if (simParameters->lesOn) {
688            iout << iINFO << molecule->numFepInitial <<
689                " LOCALLY ENHANCED ATOMS ENABLED\n";
690         }
691
692         if (simParameters->soluteScalingOn) {
693            iout << iINFO << " SOLUTE SCALING ENABLED\n";
694         }
695        
696         if (simParameters->pairInteractionOn) {
697            iout << iINFO << "PAIR INTERACTION GROUP 1 CONTAINS "
698                 <<  molecule->numFepInitial << " ATOMS\n";
699            if (!simParameters->pairInteractionSelf) {
700              iout << iINFO << "PAIR INTERACTION GROUP 2 CONTAINS "
701                   <<  molecule->numFepFinal << " ATOMS\n";
702            }
703         }
704            
705 #if 1
706         if (molecule->numLonepairs != 0) {
707           iout << iINFO << molecule->numLonepairs << " LONE PAIRS\n";
708         }
709         if (molecule->numDrudeAtoms != 0) {
710           iout << iINFO << molecule->numDrudeAtoms << " DRUDE ATOMS\n";
711         }
712         iout << iINFO << molecule->num_deg_freedom(1)
713              << " DEGREES OF FREEDOM\n";
714         if (simParameters->drudeOn) {
715           int g_bond = 3 * molecule->numDrudeAtoms;
716           int g_com = molecule->num_deg_freedom(1) - g_bond;
717           iout << iINFO << g_com << " DRUDE COM DEGREES OF FREEDOM\n";
718           iout << iINFO << g_bond << " DRUDE BOND DEGREES OF FREEDOM\n";
719         }
720 #endif
721 #if 0
722         {
723           // Copied from Controller::printEnergies()
724           int64 numAtoms = molecule->numAtoms;
725           int64 numDegFreedom = 3 * numAtoms;
726     int numLonepairs = molecule->numLonepairs;
727           int numFixedAtoms = molecule->numFixedAtoms;
728           if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
729           if ( ! ( numFixedAtoms || molecule->numConstraints
730                 || simParameters->comMove || simParameters->langevinOn ) ) {
731             numDegFreedom -= 3;
732           }
733     if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
734           int numRigidBonds = molecule->numRigidBonds;
735           int numFixedRigidBonds = molecule->numFixedRigidBonds;
736     // numLonepairs is subtracted here because all lonepairs have a rigid bond
737     // to oxygen, but all of the LP degrees of freedom are dealt with above
738     numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
739           iout << iINFO << numDegFreedom << " DEGREES OF FREEDOM\n";
740         }
741 #endif
742
743         iout << iINFO << molecule->numHydrogenGroups << " HYDROGEN GROUPS\n";
744         iout << iINFO << molecule->maxHydrogenGroupSize
745                 << " ATOMS IN LARGEST HYDROGEN GROUP\n";
746         iout << iINFO << molecule->numMigrationGroups << " MIGRATION GROUPS\n";
747         iout << iINFO << molecule->maxMigrationGroupSize
748                 << " ATOMS IN LARGEST MIGRATION GROUP\n";
749         if (simParameters->fixedAtomsOn)
750         {
751            iout << iINFO << molecule->numFixedGroups <<
752                         " HYDROGEN GROUPS WITH ALL ATOMS FIXED\n";
753         }
754
755         {
756           BigReal totalMass = 0;
757           BigReal totalCharge = 0;
758           int i;
759           for ( i = 0; i < molecule->numAtoms; ++i ) {
760             totalMass += molecule->atommass(i);
761             totalCharge += molecule->atomcharge(i);
762           }
763           iout << iINFO << "TOTAL MASS = " << totalMass << " amu\n"; 
764           iout << iINFO << "TOTAL CHARGE = " << totalCharge << " e\n"; 
765
766           BigReal volume = lattice.volume();
767           if ( volume ) {
768             iout << iINFO << "MASS DENSITY = "
769               << ((totalMass/volume) / 0.6022) << " g/cm^3\n";
770             iout << iINFO << "ATOM DENSITY = "
771               << (molecule->numAtoms/volume) << " atoms/A^3\n";
772           }
773         }
774
775         iout << iINFO << "*****************************\n";
776         iout << endi;
777         fflush(stdout);
778
779   StringList *binCoordinateFilename = configList->find("bincoordinates");
780   if ( binCoordinateFilename && ! reload ) {
781     read_binary_coors(binCoordinateFilename->data, pdb);
782   }
783
784   DebugM(4, "::configFileInit() - printing Molecule Information\n");
785
786   molecule->print_atoms(parameters);
787   molecule->print_bonds(parameters);
788   molecule->print_exclusions();
789   fflush(stdout);
790 #endif
791
792   DebugM(4, "::configFileInit() - done printing Molecule Information\n");
793   DebugM(1, "::configFileInit() - done\n");
794
795   return(0);
796 }
797