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 NamdState::NamdState()
00031 {
00032 configList = NULL;
00033 simParameters = NULL;
00034 parameters = NULL;
00035 molecule = NULL;
00036 pdb = NULL;
00037 }
00038
00039 int
00040 NamdState::status()
00041 {
00042 int ret=0;
00043 if (configList != NULL) {
00044 DebugM(1, "Config List exists\n");
00045 } else ret++;
00046
00047 if (simParameters != NULL) {
00048 DebugM(1, "SimParameters exists\n");
00049 }
00050 else ret++;
00051
00052 if (parameters != NULL) {
00053 DebugM(1,"Parameters exists\n");
00054 }
00055 else ret++;
00056
00057 if (molecule != NULL) {
00058 DebugM(1, "Molecule exists\n");
00059 }
00060 else ret++;
00061
00062 if (pdb != NULL) {
00063 DebugM(1,"PDB exists \n");
00064 }
00065 else ret++;
00066
00067 return(ret);
00068 }
00069
00070 void NamdState:: useController(Controller *controllerPtr)
00071 {
00072 controller=controllerPtr;
00073 }
00074
00075 void NamdState:: runController(void)
00076 {
00077 controller->run();
00078 }
00079
00080 extern void read_binary_coors(char *fname, PDB *pdbobj);
00081
00082 int NamdState::configListInit(ConfigList *cfgList) {
00083 configList = cfgList;
00084 if (!configList->okay()) {
00085 NAMD_die("Simulation config file is incomplete or contains errors.");
00086 }
00087 DebugM(1,"NamdState::configFileInit configList okay\n");
00088
00089 char *currentdir = 0;
00090 simParameters = new SimParameters(configList,currentdir);
00091 fflush(stdout);
00092 lattice = simParameters->lattice;
00093
00094
00095
00096
00097
00098 if (simParameters->amberOn)
00099 { StringList *parmFilename = configList->find("parmfile");
00100 StringList *coorFilename = configList->find("ambercoor");
00101
00102
00103
00104 Ambertoppar *amber;
00105 amber = new Ambertoppar;
00106 if (amber->readparm(parmFilename->data))
00107 { parameters = new Parameters(amber, simParameters->vdwscale14);
00108 molecule = new Molecule(simParameters, parameters, amber);
00109 if (coorFilename != NULL)
00110 pdb = new PDB(coorFilename->data,amber);
00111 delete amber;
00112 }
00113 else
00114 NAMD_die("Failed to read AMBER parm file!");
00115 parameters->print_param_summary();
00116 }
00117 else if (simParameters->gromacsOn) {
00118 StringList *topFilename = configList->find("grotopfile");
00119 StringList *coorFilename = configList->find("grocoorfile");
00120
00121
00122
00123 GromacsTopFile *gromacsFile;
00124 gromacsFile = new GromacsTopFile(topFilename->data);
00125 parameters = new Parameters(gromacsFile,simParameters->minimizeCGOn);
00126 if (coorFilename != NULL)
00127 pdb = new PDB(coorFilename->data,gromacsFile);
00128
00129 molecule = new Molecule(simParameters, parameters, gromacsFile);
00130
00131
00132 delete gromacsFile;
00133
00134
00135
00136
00137
00138 parameters->print_param_summary();
00139 }
00140 else
00141 { StringList *moleculeFilename = configList->find("structure");
00142 StringList *parameterFilename = configList->find("parameters");
00143
00144
00145 parameters = new Parameters(simParameters, parameterFilename);
00146
00147
00148 if(simParameters->genCompressedPsf){
00149 molecule = new Molecule(simParameters, parameters, moleculeFilename->data, configList);
00150 iout << "Finished Compressing psf file!\n" <<endi;
00151 CkExit();
00152 }
00153 else{
00154 parameters->print_param_summary();
00155 double fileReadTime = CmiWallTimer();
00156 molecule = new Molecule(simParameters, parameters, moleculeFilename->data);
00157 iout << iINFO << "TIME FOR READING PSF FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00158 }
00159 }
00160 fflush(stdout);
00161
00162 StringList *coordinateFilename = configList->find("coordinates");
00163
00164 double fileReadTime = CmiWallTimer();
00165 #ifdef MEM_OPT_VERSION
00166 if (coordinateFilename != NULL)
00167 pdb = new PDB(coordinateFilename->data, molecule->numAtoms);
00168 #else
00169 if (coordinateFilename != NULL)
00170 pdb = new PDB(coordinateFilename->data);
00171 if (pdb->num_atoms() != molecule->numAtoms) {
00172 NAMD_die("Number of pdb and psf atoms are not the same!");
00173 }
00174 #endif
00175 iout << iINFO << "TIME FOR READING PDB FILE: " << CmiWallTimer() - fileReadTime << "\n" << endi;
00176 iout << iINFO << "\n" << endi;
00177
00178 StringList *binCoordinateFilename = configList->find("bincoordinates");
00179 if ( binCoordinateFilename ) {
00180 read_binary_coors(binCoordinateFilename->data, pdb);
00181 }
00182
00183
00184 if (simParameters->constraintsOn)
00185 {
00186 StringList *consRefFile = configList->find("consref");
00187 StringList *consKFile = configList->find("conskfile");
00188
00189 if (coordinateFilename != NULL) {
00190 if(strcasecmp(coordinateFilename->data, consRefFile->data)==0)
00191 consRefFile = NULL;
00192 if(strcasecmp(coordinateFilename->data, consKFile->data)==0)
00193 consKFile = NULL;
00194 }
00195
00196 molecule->build_constraint_params(consRefFile, consKFile,
00197 configList->find("conskcol"),
00198 pdb,
00199 NULL);
00200 }
00201
00202 if (simParameters->stirOn)
00203 {
00204
00205
00206 molecule->build_stirred_atoms(configList->find("stirFilename"),
00207 configList->find("stirredAtomsCol"),
00208 pdb,
00209 NULL);
00210 }
00211
00212 if (simParameters->extraBondsOn) {
00213 #ifdef MEM_OPT_VERSION
00214 iout << "Information about the extra bonds have been already integrated into the compressed psf file!\n" <<endi;
00215 #else
00216 molecule->build_extra_bonds(parameters,
00217 configList->find("extraBondsFile"));
00218 #endif
00219 }
00220
00221 if (simParameters->fixedAtomsOn)
00222 {
00223 molecule->build_fixed_atoms(configList->find("fixedatomsfile"),
00224 configList->find("fixedatomscol"),
00225 pdb,
00226 NULL);
00227 }
00228
00229
00230 if (simParameters->gridforceOn)
00231 {
00232 molecule->build_gridforce_params(configList->find("gridforcefile"),
00233 configList->find("gridforcecol"),
00234 configList->find("gridforceqcol"),
00235 configList->find("gridforcevfile"),
00236 pdb,
00237 NULL);
00238 }
00239
00240
00241
00242 if (simParameters->consForceOn) {
00243 char *filename = NULL;
00244 if (configList->find("consforcefile"))
00245 filename = configList->find("consforcefile")->data;
00246 molecule->build_constant_forces(filename);
00247 }
00248
00249 if (simParameters->excludeFromPressure) {
00250 molecule->build_exPressure_atoms(
00251 configList->find("excludeFromPressureFile"),
00252 configList->find("excludeFromPressureCol"),
00253 pdb, NULL);
00254 }
00255
00256
00257 if (simParameters->movDragOn) {
00258 molecule->build_movdrag_params(configList->find("movDragFile"),
00259 configList->find("movDragCol"),
00260 configList->find("movDragVelFile"),
00261 pdb,
00262 NULL);
00263 }
00264
00265
00266 if (simParameters->rotDragOn) {
00267 molecule->build_rotdrag_params(configList->find("rotDragFile"),
00268 configList->find("rotDragCol"),
00269 configList->find("rotDragAxisFile"),
00270 configList->find("rotDragPivotFile"),
00271 configList->find("rotDragVelFile"),
00272 configList->find("rotDragVelCol"),
00273 pdb,
00274 NULL);
00275 }
00276
00277
00278 if (simParameters->consTorqueOn) {
00279 molecule->build_constorque_params(configList->find("consTorqueFile"),
00280 configList->find("consTorqueCol"),
00281 configList->find("consTorqueAxisFile"),
00282 configList->find("consTorquePivotFile"),
00283 configList->find("consTorqueValFile"),
00284 configList->find("consTorqueValCol"),
00285 pdb,
00286 NULL);
00287 }
00288
00289
00290
00291 if (simParameters->langevinOn)
00292 {
00293 if (simParameters->langevinDamping == 0.0) {
00294 molecule->build_langevin_params(configList->find("langevinfile"),
00295 configList->find("langevincol"),
00296 pdb,
00297 NULL);
00298 } else {
00299 molecule->build_langevin_params(simParameters->langevinDamping,
00300 simParameters->langevinHydrogen);
00301 }
00302 }
00303 else if (simParameters->tCoupleOn)
00304 {
00305
00306
00307 molecule->build_langevin_params(configList->find("tcouplefile"),
00308 configList->find("tcouplecol"),
00309 pdb,
00310 NULL);
00311 }
00312
00313
00314
00315
00316 if (simParameters->fepOn) {
00317 molecule->build_fep_flags(configList->find("fepfile"),
00318 configList->find("fepcol"), pdb, NULL);
00319 }
00320 if (simParameters->thermInt) {
00321 molecule->build_fep_flags(configList->find("tifile"),
00322 configList->find("ticol"), pdb, NULL);
00323 }
00324
00325 if (simParameters->lesOn) {
00326 if (simParameters->fepOn || simParameters->thermInt) NAMD_bug("FEP/TI and LES are incompatible!");
00327 molecule->build_fep_flags(configList->find("lesfile"),
00328 configList->find("lescol"), pdb, NULL);
00329 }
00330 if (simParameters->pairInteractionOn) {
00331 molecule->build_fep_flags(configList->find("pairInteractionFile"),
00332 configList->find("pairInteractionCol"), pdb, NULL);
00333 }
00334 if (simParameters->pressureProfileAtomTypes > 1) {
00335 molecule->build_fep_flags(configList->find("pressureProfileAtomTypesFile"),
00336 configList->find("pressureProfileAtomTypesCol"), pdb, NULL);
00337 }
00338
00339 iout << iINFO << "****************************\n";
00340 iout << iINFO << "STRUCTURE SUMMARY:\n";
00341 iout << iINFO << molecule->numAtoms << " ATOMS\n";
00342 iout << iINFO << molecule->numBonds << " BONDS\n";
00343 iout << iINFO << molecule->numAngles << " ANGLES\n";
00344 iout << iINFO << molecule->numDihedrals << " DIHEDRALS\n";
00345 iout << iINFO << molecule->numImpropers << " IMPROPERS\n";
00346 iout << iINFO << molecule->numCrossterms << " CROSSTERMS\n";
00347 iout << iINFO << molecule->numExclusions << " EXCLUSIONS\n";
00348
00349
00350 if ((molecule->numMultipleDihedrals) && (simParameters->paraTypeXplorOn))
00351 {
00352 iout << iINFO << molecule->numMultipleDihedrals
00353 << " DIHEDRALS WITH MULTIPLE PERIODICITY (BASED ON PSF FILE)\n";
00354 }
00355 if ((molecule->numMultipleDihedrals) && (simParameters->paraTypeCharmmOn))
00356 {
00357 iout << iINFO << molecule->numMultipleDihedrals
00358 << " DIHEDRALS WITH MULTIPLE PERIODICITY IGNORED (BASED ON PSF FILE) \n";
00359 iout << iINFO
00360 << " CHARMM MULTIPLICITIES BASED ON PARAMETER FILE INFO! \n";
00361 }
00362
00363
00364 if (molecule->numMultipleImpropers)
00365 {
00366 iout << iINFO << molecule->numMultipleImpropers
00367 << " IMPROPERS WITH MULTIPLE PERIODICITY\n";
00368 }
00369
00370 if (simParameters->constraintsOn)
00371 {
00372 iout << iINFO << molecule->numConstraints << " CONSTRAINTS\n";
00373 }
00374
00375 if (simParameters->consForceOn)
00376 iout << iINFO << molecule->numConsForce << " CONSTANT FORCES\n";
00377
00378 if (simParameters->stirOn)
00379 iout << iINFO << molecule->numStirredAtoms << " STIRRED ATOMS\n";
00380
00381 if (simParameters->fixedAtomsOn)
00382 {
00383 iout << iINFO << molecule->numFixedAtoms << " FIXED ATOMS\n";
00384 }
00385
00386 if (simParameters->rigidBonds)
00387 {
00388 iout << iINFO << molecule->numRigidBonds << " RIGID BONDS\n";
00389 }
00390
00391 if (simParameters->fixedAtomsOn && simParameters->rigidBonds)
00392 {
00393 iout << iINFO << molecule->numFixedRigidBonds <<
00394 " RIGID BONDS BETWEEN FIXED ATOMS\n";
00395 }
00396
00397
00398 if (simParameters->gridforceOn)
00399 {
00400 iout << iINFO << molecule->numGridforces << " GRID FORCES\n";
00401 }
00402
00403
00404
00405
00406 if (simParameters->fepOn || simParameters->thermInt) {
00407 iout << iINFO << molecule->numFepInitial <<
00408 " ATOMS TO DISAPPEAR IN FINAL STATE\n";
00409 iout << iINFO << molecule->numFepFinal <<
00410 " ATOMS TO APPEAR IN FINAL STATE\n";
00411 }
00412
00413
00414 if (simParameters->lesOn) {
00415 iout << iINFO << molecule->numFepInitial <<
00416 " LOCALLY ENHANCED ATOMS ENABLED\n";
00417 }
00418
00419 if (simParameters->pairInteractionOn) {
00420 iout << iINFO << "PAIR INTERACTION GROUP 1 CONTAINS "
00421 << molecule->numFepInitial << " ATOMS\n";
00422 if (!simParameters->pairInteractionSelf) {
00423 iout << iINFO << "PAIR INTERACTION GROUP 2 CONTAINS "
00424 << molecule->numFepFinal << " ATOMS\n";
00425 }
00426 }
00427
00428 #if 1
00429 if (molecule->numLonepairs != 0) {
00430 iout << iINFO << molecule->numLonepairs << " LONE PAIRS\n";
00431 }
00432 if (molecule->numDrudeAtoms != 0) {
00433 iout << iINFO << molecule->numDrudeAtoms << " DRUDE ATOMS\n";
00434 }
00435 iout << iINFO << molecule->num_deg_freedom(1)
00436 << " DEGREES OF FREEDOM\n";
00437 if (simParameters->drudeOn) {
00438 int g_bond = 3 * molecule->numDrudeAtoms;
00439 int g_com = molecule->num_deg_freedom(1) - g_bond;
00440 iout << iINFO << g_com << " DRUDE COM DEGREES OF FREEDOM\n";
00441 iout << iINFO << g_bond << " DRUDE BOND DEGREES OF FREEDOM\n";
00442 }
00443 #endif
00444 #if 0
00445 {
00446
00447 int numAtoms = molecule->numAtoms;
00448 int numDegFreedom = 3 * numAtoms;
00449 int numLonepairs = molecule->numLonepairs;
00450 int numFixedAtoms = molecule->numFixedAtoms;
00451 if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
00452 if ( ! ( numFixedAtoms || molecule->numConstraints
00453 || simParameters->comMove || simParameters->langevinOn ) ) {
00454 numDegFreedom -= 3;
00455 }
00456 if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
00457 int numRigidBonds = molecule->numRigidBonds;
00458 int numFixedRigidBonds = molecule->numFixedRigidBonds;
00459
00460
00461 numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
00462 iout << iINFO << numDegFreedom << " DEGREES OF FREEDOM\n";
00463 }
00464 #endif
00465
00466 iout << iINFO << molecule->numHydrogenGroups << " HYDROGEN GROUPS\n";
00467 if (simParameters->fixedAtomsOn)
00468 {
00469 iout << iINFO << molecule->numFixedGroups <<
00470 " HYDROGEN GROUPS WITH ALL ATOMS FIXED\n";
00471 }
00472
00473 {
00474 BigReal totalMass = 0;
00475 BigReal totalCharge = 0;
00476 int i;
00477 for ( i = 0; i < molecule->numAtoms; ++i ) {
00478 totalMass += molecule->atommass(i);
00479 totalCharge += molecule->atomcharge(i);
00480 }
00481 iout << iINFO << "TOTAL MASS = " << totalMass << " amu\n";
00482 iout << iINFO << "TOTAL CHARGE = " << totalCharge << " e\n";
00483 }
00484
00485 iout << iINFO << "*****************************\n";
00486 iout << endi;
00487 fflush(stdout);
00488
00489 DebugM(4, "::configFileInit() - printing Molecule Information\n");
00490
00491 molecule->print_atoms(parameters);
00492 molecule->print_bonds(parameters);
00493 molecule->print_exclusions();
00494 fflush(stdout);
00495
00496 DebugM(4, "::configFileInit() - done printing Molecule Information\n");
00497 DebugM(1, "::configFileInit() - done\n");
00498
00499 return(0);
00500 }
00501