NamdState.C

Go to the documentation of this file.
00001 /*
00002 ***  Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
00003 ***  The Board of Trustees of the University of Illinois.
00004 ***  All rights reserved.
00005 **/
00006 
00007 /*
00008    Holds pointers to large molecule data structure, simulation parameters...
00009 */
00010 
00011 #include "InfoStream.h"
00012 #include "common.h"
00013 #include "Molecule.h"
00014 #include "Parameters.h"
00015 #include "SimParameters.h"
00016 #include "ConfigList.h"
00017 #include "PDB.h"
00018 #include "NamdState.h"
00019 #include "Controller.h"
00020 #include "ScriptTcl.h"
00021 #ifndef WIN32
00022 #include <unistd.h>
00023 #endif
00024 #include <sys/stat.h>
00025 #include "parm.h"
00026 
00027 //#define DEBUGM
00028 #include "Debug.h"
00029 
00030 #include "CompressPsf.h"
00031 #include "PluginIOMgr.h"
00032 #include "BackEnd.h"
00033 
00034 NamdState::NamdState()
00035 {
00036     configList = NULL;
00037     simParameters = NULL;
00038     parameters = NULL;
00039     molecule = NULL;
00040     pdb = NULL;
00041 }
00042 
00043 int
00044 NamdState::status()
00045 {
00046     int ret=0;
00047     if (configList != NULL) {
00048       DebugM(1, "Config List exists\n");
00049     } else ret++;
00050 
00051     if (simParameters != NULL) {
00052       DebugM(1, "SimParameters exists\n");
00053     }
00054     else ret++;
00055 
00056     if (parameters != NULL) {
00057       DebugM(1,"Parameters exists\n");
00058     }
00059     else ret++;
00060 
00061     if (molecule != NULL) {
00062       DebugM(1, "Molecule exists\n");
00063     }
00064     else ret++;
00065 
00066     if (pdb != NULL) {
00067       DebugM(1,"PDB exists \n");
00068     }
00069     else ret++;
00070     
00071     return(ret);
00072 }
00073 
00074 void NamdState:: useController(Controller *controllerPtr)
00075 {
00076   controller=controllerPtr;
00077 }
00078 
00079 void NamdState:: runController(void)
00080 {
00081   controller->run();
00082 }
00083 
00084 extern void read_binary_coors(char *fname, PDB *pdbobj);
00085 
00086 #ifdef MEM_OPT_VERSION
00087  //Check features that are not supported in the memory optimized 
00088  //version.  --Chao Mei
00089 void NamdState::checkMemOptCompatibility(){
00090     if(simParameters->genCompressedPsf) {
00091         iout << "In the memory optimized version, the compression of molecule information is not supported! "
00092                 << "Please use the non-memory optimized version.\n" <<endi;
00093         NAMD_die("MEMOPT: unsupported usage");
00094     }
00095 
00096     if(simParameters->langevinOn && simParameters->langevinDamping == 0.0) {
00097         iout << iWARN << "langevinDamping MUST NOT BE 0.0 IF LANGEVIN IS"
00098             <<" TURNED ON IN MEMORY OPTIMIZED VERSION\n" <<endi;
00099         NAMD_die("MEMOPT: unsupported feature");
00100     }
00101     if(simParameters->tCoupleOn)
00102         NAMD_die("MEMOPT: tCouple is not supported in memory optimized version");
00103     if(simParameters->pairInteractionOn)
00104         NAMD_die("MEMOPT: pairInteractionOn could not be enabled in memory optimized version");
00105     if(simParameters->alchFepOn || simParameters->alchOn){
00106         iout << iWARN << "ALCH: AUTOMATIC DELETION OF BONDED INTERACTIONS "
00107         << "BETWEEN INITIAL AND FINAL GROUPS IS NOT SUPPORTED IN MEMORY "
00108         << "OPTIMISED VERSION - MANUAL PROCESSING IS NECESSARY\n" << endi;
00109         NAMD_die("MEMOPT: unsupported feature");
00110     }
00111     if(simParameters->alchThermIntOn)
00112         NAMD_die("MEMOPT: alchThermIntOn could not be enabled in memory optimized version");
00113     if(simParameters->lesOn)
00114         NAMD_die("MEMOPT: lesOn could not be enabled in memory optimized version");
00115     if(simParameters->lonepairs) {
00116         NAMD_die("MEMOPT: lonepairs could not be enabled in memory optimized version");
00117     }
00118 }
00119 #endif
00120 
00121 int NamdState::configListInit(ConfigList *cfgList) {
00122   configList = cfgList;
00123   if (!configList->okay()) {
00124     NAMD_die("Simulation config file is incomplete or contains errors.");
00125   }
00126   DebugM(1,"NamdState::configFileInit configList okay\n");
00127 
00128   char *currentdir = 0;
00129   simParameters =  new SimParameters(configList,currentdir);
00130   fflush(stdout);
00131   lattice = simParameters->lattice;
00132 
00133  //Check features that are not supported in the memory optimized 
00134  //version.  --Chao Mei
00135 #ifdef MEM_OPT_VERSION
00136   checkMemOptCompatibility();
00137 #endif
00138 
00139   //Check rigidBonds type when generating the compressed psf files.
00140   if(simParameters->genCompressedPsf) {
00141       if(simParameters->rigidBonds == RIGID_NONE){
00142           //default to RIGID_ALL
00143           simParameters->rigidBonds = RIGID_ALL;
00144       }
00145   }
00146 
00147   return loadStructure(0,0,0);
00148 }
00149 
00150 int NamdState::loadStructure(const char *molFilename, const char *pdbFilename, int reload) {
00151 
00152   StringList *molInfoFilename;
00153   // If it's AMBER force field, read the AMBER style files;
00154   // if it's GROMACS, read the GROMACS files;
00155   // Otherwise read the CHARMM style files
00156 
00157   if (simParameters->amberOn) {
00158     if ( reload ) NAMD_die("Molecular structure reloading not supported for Amber input files.\n");
00159     StringList *parmFilename = configList->find("parmfile");
00160     molInfoFilename = parmFilename;
00161     StringList *coorFilename = configList->find("ambercoor");
00162     // "amber" is a temporary data structure, which records all
00163     // the data from the parm file. After copying them into
00164     // molecule, parameter and pdb structures, it will be deleted.
00165     Ambertoppar *amber;
00166     amber = new Ambertoppar;
00167     if (amber->readparm(parmFilename->data))
00168     { parameters = new Parameters(amber, simParameters->vdwscale14);
00169       molecule = new Molecule(simParameters, parameters, amber);
00170       if (coorFilename != NULL)
00171         pdb = new PDB(coorFilename->data,amber);
00172       delete amber;
00173     }
00174     else
00175       NAMD_die("Failed to read AMBER parm file!");
00176     parameters->print_param_summary();
00177   }
00178   else if (simParameters->gromacsOn) {
00179     if ( reload ) NAMD_die("Molecular structure reloading not supported for Gromacs input files.\n");
00180     StringList *topFilename = configList->find("grotopfile");
00181     molInfoFilename = topFilename;
00182     StringList *coorFilename = configList->find("grocoorfile");
00183     // "gromacsFile" is a temporary data structure, which records all
00184     // the data from the topology file. After copying it into the
00185     // molecule and parameter and pdb, it will be deleted.
00186     GromacsTopFile *gromacsFile;
00187     gromacsFile = new GromacsTopFile(topFilename->data);
00188     parameters = new Parameters(gromacsFile,simParameters->minimizeCGOn);
00189     if (coorFilename != NULL)
00190       pdb = new PDB(coorFilename->data,gromacsFile);
00191 
00192     molecule = new Molecule(simParameters, parameters, gromacsFile);
00193     // XXX does Molecule(needAll,these,arguments)?
00194 
00195     delete gromacsFile; // XXX unimplemented
00196 
00197     // XXX add error handling when the file doesn't exist
00198     // XXX make sure the right things happen when the parameters are
00199     // not even specified.
00200     // NAMD_die("Failed to read AMBER parm file!");
00201     parameters->print_param_summary();
00202   }
00203   else if (simParameters->usePluginIO){
00204 #ifdef MEM_OPT_VERSION          
00205         NAMD_die("Using plugin IO is not supported in memory optimized version!");
00206 #else    
00207     if ( pdbFilename ) {
00208       NAMD_bug("NamdState::loadStructure pdbFilename non-null with usePluginIO\n");
00209     }
00210 
00211     PluginIOMgr *pIOMgr = new PluginIOMgr();
00212     
00213     iout << iWARN << "Plugin-based I/O is still in development and may still have bugs\n" << endi;
00214 
00215     molfile_plugin_t *pIOHandle = pIOMgr->getPlugin();
00216     if (pIOHandle == NULL) {
00217         NAMD_die("ERROR: Failed to match requested plugin type");
00218     }
00219     if ( pIOHandle->open_file_read == NULL )
00220        NAMD_die("ERROR: Selected plugin type cannot open files"); 
00221     if ( pIOHandle->read_structure == NULL )
00222        NAMD_die("ERROR: Selected plugin type cannot read structures"); 
00223     if ( pIOHandle->read_next_timestep == NULL )
00224        NAMD_die("ERROR: Selected plugin type cannot read coordinates"); 
00225 
00226     StringList *moleculeFilename = configList->find("structure");
00227     molInfoFilename = moleculeFilename;
00228     if ( ! molFilename ) molFilename = moleculeFilename->data;
00229   if ( ! reload ) {
00230     StringList *parameterFilename = configList->find("parameters");
00231     //****** BEGIN CHARMM/XPLOR type changes
00232     // For AMBER use different constructor based on parm_struct!!!  -JCP
00233     parameters = new Parameters(simParameters, parameterFilename);
00234     parameters->print_param_summary();
00235   }
00236 
00237     int numAtoms = 0;
00238     //TODO: not sure about the name field in the handler
00239     void *plgFile = pIOHandle->open_file_read(molFilename, 
00240                                               pIOHandle->name, &numAtoms);
00241     if(plgFile ==  NULL) {
00242         NAMD_die("ERROR: Opening structure file failed!");
00243     }
00244 
00245     double fileReadTime = CmiWallTimer();
00246     molecule = new Molecule(simParameters, parameters, pIOHandle, plgFile, numAtoms);
00247     iout << iINFO << "TIME FOR LOAD MOLECULE STRUCTURE INFORMATION: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00248 
00249     /* If we are only generating compressed molecule information, the PDB object is not needed */
00250     if(!simParameters->genCompressedPsf) {        
00251         fileReadTime = CmiWallTimer();
00252         //get the occupancy data from the Molecule object and then free it
00253         //as it is stored in the Molecule object.
00254         pdb = new PDB(pIOHandle, plgFile, molecule->numAtoms, molecule->getOccupancyData(), molecule->getBFactorData());
00255         molecule->freeOccupancyData();
00256         molecule->freeBFactorData();
00257         iout << iINFO << "TIME FOR LOADING ATOMS' COORDINATES INFORMATION: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00258     }
00259 
00260     pIOHandle->close_file_read(plgFile);
00261     delete pIOMgr;
00262 #endif    
00263   }
00264   else
00265   { 
00266     StringList *moleculeFilename = configList->find("structure");
00267     molInfoFilename = moleculeFilename; 
00268     if ( ! molFilename ) molFilename = moleculeFilename->data;
00269   if ( ! reload ) {
00270     StringList *parameterFilename = configList->find("parameters");
00271     //****** BEGIN CHARMM/XPLOR type changes
00272     // For AMBER use different constructor based on parm_struct!!!  -JCP
00273     parameters = new Parameters(simParameters, parameterFilename);
00274     //****** END CHARMM/XPLOR type changes    
00275 
00276     parameters->print_param_summary();
00277   }
00278 
00279     double fileReadTime = CmiWallTimer();
00280     molecule = new Molecule(simParameters, parameters, (char*)molFilename, configList);
00281     iout << iINFO << "TIME FOR READING PSF FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00282 }
00283 
00284   fflush(stdout);
00285 
00286 #ifdef MEM_OPT_VERSION
00287   //upon knowing the number of atoms, it's good time to estimate the number of
00288   //input/output processors if their value is not set
00289   if(simParameters->numinputprocs==0){
00290     int numatoms = molecule->numAtoms;
00291     long estval = (sizeof(InputAtom)+2*sizeof(int)+1)*((long)(numatoms));
00292     int numprocs = estval>>26; //considering every input proc consumes about 64M.
00293     if(numprocs==0){
00294         numprocs=1;
00295     }else if(numprocs>CkNumPes()){
00296         numprocs=CkNumPes();
00297     }
00298     simParameters->numinputprocs=numprocs;
00299   }
00300   if(simParameters->numoutputprocs==0){
00301     int numatoms = molecule->numAtoms;
00302     long estval = (sizeof(Vector)*2)*((long)(numatoms));
00303     int numprocs = estval>>26; //considering every input proc consumes about 64M.
00304     if(numprocs==0){
00305       numprocs=1;
00306     }else if(numprocs>CkNumPes()){
00307       numprocs=CkNumPes();
00308     }
00309     simParameters->numoutputprocs=numprocs;    
00310   }
00311   //check the number of output procs that simultaneously write to a file
00312   if(simParameters->numoutputwrts > simParameters->numoutputprocs) {
00313       simParameters->numoutputwrts = simParameters->numoutputprocs;
00314   }
00315 
00316   if (simParameters->fixedAtomsOn){
00317       double fileReadTime = CmiWallTimer();
00318       molecule->load_fixed_atoms(configList->find("fixedAtomListFile"));
00319       iout << iINFO << "TIME FOR READING FIXED ATOMS FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00320   }
00321 
00322   if (simParameters->constraintsOn){
00323       double fileReadTime = CmiWallTimer();
00324       molecule->load_constrained_atoms(configList->find("consAtomListFile"));
00325       iout << iINFO << "TIME FOR READING CONSTRAINED ATOMS FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00326   }
00327 #else
00328   if (simParameters->extraBondsOn) {        
00329     //The extra bonds building will be executed in read_compressed_psf in
00330     //the memory optimized version, so avoid calling this function in the 
00331     //memory optimized run.
00332     if(!simParameters->useCompressedPsf)
00333       molecule->build_extra_bonds(parameters, configList->find("extraBondsFile"));         
00334   }
00335   if(simParameters->genCompressedPsf) {
00336       double fileReadTime = CmiWallTimer();
00337       compress_molecule_info(molecule, molInfoFilename->data, parameters, simParameters, configList);
00338       iout << "Finished compressing molecule information, which takes " << CmiWallTimer()-fileReadTime <<"(s)\n"<<endi;
00339       BackEnd::exit();
00340   }
00341 
00342   //If using plugin-based IO, the PDB object is already created!
00343   StringList *coordinateFilename = NULL;
00344   if(!simParameters->usePluginIO) {
00345       //In the memory opt version, the coordinates of atoms
00346       //are read during startup in parallel with a bincoordinates input
00347       //-Chao Mei
00348       double fileReadTime = CmiWallTimer();
00349       if ( pdbFilename ) {
00350         iout << iINFO << "Reading pdb file " << pdbFilename << "\n" << endi;
00351         pdb = new PDB(pdbFilename);
00352       } else {
00353         coordinateFilename = configList->find("coordinates");    
00354         if (coordinateFilename != NULL) {
00355           iout << iINFO << "Reading pdb file " << coordinateFilename->data << "\n" << endi;
00356           pdb = new PDB(coordinateFilename->data);
00357         }
00358       }
00359       if (pdb->num_atoms() != molecule->numAtoms) {
00360         NAMD_die("Number of pdb and psf atoms are not the same!");
00361       }
00362       iout << iINFO << "TIME FOR READING PDB FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00363       iout << iINFO << "\n" << endi;
00364   }
00365 
00366         //  If constraints are active, build the parameters necessary
00367         if (simParameters->constraintsOn)
00368         {
00369            StringList *consRefFile = configList->find("consref");
00370            StringList *consKFile = configList->find("conskfile");
00371 
00372           if (coordinateFilename != NULL) {
00373            if(strcasecmp(coordinateFilename->data, consRefFile->data)==0)
00374                 consRefFile = NULL;
00375            if(strcasecmp(coordinateFilename->data, consKFile->data)==0)
00376                 consKFile = NULL;
00377           }
00378 
00379            molecule->build_constraint_params(consRefFile, consKFile,
00380                                              configList->find("conskcol"),
00381                                              pdb,
00382                                              NULL);
00383         }
00384 #endif
00385         //CkPrintf ("DEBUG--check if StirOn to build stir params..\n");
00386 
00387         if (simParameters->stirOn)
00388         {       
00389         //CkPrintf ("DEBUG--now to build stir params..\n");
00390           
00391            molecule->build_stirred_atoms(configList->find("stirFilename"),
00392                                        configList->find("stirredAtomsCol"),
00393                                        pdb,
00394                                        NULL);
00395         }
00396 
00397 
00398 #ifndef MEM_OPT_VERSION
00399         if (simParameters->fixedAtomsOn)
00400         {
00401            molecule->build_fixed_atoms(configList->find("fixedatomsfile"),
00402                                         configList->find("fixedatomscol"),
00403                                         pdb,
00404                                         NULL);
00405         }
00406 #endif
00407         
00408         /* BEGIN gf */
00409         if (simParameters->mgridforceOn)
00410         {
00411             molecule->build_gridforce_params(configList->find("gridforcefile"),
00412                                              configList->find("gridforcecol"),
00413                                              configList->find("gridforcechargecol"),
00414                                              configList->find("gridforcepotfile"),
00415                                              pdb,
00416                                              NULL);
00417         }
00418         /* END gf */
00419 
00420         // If constant forces are active, read the forces necessary
00421         if (simParameters->consForceOn) {
00422     char *filename = NULL;
00423     if (configList->find("consforcefile"))
00424       filename = configList->find("consforcefile")->data;
00425     molecule->build_constant_forces(filename);
00426   }
00427 
00428         if (simParameters->excludeFromPressure) {
00429            molecule->build_exPressure_atoms(
00430              configList->find("excludeFromPressureFile"),
00431              configList->find("excludeFromPressureCol"),
00432              pdb, NULL);
00433         }
00434 
00435         // If moving drag is active, build the parameters necessary
00436         if (simParameters->movDragOn) {
00437           molecule->build_movdrag_params(configList->find("movDragFile"),
00438                                          configList->find("movDragCol"),
00439                                          configList->find("movDragVelFile"),
00440                                          pdb,
00441                                          NULL);
00442         }
00443 
00444         // If rotating drag is active, build the parameters necessary
00445         if (simParameters->rotDragOn) {
00446           molecule->build_rotdrag_params(configList->find("rotDragFile"),
00447                                          configList->find("rotDragCol"),
00448                                          configList->find("rotDragAxisFile"),
00449                                          configList->find("rotDragPivotFile"),
00450                                          configList->find("rotDragVelFile"),
00451                                          configList->find("rotDragVelCol"),
00452                                          pdb,
00453                                          NULL);
00454         }
00455 
00456         // If "constant" torque is active, build the parameters necessary
00457         if (simParameters->consTorqueOn) {
00458           molecule->build_constorque_params(configList->find("consTorqueFile"),
00459                                        configList->find("consTorqueCol"),
00460                                        configList->find("consTorqueAxisFile"),
00461                                        configList->find("consTorquePivotFile"),
00462                                        configList->find("consTorqueValFile"),
00463                                        configList->find("consTorqueValCol"),
00464                                        pdb,
00465                                        NULL);
00466         }
00467 
00468 #ifndef MEM_OPT_VERSION
00469         //  If langevin dynamics or temperature coupling are active, build 
00470         //  the parameters necessary
00471         if (simParameters->langevinOn)
00472         {
00473           if (simParameters->langevinDamping == 0.0) {
00474             molecule->build_langevin_params(configList->find("langevinfile"),
00475                                             configList->find("langevincol"),
00476                                             pdb,
00477                                             NULL);
00478           } else {
00479             molecule->build_langevin_params(simParameters->langevinDamping,
00480                                             simParameters->drudeDamping,
00481                                             simParameters->langevinHydrogen);
00482           }
00483         }
00484         else if (simParameters->tCoupleOn)
00485         {
00486            //  Temperature coupling uses the same parameters, but with different
00487            //  names . . .
00488            molecule->build_langevin_params(configList->find("tcouplefile"),
00489                                             configList->find("tcouplecol"),
00490                                             pdb,
00491                                             NULL);
00492         }
00493 
00494 //Modifications for alchemical fep
00495      //identify the mutant atoms for fep simulation
00496         if (simParameters->alchOn) {
00497            molecule->build_fep_flags(configList->find("alchfile"),
00498                 configList->find("alchcol"), pdb, NULL, "alch" );
00499            molecule->delete_alch_bonded();
00500         }
00501 //fepe
00502         if (simParameters->lesOn) {
00503            if (simParameters->alchOn) NAMD_bug("FEP/TI and LES are incompatible!");
00504            molecule->build_fep_flags(configList->find("lesfile"),
00505                 configList->find("lescol"), pdb, NULL, "les");
00506         }
00507         if (simParameters->pairInteractionOn) {
00508            molecule->build_fep_flags(configList->find("pairInteractionFile"),
00509                 configList->find("pairInteractionCol"), pdb, NULL, "pairInteraction");
00510         }      
00511         if (simParameters->pressureProfileAtomTypes > 1) {
00512           molecule->build_fep_flags(configList->find("pressureProfileAtomTypesFile"),
00513                 configList->find("pressureProfileAtomTypesCol"), pdb, NULL, "pressureProfileAtomTypes");
00514         }
00515        
00516         #ifdef OPENATOM_VERSION
00517         if (simParameters->openatomOn) {
00518           molecules->build_qmmm_flags(configList->find("openatomPdbFile",
00519                 configList->find("openatomPdbCol"), pdb, NULL, "openatomPdb")
00520         }
00521         #endif // OPENATOM_VERSION
00522 
00523         if (simParameters->qmForcesOn){
00524             
00525 #ifdef MEM_OPT_VERSION
00526             NAMD_die("QM forces are not supported in memory-optimized builds.");
00527 #endif
00528             
00529 #ifdef NAMD_CUDA
00530            NAMD_die("QM forces are not compatible with CUDA at this time");
00531 #endif
00532             
00533             molecule->set_qm_replaceAll(simParameters->qmReplaceAll);
00534             
00535             if (simParameters->qmParamPDBDefined)
00536                 molecule->prepare_qm(simParameters->qmParamPDB,
00537                                           parameters, configList);
00538             else if (pdbFilename)
00539                 molecule->prepare_qm(pdbFilename,
00540                                           parameters, configList);
00541             else
00542                 molecule->prepare_qm(configList->find("coordinates")->data,
00543                                           parameters, configList);
00544             
00545         }
00546         
00547         
00548         if (simParameters->LJcorrection) {
00549           molecule->compute_LJcorrection();
00550         }
00551 #endif
00552 
00553         // JLai checks to see if Go Forces are turned on
00554         if (simParameters->goForcesOn) {
00555 #ifdef MEM_OPT_VERSION
00556           NAMD_die("Go forces are not supported in memory-optimized builds.");
00557 #else
00558           StringList *moleculeFilename = configList->find("structure");
00559           StringList *parameterFilename = configList->find("parameters");
00560           StringList *goFilename = configList->find("goParameters");
00561           StringList *goStructureFilename = configList->find("goCoordinates");
00562           
00563           // Added by JLai -- 1.10.12 -- Code to build the Go parameters (from within the Go molecule instead of the parameters object)
00564           molecule->build_go_params(goFilename);
00565           // Added by JLai -- 6.3.11 -- flag to switch between the different goMethodologies
00566           int goMethod = simParameters->goMethod;
00567           if (goMethod == 1) { // should probably replace with switch statement
00568             iout << iINFO << "Using Go method matrix\n" << endi;
00569             molecule->build_go_sigmas(goStructureFilename, NULL);
00570           } else if (goMethod == 2) {
00571             iout << iINFO << "Using Go method jump table\n" << endi;
00572             molecule->build_go_sigmas2(goStructureFilename, NULL);
00573           } else if (goMethod == 3) {
00574             iout << iINFO << "Using Go method lowmem\n" << endi;
00575             molecule->build_go_arrays(goStructureFilename, NULL);
00576           } else {
00577             NAMD_die("Failed to read goMethod variable in NamdState.C");
00578           }
00579 #endif
00580         }
00581         // End of Go code -- JLai
00582 
00583 #ifndef MEM_OPT_VERSION
00584         iout << iINFO << "****************************\n";
00585         iout << iINFO << "STRUCTURE SUMMARY:\n";
00586         iout << iINFO << molecule->numAtoms << " ATOMS\n";
00587         iout << iINFO << molecule->numBonds << " BONDS\n";
00588         iout << iINFO << molecule->numAngles << " ANGLES\n";
00589         iout << iINFO << molecule->numDihedrals << " DIHEDRALS\n";
00590         iout << iINFO << molecule->numImpropers << " IMPROPERS\n";
00591         iout << iINFO << molecule->numCrossterms << " CROSSTERMS\n";
00592         iout << iINFO << molecule->numExclusions << " EXCLUSIONS\n";
00593 
00594         //****** BEGIN CHARMM/XPLOR type changes
00595         if ((molecule->numMultipleDihedrals) && (simParameters->paraTypeXplorOn))
00596         {
00597                 iout << iINFO << molecule->numMultipleDihedrals 
00598              << " DIHEDRALS WITH MULTIPLE PERIODICITY (BASED ON PSF FILE)\n";
00599         }
00600         if ((molecule->numMultipleDihedrals) && (simParameters->paraTypeCharmmOn))
00601         {
00602                 iout << iINFO << molecule->numMultipleDihedrals 
00603          << " DIHEDRALS WITH MULTIPLE PERIODICITY IGNORED (BASED ON PSF FILE) \n";
00604                 iout << iINFO  
00605          << " CHARMM MULTIPLICITIES BASED ON PARAMETER FILE INFO! \n";
00606         }
00607         //****** END CHARMM/XPLOR type changes
00608 
00609         if (molecule->numMultipleImpropers)
00610         {
00611                 iout << iINFO << molecule->numMultipleImpropers 
00612                          << " IMPROPERS WITH MULTIPLE PERIODICITY\n";
00613         }
00614         
00615         if (simParameters->constraintsOn)
00616         {
00617            iout << iINFO << molecule->numConstraints << " CONSTRAINTS\n";
00618         }
00619 
00620         if (simParameters->consForceOn)
00621           iout << iINFO << molecule->numConsForce << " CONSTANT FORCES\n";
00622 
00623         if (simParameters->stirOn)
00624           iout << iINFO << molecule->numStirredAtoms << " STIRRED ATOMS\n";
00625 
00626         if (simParameters->fixedAtomsOn)
00627         {
00628            iout << iINFO << molecule->numFixedAtoms << " FIXED ATOMS\n";
00629         }
00630 
00631         if (simParameters->rigidBonds)
00632         {
00633            iout << iINFO << molecule->numRigidBonds << " RIGID BONDS\n";
00634         }
00635 
00636         if (simParameters->fixedAtomsOn && simParameters->rigidBonds)
00637         {
00638            iout << iINFO << molecule->numFixedRigidBonds <<
00639                         " RIGID BONDS BETWEEN FIXED ATOMS\n";
00640         }
00641         
00642         /* BEGIN gf */
00643         if (simParameters->mgridforceOn)
00644         {
00645             int i;
00646             iout << iINFO << molecule->numGridforceGrids 
00647                  << " GRIDS ACTIVE\n";
00648         }
00649         /* END gf */
00650 
00651 //Modifications for alchemical fep
00652         if (simParameters->alchOn) {
00653           iout << iINFO << "ALCH: " 
00654                << molecule->numFepInitial <<
00655                " ATOMS TO DISAPPEAR IN FINAL STATE\n";
00656            iout << iINFO << "ALCH: " 
00657                <<  molecule->numFepFinal <<
00658                " ATOMS TO APPEAR IN FINAL STATE\n";
00659            if (molecule->suspiciousAlchBonds) {
00660              iout << iWARN << "ALCH: SUSPICIOUS BONDS BETWEEN INITIAL AND " <<
00661              "FINAL GROUPS WERE FOUND" << "\n" << endi;
00662            }
00663            if (molecule->alchDroppedAngles) {
00664              iout << iINFO << "ALCH: " 
00665                  << molecule->alchDroppedAngles <<
00666                  " ANGLES LINKING INITIAL AND FINAL ATOMS DELETED\n";
00667            }
00668            if (molecule->alchDroppedDihedrals) {
00669              iout << iINFO << "ALCH: "
00670                  << molecule->alchDroppedDihedrals <<
00671                  " DIHEDRALS LINKING INITIAL AND FINAL ATOMS DELETED\n";
00672            }
00673            if (molecule->alchDroppedImpropers) {
00674              iout << iINFO << "ALCH: "
00675                  << molecule->alchDroppedImpropers <<
00676                  " IMPROPERS LINKING INITIAL AND FINAL ATOMS DELETED\n";
00677            }
00678         }
00679 //fepe
00680 
00681         if (simParameters->lesOn) {
00682            iout << iINFO << molecule->numFepInitial <<
00683                " LOCALLY ENHANCED ATOMS ENABLED\n";
00684         }
00685        
00686         if (simParameters->pairInteractionOn) {
00687            iout << iINFO << "PAIR INTERACTION GROUP 1 CONTAINS "
00688                 <<  molecule->numFepInitial << " ATOMS\n";
00689            if (!simParameters->pairInteractionSelf) {
00690              iout << iINFO << "PAIR INTERACTION GROUP 2 CONTAINS "
00691                   <<  molecule->numFepFinal << " ATOMS\n";
00692            }
00693         }
00694            
00695 #if 1
00696         if (molecule->numLonepairs != 0) {
00697           iout << iINFO << molecule->numLonepairs << " LONE PAIRS\n";
00698         }
00699         if (molecule->numDrudeAtoms != 0) {
00700           iout << iINFO << molecule->numDrudeAtoms << " DRUDE ATOMS\n";
00701         }
00702         iout << iINFO << molecule->num_deg_freedom(1)
00703              << " DEGREES OF FREEDOM\n";
00704         if (simParameters->drudeOn) {
00705           int g_bond = 3 * molecule->numDrudeAtoms;
00706           int g_com = molecule->num_deg_freedom(1) - g_bond;
00707           iout << iINFO << g_com << " DRUDE COM DEGREES OF FREEDOM\n";
00708           iout << iINFO << g_bond << " DRUDE BOND DEGREES OF FREEDOM\n";
00709         }
00710 #endif
00711 #if 0
00712         {
00713           // Copied from Controller::printEnergies()
00714           int64 numAtoms = molecule->numAtoms;
00715           int64 numDegFreedom = 3 * numAtoms;
00716     int numLonepairs = molecule->numLonepairs;
00717           int numFixedAtoms = molecule->numFixedAtoms;
00718           if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
00719           if ( ! ( numFixedAtoms || molecule->numConstraints
00720                 || simParameters->comMove || simParameters->langevinOn ) ) {
00721             numDegFreedom -= 3;
00722           }
00723     if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
00724           int numRigidBonds = molecule->numRigidBonds;
00725           int numFixedRigidBonds = molecule->numFixedRigidBonds;
00726     // numLonepairs is subtracted here because all lonepairs have a rigid bond
00727     // to oxygen, but all of the LP degrees of freedom are dealt with above
00728     numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
00729           iout << iINFO << numDegFreedom << " DEGREES OF FREEDOM\n";
00730         }
00731 #endif
00732 
00733         iout << iINFO << molecule->numHydrogenGroups << " HYDROGEN GROUPS\n";
00734         iout << iINFO << molecule->maxHydrogenGroupSize
00735                 << " ATOMS IN LARGEST HYDROGEN GROUP\n";
00736         iout << iINFO << molecule->numMigrationGroups << " MIGRATION GROUPS\n";
00737         iout << iINFO << molecule->maxMigrationGroupSize
00738                 << " ATOMS IN LARGEST MIGRATION GROUP\n";
00739         if (simParameters->fixedAtomsOn)
00740         {
00741            iout << iINFO << molecule->numFixedGroups <<
00742                         " HYDROGEN GROUPS WITH ALL ATOMS FIXED\n";
00743         }
00744 
00745         {
00746           BigReal totalMass = 0;
00747           BigReal totalCharge = 0;
00748           int i;
00749           for ( i = 0; i < molecule->numAtoms; ++i ) {
00750             totalMass += molecule->atommass(i);
00751             totalCharge += molecule->atomcharge(i);
00752           }
00753           iout << iINFO << "TOTAL MASS = " << totalMass << " amu\n"; 
00754           iout << iINFO << "TOTAL CHARGE = " << totalCharge << " e\n"; 
00755 
00756           BigReal volume = lattice.volume();
00757           if ( volume ) {
00758             iout << iINFO << "MASS DENSITY = "
00759               << ((totalMass/volume) / 0.6022) << " g/cm^3\n";
00760             iout << iINFO << "ATOM DENSITY = "
00761               << (molecule->numAtoms/volume) << " atoms/A^3\n";
00762           }
00763         }
00764 
00765         iout << iINFO << "*****************************\n";
00766         iout << endi;
00767         fflush(stdout);
00768 
00769   StringList *binCoordinateFilename = configList->find("bincoordinates");
00770   if ( binCoordinateFilename && ! reload ) {
00771     read_binary_coors(binCoordinateFilename->data, pdb);
00772   }
00773 
00774   DebugM(4, "::configFileInit() - printing Molecule Information\n");
00775 
00776   molecule->print_atoms(parameters);
00777   molecule->print_bonds(parameters);
00778   molecule->print_exclusions();
00779   fflush(stdout);
00780 #endif
00781 
00782   DebugM(4, "::configFileInit() - done printing Molecule Information\n");
00783   DebugM(1, "::configFileInit() - done\n");
00784 
00785   return(0);
00786 }
00787 

Generated on Mon Nov 20 01:17:13 2017 for NAMD by  doxygen 1.4.7