00001
00002
00003
00004
00005
00006
00007
00008
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
00028 #include "Debug.h"
00029
00030 #include "CompressPsf.h"
00031 #include "PluginIOMgr.h"
00032
00033 NamdState::NamdState()
00034 {
00035 configList = NULL;
00036 simParameters = NULL;
00037 parameters = NULL;
00038 molecule = NULL;
00039 pdb = NULL;
00040 }
00041
00042 int
00043 NamdState::status()
00044 {
00045 int ret=0;
00046 if (configList != NULL) {
00047 DebugM(1, "Config List exists\n");
00048 } else ret++;
00049
00050 if (simParameters != NULL) {
00051 DebugM(1, "SimParameters exists\n");
00052 }
00053 else ret++;
00054
00055 if (parameters != NULL) {
00056 DebugM(1,"Parameters exists\n");
00057 }
00058 else ret++;
00059
00060 if (molecule != NULL) {
00061 DebugM(1, "Molecule exists\n");
00062 }
00063 else ret++;
00064
00065 if (pdb != NULL) {
00066 DebugM(1,"PDB exists \n");
00067 }
00068 else ret++;
00069
00070 return(ret);
00071 }
00072
00073 void NamdState:: useController(Controller *controllerPtr)
00074 {
00075 controller=controllerPtr;
00076 }
00077
00078 void NamdState:: runController(void)
00079 {
00080 controller->run();
00081 }
00082
00083 extern void read_binary_coors(char *fname, PDB *pdbobj);
00084
00085 #ifdef MEM_OPT_VERSION
00086
00087
00088 void NamdState::checkMemOptCompatibility(){
00089 if(simParameters->genCompressedPsf) {
00090 iout << "In the memory optimized version, the compression of molecule information is not supported! "
00091 << "Please use the non-memory optimized version.\n" <<endi;
00092 NAMD_die("MEMOPT: unsupported usage");
00093 }
00094
00095 if(simParameters->langevinOn && simParameters->langevinDamping == 0.0) {
00096 iout << iWARN << "langevinDamping MUST NOT BE 0.0 IF LANGEVIN IS"
00097 <<" TURNED ON IN MEMORY OPTIMIZED VERSION\n" <<endi;
00098 NAMD_die("MEMOPT: unsupported feature");
00099 }
00100 if(simParameters->tCoupleOn)
00101 NAMD_die("MEMOPT: tCouple is not supported in memory optimized version");
00102 if(simParameters->pairInteractionOn)
00103 NAMD_die("MEMOPT: pairInteractionOn could not be enabled in memory optimized version");
00104 if(simParameters->alchFepOn || simParameters->alchOn){
00105 iout << iWARN << "ALCH: AUTOMATIC DELETION OF BONDED INTERACTIONS "
00106 << "BETWEEN INITIAL AND FINAL GROUPS IS NOT SUPPORTED IN MEMORY "
00107 << "OPTIMISED VERSION - MANUAL PROCESSING IS NECESSARY\n" << endi;
00108 NAMD_die("MEMOPT: unsupported feature");
00109 }
00110 if(simParameters->alchThermIntOn)
00111 NAMD_die("MEMOPT: alchThermIntOn could not be enabled in memory optimized version");
00112 if(simParameters->lesOn)
00113 NAMD_die("MEMOPT: lesOn could not be enabled in memory optimized version");
00114 if(simParameters->lonepairs) {
00115 NAMD_die("MEMOPT: lonepairs could not be enabled in memory optimized version");
00116 }
00117 }
00118 #endif
00119
00120 int NamdState::configListInit(ConfigList *cfgList) {
00121 configList = cfgList;
00122 if (!configList->okay()) {
00123 NAMD_die("Simulation config file is incomplete or contains errors.");
00124 }
00125 DebugM(1,"NamdState::configFileInit configList okay\n");
00126
00127 char *currentdir = 0;
00128 simParameters = new SimParameters(configList,currentdir);
00129 fflush(stdout);
00130 lattice = simParameters->lattice;
00131
00132
00133
00134 #ifdef MEM_OPT_VERSION
00135 checkMemOptCompatibility();
00136 #endif
00137
00138
00139 if(simParameters->genCompressedPsf) {
00140 if(simParameters->rigidBonds == RIGID_NONE){
00141
00142 simParameters->rigidBonds = RIGID_ALL;
00143 }
00144 }
00145
00146 StringList *molInfoFilename;
00147
00148
00149
00150
00151 if (simParameters->amberOn)
00152 { StringList *parmFilename = configList->find("parmfile");
00153 molInfoFilename = parmFilename;
00154 StringList *coorFilename = configList->find("ambercoor");
00155
00156
00157
00158 Ambertoppar *amber;
00159 amber = new Ambertoppar;
00160 if (amber->readparm(parmFilename->data))
00161 { parameters = new Parameters(amber, simParameters->vdwscale14);
00162 molecule = new Molecule(simParameters, parameters, amber);
00163 if (coorFilename != NULL)
00164 pdb = new PDB(coorFilename->data,amber);
00165 delete amber;
00166 }
00167 else
00168 NAMD_die("Failed to read AMBER parm file!");
00169 parameters->print_param_summary();
00170 }
00171 else if (simParameters->gromacsOn) {
00172 StringList *topFilename = configList->find("grotopfile");
00173 molInfoFilename = topFilename;
00174 StringList *coorFilename = configList->find("grocoorfile");
00175
00176
00177
00178 GromacsTopFile *gromacsFile;
00179 gromacsFile = new GromacsTopFile(topFilename->data);
00180 parameters = new Parameters(gromacsFile,simParameters->minimizeCGOn);
00181 if (coorFilename != NULL)
00182 pdb = new PDB(coorFilename->data,gromacsFile);
00183
00184 molecule = new Molecule(simParameters, parameters, gromacsFile);
00185
00186
00187 delete gromacsFile;
00188
00189
00190
00191
00192
00193 parameters->print_param_summary();
00194 }
00195 else if (simParameters->usePluginIO){
00196 #ifdef MEM_OPT_VERSION
00197 NAMD_die("Using plugin IO is not supported in memory optimized version!");
00198 #else
00199 PluginIOMgr *pIOMgr = new PluginIOMgr();
00200
00201 iout << iWARN << "Plugin-based I/O is still in development and may still have bugs\n" << endi;
00202
00203 molfile_plugin_t *pIOHandle = pIOMgr->getPlugin();
00204 if (pIOHandle == NULL) {
00205 NAMD_die("ERROR: Failed to match requested plugin type");
00206 }
00207 if ( pIOHandle->open_file_read == NULL )
00208 NAMD_die("ERROR: Selected plugin type cannot open files");
00209 if ( pIOHandle->read_structure == NULL )
00210 NAMD_die("ERROR: Selected plugin type cannot read structures");
00211 if ( pIOHandle->read_next_timestep == NULL )
00212 NAMD_die("ERROR: Selected plugin type cannot read coordinates");
00213
00214 StringList *moleculeFilename = configList->find("structure");
00215 molInfoFilename = moleculeFilename;
00216 StringList *parameterFilename = configList->find("parameters");
00217
00218
00219 parameters = new Parameters(simParameters, parameterFilename);
00220 parameters->print_param_summary();
00221
00222 int numAtoms = 0;
00223
00224 void *plgFile = pIOHandle->open_file_read(moleculeFilename->data,
00225 pIOHandle->name, &numAtoms);
00226 if(plgFile == NULL) {
00227 NAMD_die("ERROR: Opening structure file failed!");
00228 }
00229
00230 double fileReadTime = CmiWallTimer();
00231 molecule = new Molecule(simParameters, parameters, pIOHandle, plgFile, numAtoms);
00232 iout << iINFO << "TIME FOR LOAD MOLECULE STRUCTURE INFORMATION: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00233
00234
00235 if(!simParameters->genCompressedPsf) {
00236 fileReadTime = CmiWallTimer();
00237
00238
00239 pdb = new PDB(pIOHandle, plgFile, molecule->numAtoms, molecule->getOccupancyData(), molecule->getBFactorData());
00240 molecule->freeOccupancyData();
00241 molecule->freeBFactorData();
00242 iout << iINFO << "TIME FOR LOADING ATOMS' COORDINATES INFORMATION: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00243 }
00244
00245 pIOHandle->close_file_read(plgFile);
00246 delete pIOMgr;
00247 #endif
00248 }
00249 else
00250 {
00251 StringList *moleculeFilename = configList->find("structure");
00252 molInfoFilename = moleculeFilename;
00253 StringList *parameterFilename = configList->find("parameters");
00254
00255
00256 parameters = new Parameters(simParameters, parameterFilename);
00257
00258
00259 parameters->print_param_summary();
00260
00261 double fileReadTime = CmiWallTimer();
00262 molecule = new Molecule(simParameters, parameters, moleculeFilename->data, configList);
00263 iout << iINFO << "TIME FOR READING PSF FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00264 }
00265
00266 fflush(stdout);
00267
00268 #ifdef MEM_OPT_VERSION
00269
00270
00271 if(simParameters->numinputprocs==0){
00272 int numatoms = molecule->numAtoms;
00273 long estval = (sizeof(InputAtom)+2*sizeof(int)+1)*((long)(numatoms));
00274 int numprocs = estval>>26;
00275 if(numprocs==0){
00276 numprocs=1;
00277 }else if(numprocs>CkNumPes()){
00278 numprocs=CkNumPes();
00279 }
00280 simParameters->numinputprocs=numprocs;
00281 }
00282 if(simParameters->numoutputprocs==0){
00283 int numatoms = molecule->numAtoms;
00284 long estval = (sizeof(Vector)*2)*((long)(numatoms));
00285 int numprocs = estval>>26;
00286 if(numprocs==0){
00287 numprocs=1;
00288 }else if(numprocs>CkNumPes()){
00289 numprocs=CkNumPes();
00290 }
00291 simParameters->numoutputprocs=numprocs;
00292 }
00293
00294 if(simParameters->numoutputwrts > simParameters->numoutputprocs) {
00295 simParameters->numoutputwrts = simParameters->numoutputprocs;
00296 }
00297
00298 if (simParameters->fixedAtomsOn){
00299 double fileReadTime = CmiWallTimer();
00300 molecule->load_fixed_atoms(configList->find("fixedAtomListFile"));
00301 iout << iINFO << "TIME FOR READING FIXED ATOMS FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00302 }
00303 #else
00304 if (simParameters->extraBondsOn) {
00305
00306
00307
00308 if(!simParameters->useCompressedPsf)
00309 molecule->build_extra_bonds(parameters, configList->find("extraBondsFile"));
00310 }
00311 if(simParameters->genCompressedPsf) {
00312 double fileReadTime = CmiWallTimer();
00313 compress_molecule_info(molecule, molInfoFilename->data, parameters, simParameters, configList);
00314 iout << "Finished compressing molecule information, which takes " << CmiWallTimer()-fileReadTime <<"(s)\n"<<endi;
00315 CkExit();
00316 }
00317
00318
00319 StringList *coordinateFilename = NULL;
00320 if(!simParameters->usePluginIO) {
00321
00322
00323
00324 coordinateFilename = configList->find("coordinates");
00325 double fileReadTime = CmiWallTimer();
00326 if (coordinateFilename != NULL)
00327 pdb = new PDB(coordinateFilename->data);
00328 if (pdb->num_atoms() != molecule->numAtoms) {
00329 NAMD_die("Number of pdb and psf atoms are not the same!");
00330 }
00331 iout << iINFO << "TIME FOR READING PDB FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00332 iout << iINFO << "\n" << endi;
00333 }
00334
00335
00336 if (simParameters->constraintsOn)
00337 {
00338 StringList *consRefFile = configList->find("consref");
00339 StringList *consKFile = configList->find("conskfile");
00340
00341 if (coordinateFilename != NULL) {
00342 if(strcasecmp(coordinateFilename->data, consRefFile->data)==0)
00343 consRefFile = NULL;
00344 if(strcasecmp(coordinateFilename->data, consKFile->data)==0)
00345 consKFile = NULL;
00346 }
00347
00348 molecule->build_constraint_params(consRefFile, consKFile,
00349 configList->find("conskcol"),
00350 pdb,
00351 NULL);
00352 }
00353
00354
00355 if (simParameters->stirOn)
00356 {
00357
00358
00359 molecule->build_stirred_atoms(configList->find("stirFilename"),
00360 configList->find("stirredAtomsCol"),
00361 pdb,
00362 NULL);
00363 }
00364
00365
00366 if (simParameters->fixedAtomsOn)
00367 {
00368 molecule->build_fixed_atoms(configList->find("fixedatomsfile"),
00369 configList->find("fixedatomscol"),
00370 pdb,
00371 NULL);
00372 }
00373
00374
00375 if (simParameters->mgridforceOn)
00376 {
00377 molecule->build_gridforce_params(configList->find("gridforcefile"),
00378 configList->find("gridforcecol"),
00379 configList->find("gridforcechargecol"),
00380 configList->find("gridforcepotfile"),
00381 pdb,
00382 NULL);
00383 }
00384
00385
00386
00387 if (simParameters->consForceOn) {
00388 char *filename = NULL;
00389 if (configList->find("consforcefile"))
00390 filename = configList->find("consforcefile")->data;
00391 molecule->build_constant_forces(filename);
00392 }
00393
00394 if (simParameters->excludeFromPressure) {
00395 molecule->build_exPressure_atoms(
00396 configList->find("excludeFromPressureFile"),
00397 configList->find("excludeFromPressureCol"),
00398 pdb, NULL);
00399 }
00400
00401
00402 if (simParameters->movDragOn) {
00403 molecule->build_movdrag_params(configList->find("movDragFile"),
00404 configList->find("movDragCol"),
00405 configList->find("movDragVelFile"),
00406 pdb,
00407 NULL);
00408 }
00409
00410
00411 if (simParameters->rotDragOn) {
00412 molecule->build_rotdrag_params(configList->find("rotDragFile"),
00413 configList->find("rotDragCol"),
00414 configList->find("rotDragAxisFile"),
00415 configList->find("rotDragPivotFile"),
00416 configList->find("rotDragVelFile"),
00417 configList->find("rotDragVelCol"),
00418 pdb,
00419 NULL);
00420 }
00421
00422
00423 if (simParameters->consTorqueOn) {
00424 molecule->build_constorque_params(configList->find("consTorqueFile"),
00425 configList->find("consTorqueCol"),
00426 configList->find("consTorqueAxisFile"),
00427 configList->find("consTorquePivotFile"),
00428 configList->find("consTorqueValFile"),
00429 configList->find("consTorqueValCol"),
00430 pdb,
00431 NULL);
00432 }
00433
00434
00435
00436 if (simParameters->langevinOn)
00437 {
00438 if (simParameters->langevinDamping == 0.0) {
00439 molecule->build_langevin_params(configList->find("langevinfile"),
00440 configList->find("langevincol"),
00441 pdb,
00442 NULL);
00443 } else {
00444 molecule->build_langevin_params(simParameters->langevinDamping,
00445 simParameters->drudeDamping,
00446 simParameters->langevinHydrogen);
00447 }
00448 }
00449 else if (simParameters->tCoupleOn)
00450 {
00451
00452
00453 molecule->build_langevin_params(configList->find("tcouplefile"),
00454 configList->find("tcouplecol"),
00455 pdb,
00456 NULL);
00457 }
00458
00459
00460
00461 if (simParameters->alchOn) {
00462 molecule->build_fep_flags(configList->find("alchfile"),
00463 configList->find("alchcol"), pdb, NULL, "alch" );
00464 molecule->delete_alch_bonded();
00465 }
00466
00467 if (simParameters->lesOn) {
00468 if (simParameters->alchFepOn || simParameters->alchThermIntOn) NAMD_bug("FEP/TI and LES are incompatible!");
00469 molecule->build_fep_flags(configList->find("lesfile"),
00470 configList->find("lescol"), pdb, NULL, "les");
00471 }
00472 if (simParameters->pairInteractionOn) {
00473 molecule->build_fep_flags(configList->find("pairInteractionFile"),
00474 configList->find("pairInteractionCol"), pdb, NULL, "pairInteraction");
00475 }
00476 if (simParameters->pressureProfileAtomTypes > 1) {
00477 molecule->build_fep_flags(configList->find("pressureProfileAtomTypesFile"),
00478 configList->find("pressureProfileAtomTypesCol"), pdb, NULL, "pressureProfileAtomTypes");
00479 }
00480
00481 #ifdef OPENATOM_VERSION
00482 if (simParameters->openatomOn) {
00483 molecules->build_qmmm_flags(configList->find("openatomPdbFile",
00484 configList->find("openatomPdbCol"), pdb, NULL, "openatomPdb")
00485 }
00486 #endif // OPENATOM_VERSION
00487
00488
00489 if (simParameters->goForcesOn) {
00490 #ifdef MEM_OPT_VERSION
00491 NAMD_die("Go forces are not supported in memory-optimized builds.");
00492 #else
00493 StringList *moleculeFilename = configList->find("structure");
00494 StringList *parameterFilename = configList->find("parameters");
00495 StringList *goFilename = configList->find("goParameters");
00496 StringList *goStructureFilename = configList->find("goCoordinates");
00497
00498
00499 molecule->build_go_params(goFilename);
00500
00501 int goMethod = simParameters->goMethod;
00502 if (goMethod == 1) {
00503 iout << iINFO << "Using Go method matrix\n" << endi;
00504 molecule->build_go_sigmas(goStructureFilename, NULL);
00505 } else if (goMethod == 2) {
00506 iout << iINFO << "Using Go method sparse matrix\n" << endi;
00507 NAMD_die("Go method sparse not yet implemented. Crashing");
00508
00509 } else if (goMethod == 3) {
00510 iout << iINFO << "Using Go method lowmem\n" << endi;
00511 molecule->build_go_arrays(goStructureFilename, NULL);
00512 } else {
00513 NAMD_die("Failed to read goMethod variable in NamdState.C");
00514 }
00515 #endif
00516 }
00517
00518
00519 iout << iINFO << "****************************\n";
00520 iout << iINFO << "STRUCTURE SUMMARY:\n";
00521 iout << iINFO << molecule->numAtoms << " ATOMS\n";
00522 iout << iINFO << molecule->numBonds << " BONDS\n";
00523 iout << iINFO << molecule->numAngles << " ANGLES\n";
00524 iout << iINFO << molecule->numDihedrals << " DIHEDRALS\n";
00525 iout << iINFO << molecule->numImpropers << " IMPROPERS\n";
00526 iout << iINFO << molecule->numCrossterms << " CROSSTERMS\n";
00527 iout << iINFO << molecule->numExclusions << " EXCLUSIONS\n";
00528
00529
00530 if ((molecule->numMultipleDihedrals) && (simParameters->paraTypeXplorOn))
00531 {
00532 iout << iINFO << molecule->numMultipleDihedrals
00533 << " DIHEDRALS WITH MULTIPLE PERIODICITY (BASED ON PSF FILE)\n";
00534 }
00535 if ((molecule->numMultipleDihedrals) && (simParameters->paraTypeCharmmOn))
00536 {
00537 iout << iINFO << molecule->numMultipleDihedrals
00538 << " DIHEDRALS WITH MULTIPLE PERIODICITY IGNORED (BASED ON PSF FILE) \n";
00539 iout << iINFO
00540 << " CHARMM MULTIPLICITIES BASED ON PARAMETER FILE INFO! \n";
00541 }
00542
00543
00544 if (molecule->numMultipleImpropers)
00545 {
00546 iout << iINFO << molecule->numMultipleImpropers
00547 << " IMPROPERS WITH MULTIPLE PERIODICITY\n";
00548 }
00549
00550 if (simParameters->constraintsOn)
00551 {
00552 iout << iINFO << molecule->numConstraints << " CONSTRAINTS\n";
00553 }
00554
00555 if (simParameters->consForceOn)
00556 iout << iINFO << molecule->numConsForce << " CONSTANT FORCES\n";
00557
00558 if (simParameters->stirOn)
00559 iout << iINFO << molecule->numStirredAtoms << " STIRRED ATOMS\n";
00560
00561 if (simParameters->fixedAtomsOn)
00562 {
00563 iout << iINFO << molecule->numFixedAtoms << " FIXED ATOMS\n";
00564 }
00565
00566 if (simParameters->rigidBonds)
00567 {
00568 iout << iINFO << molecule->numRigidBonds << " RIGID BONDS\n";
00569 }
00570
00571 if (simParameters->fixedAtomsOn && simParameters->rigidBonds)
00572 {
00573 iout << iINFO << molecule->numFixedRigidBonds <<
00574 " RIGID BONDS BETWEEN FIXED ATOMS\n";
00575 }
00576
00577
00578 if (simParameters->mgridforceOn)
00579 {
00580 int i;
00581 iout << iINFO << molecule->numGridforceGrids
00582 << " GRIDS ACTIVE\n";
00583 }
00584
00585
00586
00587 if (simParameters->alchFepOn || simParameters->alchThermIntOn) {
00588 iout << iINFO << "ALCH: "
00589 << molecule->numFepInitial <<
00590 " ATOMS TO DISAPPEAR IN FINAL STATE\n";
00591 iout << iINFO << "ALCH: "
00592 << molecule->numFepFinal <<
00593 " ATOMS TO APPEAR IN FINAL STATE\n";
00594 if (molecule->suspiciousAlchBonds) {
00595 iout << iWARN << "ALCH: SUSPICIOUS BONDS BETWEEN INITIAL AND " <<
00596 "FINAL GROUPS WERE FOUND" << "\n" << endi;
00597 }
00598 if (molecule->alchDroppedAngles) {
00599 iout << iINFO << "ALCH: "
00600 << molecule->alchDroppedAngles <<
00601 " ANGLES LINKING INITIAL AND FINAL ATOMS DELETED\n";
00602 }
00603 if (molecule->alchDroppedDihedrals) {
00604 iout << iINFO << "ALCH: "
00605 << molecule->alchDroppedDihedrals <<
00606 " DIHEDRALS LINKING INITIAL AND FINAL ATOMS DELETED\n";
00607 }
00608 if (molecule->alchDroppedImpropers) {
00609 iout << iINFO << "ALCH: "
00610 << molecule->alchDroppedImpropers <<
00611 " IMPROPERS LINKING INITIAL AND FINAL ATOMS DELETED\n";
00612 }
00613 }
00614
00615
00616 if (simParameters->lesOn) {
00617 iout << iINFO << molecule->numFepInitial <<
00618 " LOCALLY ENHANCED ATOMS ENABLED\n";
00619 }
00620
00621 if (simParameters->pairInteractionOn) {
00622 iout << iINFO << "PAIR INTERACTION GROUP 1 CONTAINS "
00623 << molecule->numFepInitial << " ATOMS\n";
00624 if (!simParameters->pairInteractionSelf) {
00625 iout << iINFO << "PAIR INTERACTION GROUP 2 CONTAINS "
00626 << molecule->numFepFinal << " ATOMS\n";
00627 }
00628 }
00629
00630 #if 1
00631 if (molecule->numLonepairs != 0) {
00632 iout << iINFO << molecule->numLonepairs << " LONE PAIRS\n";
00633 }
00634 if (molecule->numDrudeAtoms != 0) {
00635 iout << iINFO << molecule->numDrudeAtoms << " DRUDE ATOMS\n";
00636 }
00637 iout << iINFO << molecule->num_deg_freedom(1)
00638 << " DEGREES OF FREEDOM\n";
00639 if (simParameters->drudeOn) {
00640 int g_bond = 3 * molecule->numDrudeAtoms;
00641 int g_com = molecule->num_deg_freedom(1) - g_bond;
00642 iout << iINFO << g_com << " DRUDE COM DEGREES OF FREEDOM\n";
00643 iout << iINFO << g_bond << " DRUDE BOND DEGREES OF FREEDOM\n";
00644 }
00645 #endif
00646 #if 0
00647 {
00648
00649 int numAtoms = molecule->numAtoms;
00650 int numDegFreedom = 3 * numAtoms;
00651 int numLonepairs = molecule->numLonepairs;
00652 int numFixedAtoms = molecule->numFixedAtoms;
00653 if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
00654 if ( ! ( numFixedAtoms || molecule->numConstraints
00655 || simParameters->comMove || simParameters->langevinOn ) ) {
00656 numDegFreedom -= 3;
00657 }
00658 if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
00659 int numRigidBonds = molecule->numRigidBonds;
00660 int numFixedRigidBonds = molecule->numFixedRigidBonds;
00661
00662
00663 numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
00664 iout << iINFO << numDegFreedom << " DEGREES OF FREEDOM\n";
00665 }
00666 #endif
00667
00668 iout << iINFO << molecule->numHydrogenGroups << " HYDROGEN GROUPS\n";
00669 iout << iINFO << molecule->maxHydrogenGroupSize
00670 << " ATOMS IN LARGEST HYDROGEN GROUP\n";
00671 iout << iINFO << molecule->numMigrationGroups << " MIGRATION GROUPS\n";
00672 iout << iINFO << molecule->maxMigrationGroupSize
00673 << " ATOMS IN LARGEST MIGRATION GROUP\n";
00674 if (simParameters->fixedAtomsOn)
00675 {
00676 iout << iINFO << molecule->numFixedGroups <<
00677 " HYDROGEN GROUPS WITH ALL ATOMS FIXED\n";
00678 }
00679
00680 {
00681 BigReal totalMass = 0;
00682 BigReal totalCharge = 0;
00683 int i;
00684 for ( i = 0; i < molecule->numAtoms; ++i ) {
00685 totalMass += molecule->atommass(i);
00686 totalCharge += molecule->atomcharge(i);
00687 }
00688 iout << iINFO << "TOTAL MASS = " << totalMass << " amu\n";
00689 iout << iINFO << "TOTAL CHARGE = " << totalCharge << " e\n";
00690
00691 BigReal volume = lattice.volume();
00692 if ( volume ) {
00693 iout << iINFO << "MASS DENSITY = "
00694 << ((totalMass/volume) / 0.6022) << " g/cm^3\n";
00695 iout << iINFO << "ATOM DENSITY = "
00696 << (molecule->numAtoms/volume) << " atoms/A^3\n";
00697 }
00698 }
00699
00700 iout << iINFO << "*****************************\n";
00701 iout << endi;
00702 fflush(stdout);
00703
00704 StringList *binCoordinateFilename = configList->find("bincoordinates");
00705 if ( binCoordinateFilename ) {
00706 read_binary_coors(binCoordinateFilename->data, pdb);
00707 }
00708
00709 DebugM(4, "::configFileInit() - printing Molecule Information\n");
00710
00711 molecule->print_atoms(parameters);
00712 molecule->print_bonds(parameters);
00713 molecule->print_exclusions();
00714 fflush(stdout);
00715 #endif
00716
00717 DebugM(4, "::configFileInit() - done printing Molecule Information\n");
00718 DebugM(1, "::configFileInit() - done\n");
00719
00720 return(0);
00721 }
00722