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

Generated on Sat Jul 21 01:17:14 2018 for NAMD by  doxygen 1.4.7