00001 #include "common.h"
00002 #include "BackEnd.h"
00003 #include "InfoStream.h"
00004 #include "Node.h"
00005 #include "Molecule.h"
00006 #include "PDB.h"
00007 #include "PDBData.h"
00008 #include "ReductionMgr.h"
00009
00010 #include "colvarmodule.h"
00011 #include "colvaratoms.h"
00012 #include "colvarproxy.h"
00013 #include "colvarproxy_namd.h"
00014
00015
00016 colvarproxy_namd::colvarproxy_namd()
00017 {
00018 first_timestep = true;
00019 system_force_requested = false;
00020
00021
00022 simparams = Node::Object()->simParameters;
00023 lattice = &(simparams->lattice);
00024
00025 if (cvm::debug())
00026 iout << "Info: initializing the colvars proxy object.\n" << endi;
00027
00028
00029 StringList *config = Node::Object()->configList->find ("colvarsConfig");
00030 if (!config)
00031 NAMD_die ("No configuration file for collective variables: exiting.\n");
00032
00033
00034 StringList *input_restart = Node::Object()->configList->find ("colvarsInput");
00035 input_prefix_str = std::string (input_restart ? input_restart->data : "");
00036 if (input_prefix_str.rfind (".colvars.state") != std::string::npos) {
00037
00038 input_prefix_str.erase (input_prefix_str.rfind (".colvars.state"),
00039 std::string (".colvars.state").size());
00040 }
00041
00042
00043 if (simparams->rescaleFreq > 0)
00044 thermostat_temperature = simparams->rescaleTemp;
00045 else if (simparams->reassignFreq > 0)
00046 thermostat_temperature = simparams->reassignTemp;
00047 else if (simparams->langevinOn)
00048 thermostat_temperature = simparams->langevinTemp;
00049 else if (simparams->tCoupleOn)
00050 thermostat_temperature = simparams->tCoupleTemp;
00051
00052
00053 else
00054 thermostat_temperature = 0.0;
00055
00056 random = Random (simparams->randomSeed);
00057
00058
00059 output_prefix_str = std::string (simparams->outputFilename);
00060 restart_output_prefix_str = std::string (simparams->restartFilename);
00061 restart_frequency_s = simparams->restartFrequency;
00062
00063
00064
00065 colvars = new colvarmodule (config->data, this);
00066
00067 Node::Object()->colvars = colvars;
00068
00069 if (simparams->firstTimestep != 0) {
00070 cvm::log ("Initializing step number as firstTimestep.\n");
00071 colvars->it = colvars->it_restart = simparams->firstTimestep;
00072 }
00073
00074 if (cvm::debug()) {
00075 cvm::log ("colvars_atoms = "+cvm::to_str (colvars_atoms)+"\n");
00076 cvm::log ("colvars_atoms_ncopies = "+cvm::to_str (colvars_atoms_ncopies)+"\n");
00077 cvm::log ("positions = "+cvm::to_str (positions)+"\n");
00078 cvm::log ("total_forces = "+cvm::to_str (total_forces)+"\n");
00079 cvm::log ("applied_forces = "+cvm::to_str (applied_forces)+"\n");
00080 cvm::log (cvm::line_marker);
00081 }
00082
00083
00084 reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00085
00086 if (cvm::debug())
00087 iout << "Info: done initializing the colvars proxy object.\n" << endi;
00088 }
00089
00090
00091 colvarproxy_namd::~colvarproxy_namd()
00092 {
00093 delete reduction;
00094 if (colvars != NULL) {
00095 delete colvars;
00096 colvars = NULL;
00097 }
00098 }
00099
00100
00101 void colvarproxy_namd::calculate()
00102 {
00103
00104 if (first_timestep) {
00105 first_timestep = false;
00106 } else {
00107
00108 if ( step - previous_NAMD_step == 1 ) {
00109 colvars->it++;
00110 }
00111
00112
00113
00114
00115 }
00116 previous_NAMD_step = step;
00117
00118 if (cvm::debug()) {
00119 cvm::log (cvm::line_marker+
00120 "colvarproxy_namd, step no. "+cvm::to_str (colvars->it)+"\n"+
00121 "Updating internal data.\n");
00122 }
00123
00124
00125
00126 modifyForcedAtoms().resize (0);
00127 modifyAppliedForces().resize (0);
00128
00129
00130
00131 for (size_t i = 0; i < colvars_atoms.size(); i++) {
00132 positions[i] = cvm::rvector (0.0, 0.0, 0.0);
00133 total_forces[i] = cvm::rvector (0.0, 0.0, 0.0);
00134 applied_forces[i] = cvm::rvector (0.0, 0.0, 0.0);
00135 }
00136
00137
00138 for (size_t i = 0; i < colvars_atoms.size(); i++) {
00139 bool found_position = false;
00140 AtomIDList::const_iterator a_i = this->getAtomIdBegin();
00141 AtomIDList::const_iterator a_e = this->getAtomIdEnd();
00142 PositionList::const_iterator p_i = this->getAtomPositionBegin();
00143 for ( ; a_i != a_e; ++a_i, ++p_i ) {
00144 if ( *a_i == colvars_atoms[i] ) {
00145 found_position = true;
00146 Position const &namd_pos = *p_i;
00147 positions[i] = cvm::rvector (namd_pos.x, namd_pos.y, namd_pos.z);
00148 break;
00149 }
00150 }
00151 if (!found_position)
00152 cvm::fatal_error ("Error: cannot find the position of atom "+
00153 cvm::to_str (colvars_atoms[i]+1)+"\n");
00154 }
00155
00156
00157 if (system_force_requested && cvm::step_relative() > 0) {
00158
00159
00160
00161 for (size_t i = 0; i < colvars_atoms.size(); i++) {
00162 bool found_total_force = false;
00163
00164 AtomIDList::const_iterator a_i = this->getForceIdBegin();
00165 AtomIDList::const_iterator a_e = this->getForceIdEnd();
00166 PositionList::const_iterator f_i = this->getTotalForce();
00167 for ( ; a_i != a_e; ++a_i, ++f_i ) {
00168 if ( *a_i == colvars_atoms[i] ) {
00169 found_total_force = true;
00170 Vector const &namd_force = *f_i;
00171 total_forces[i] = cvm::rvector (namd_force.x, namd_force.y, namd_force.z);
00172
00173
00174
00175
00176 break;
00177 }
00178 }
00179 if (!found_total_force)
00180 cvm::fatal_error ("Error: system forces were requested, but total force on atom "+
00181 cvm::to_str (colvars_atoms[i]+1) + " was not\n"
00182 "found. The most probable cause is combination of energy minimization with a\n"
00183 "biasing method that requires MD (e.g. ABF). Always run minimization\n"
00184 "and ABF separately.");
00185 }
00186
00187
00188 for (size_t i = 0; i < colvars_atoms.size(); i++) {
00189 AtomIDList::const_iterator a_i = this->getLastAtomsForcedBegin();
00190 AtomIDList::const_iterator a_e = this->getLastAtomsForcedEnd();
00191 PositionList::const_iterator f_i = this->getLastForcesBegin();
00192 for ( ; a_i != a_e; ++a_i, ++f_i ) {
00193 if ( *a_i == colvars_atoms[i] ) {
00194 Vector const &namd_force = *f_i;
00195 if (cvm::debug())
00196 cvm::log ("Found a force applied to atom "+
00197 cvm::to_str (colvars_atoms[i]+1)+": "+
00198 cvm::to_str (cvm::rvector (namd_force.x, namd_force.y, namd_force.z))+
00199 "; current total is "+
00200 cvm::to_str (applied_forces[i])+".\n");
00201 applied_forces[i] += cvm::rvector (namd_force.x, namd_force.y, namd_force.z);
00202 }
00203 }
00204 }
00205 }
00206
00207
00208 colvars->calc();
00209
00210 reduction->submit();
00211
00212
00213
00214 if (step == simparams->N) {
00215 colvars->write_output_files();
00216 }
00217 }
00218
00219
00220 void colvarproxy_namd::add_energy (cvm::real energy)
00221 {
00222 reduction->item(REDUCTION_MISC_ENERGY) += energy;
00223 }
00224
00225 void colvarproxy_namd::request_system_force (bool yesno)
00226 {
00227 system_force_requested = yesno;
00228 }
00229
00230 void colvarproxy_namd::log (std::string const &message)
00231 {
00232 std::istringstream is (message);
00233 std::string line;
00234 while (std::getline (is, line))
00235 iout << "colvars: " << line << "\n";
00236 iout << endi;
00237 }
00238
00239
00240 void colvarproxy_namd::fatal_error (std::string const &message)
00241 {
00242 cvm::log (message);
00243 if (!cvm::debug())
00244 cvm::log ("If this error message is unclear, "
00245 "try recompiling with -DCOLVARS_DEBUG.\n");
00246 NAMD_die ("Error in the collective variables module: exiting.\n");
00247 }
00248
00249
00250 void colvarproxy_namd::exit (std::string const &message)
00251 {
00252 cvm::log (message);
00253 BackEnd::exit();
00254 }
00255
00256
00257 enum e_pdb_field {
00258 e_pdb_none,
00259 e_pdb_occ,
00260 e_pdb_beta,
00261 e_pdb_x,
00262 e_pdb_y,
00263 e_pdb_z,
00264 e_pdb_ntot
00265 };
00266
00267
00268 e_pdb_field pdb_field_str2enum (std::string const &pdb_field_str)
00269 {
00270 e_pdb_field pdb_field = e_pdb_none;
00271
00272 if (colvarparse::to_lower_cppstr (pdb_field_str) ==
00273 colvarparse::to_lower_cppstr ("O")) {
00274 pdb_field = e_pdb_occ;
00275 }
00276
00277 if (colvarparse::to_lower_cppstr (pdb_field_str) ==
00278 colvarparse::to_lower_cppstr ("B")) {
00279 pdb_field = e_pdb_beta;
00280 }
00281
00282 if (colvarparse::to_lower_cppstr (pdb_field_str) ==
00283 colvarparse::to_lower_cppstr ("X")) {
00284 pdb_field = e_pdb_x;
00285 }
00286
00287 if (colvarparse::to_lower_cppstr (pdb_field_str) ==
00288 colvarparse::to_lower_cppstr ("Y")) {
00289 pdb_field = e_pdb_y;
00290 }
00291
00292 if (colvarparse::to_lower_cppstr (pdb_field_str) ==
00293 colvarparse::to_lower_cppstr ("Z")) {
00294 pdb_field = e_pdb_z;
00295 }
00296
00297 if (pdb_field == e_pdb_none) {
00298 cvm::fatal_error ("Error: unsupported PDB field, \""+
00299 pdb_field_str+"\".\n");
00300 }
00301
00302 return pdb_field;
00303 }
00304
00305
00306 void colvarproxy_namd::load_coords (char const *pdb_filename,
00307 std::vector<cvm::atom_pos> &pos,
00308 const std::vector<int> &indices,
00309 std::string const pdb_field_str,
00310 double const pdb_field_value)
00311 {
00312 if (pdb_field_str.size() == 0 && indices.size() == 0) {
00313 cvm::fatal_error ("Bug alert: either PDB field should be defined or list of "
00314 "atom IDs should be available when loading atom coordinates!\n");
00315 }
00316
00317 e_pdb_field pdb_field_index;
00318 bool const use_pdb_field = (pdb_field_str.size() > 0);
00319 if (use_pdb_field) {
00320 pdb_field_index = pdb_field_str2enum (pdb_field_str);
00321 }
00322
00323
00324 std::vector<int>::const_iterator current_index = indices.begin();
00325
00326 PDB *pdb = new PDB (pdb_filename);
00327 size_t const pdb_natoms = pdb->num_atoms();
00328
00329 if (pos.size() != pdb_natoms) {
00330
00331 bool const pos_allocated = (pos.size() > 0);
00332
00333 size_t ipos = 0, ipdb = 0;
00334 for ( ; ipdb < pdb_natoms; ipdb++) {
00335
00336 if (use_pdb_field) {
00337
00338 double atom_pdb_field_value = 0.0;
00339
00340 switch (pdb_field_index) {
00341 case e_pdb_occ:
00342 atom_pdb_field_value = (pdb->atom (ipdb))->occupancy();
00343 break;
00344 case e_pdb_beta:
00345 atom_pdb_field_value = (pdb->atom (ipdb))->temperaturefactor();
00346 break;
00347 case e_pdb_x:
00348 atom_pdb_field_value = (pdb->atom (ipdb))->xcoor();
00349 break;
00350 case e_pdb_y:
00351 atom_pdb_field_value = (pdb->atom (ipdb))->ycoor();
00352 break;
00353 case e_pdb_z:
00354 atom_pdb_field_value = (pdb->atom (ipdb))->zcoor();
00355 break;
00356 default:
00357 break;
00358 }
00359
00360 if ( (pdb_field_value) &&
00361 (atom_pdb_field_value != pdb_field_value) ) {
00362 continue;
00363 } else if (atom_pdb_field_value == 0.0) {
00364 continue;
00365 }
00366
00367 } else {
00368
00369 if (ipdb != *current_index) {
00370
00371 continue;
00372 } else {
00373 current_index++;
00374 }
00375 }
00376
00377 if (!pos_allocated) {
00378 pos.push_back (cvm::atom_pos (0.0, 0.0, 0.0));
00379 } else if (ipos >= pos.size()) {
00380 cvm::fatal_error ("Error: the PDB file \""+
00381 std::string (pdb_filename)+
00382 "\" contains coordinates for "
00383 "more atoms than needed.\n");
00384 }
00385
00386 pos[ipos] = cvm::atom_pos ((pdb->atom (ipdb))->xcoor(),
00387 (pdb->atom (ipdb))->ycoor(),
00388 (pdb->atom (ipdb))->zcoor());
00389 ipos++;
00390 if (!use_pdb_field && current_index == indices.end())
00391 break;
00392 }
00393
00394 if ((ipos < pos.size()) || (current_index != indices.end()))
00395 cvm::fatal_error ("Error: the number of records in the PDB file \""+
00396 std::string (pdb_filename)+
00397 "\" does not appear to match either the total number of atoms,"+
00398 " or the number of coordinates requested at this point ("+
00399 cvm::to_str (pos.size())+").\n");
00400
00401 } else {
00402
00403
00404
00405 for (size_t ia = 0; ia < pos.size(); ia++) {
00406 pos[ia] = cvm::atom_pos ((pdb->atom (ia))->xcoor(),
00407 (pdb->atom (ia))->ycoor(),
00408 (pdb->atom (ia))->zcoor());
00409 }
00410 }
00411
00412 delete pdb;
00413 }
00414
00415
00416 void colvarproxy_namd::load_atoms (char const *pdb_filename,
00417 std::vector<cvm::atom> &atoms,
00418 std::string const pdb_field_str,
00419 double const pdb_field_value)
00420 {
00421 if (pdb_field_str.size() == 0)
00422 cvm::fatal_error ("Error: must define which PDB field to use "
00423 "in order to define atoms from a PDB file.\n");
00424
00425 PDB *pdb = new PDB (pdb_filename);
00426 size_t const pdb_natoms = pdb->num_atoms();
00427
00428 e_pdb_field pdb_field_index = pdb_field_str2enum (pdb_field_str);
00429
00430 for (size_t ipdb = 0; ipdb < pdb_natoms; ipdb++) {
00431
00432 double atom_pdb_field_value = 0.0;
00433
00434 switch (pdb_field_index) {
00435 case e_pdb_occ:
00436 atom_pdb_field_value = (pdb->atom (ipdb))->occupancy();
00437 break;
00438 case e_pdb_beta:
00439 atom_pdb_field_value = (pdb->atom (ipdb))->temperaturefactor();
00440 break;
00441 case e_pdb_x:
00442 atom_pdb_field_value = (pdb->atom (ipdb))->xcoor();
00443 break;
00444 case e_pdb_y:
00445 atom_pdb_field_value = (pdb->atom (ipdb))->ycoor();
00446 break;
00447 case e_pdb_z:
00448 atom_pdb_field_value = (pdb->atom (ipdb))->zcoor();
00449 break;
00450 default:
00451 break;
00452 }
00453
00454 if ( (pdb_field_value) &&
00455 (atom_pdb_field_value != pdb_field_value) ) {
00456 continue;
00457 } else if (atom_pdb_field_value == 0.0) {
00458 continue;
00459 }
00460
00461 atoms.push_back (cvm::atom (ipdb+1));
00462 }
00463
00464 delete pdb;
00465 }
00466
00467
00468 void colvarproxy_namd::backup_file (char const *filename)
00469 {
00470 if (std::string (filename).rfind (std::string (".colvars.state")) != std::string::npos) {
00471 NAMD_backup_file (filename, ".old");
00472 } else {
00473 NAMD_backup_file (filename, ".BAK");
00474 }
00475 }
00476
00477
00478 size_t colvarproxy_namd::init_namd_atom (AtomID const &aid)
00479 {
00480 modifyRequestedAtoms().add (aid);
00481 for (size_t i = 0; i < colvars_atoms.size(); i++) {
00482 if (colvars_atoms[i] == aid) {
00483
00484 colvars_atoms_ncopies[i] += 1;
00485 return i;
00486 }
00487 }
00488
00489
00490 colvars_atoms_ncopies.push_back (1);
00491 colvars_atoms.push_back (aid);
00492 positions.push_back (cvm::rvector());
00493 total_forces.push_back (cvm::rvector());
00494 applied_forces.push_back (cvm::rvector());
00495
00496 return (colvars_atoms.size()-1);
00497 }
00498
00499
00500
00501 cvm::atom::atom (int const &atom_number)
00502 {
00503
00504 AtomID const aid (atom_number-1);
00505
00506 if (cvm::debug())
00507 cvm::log ("Adding atom "+cvm::to_str (aid+1)+
00508 " for collective variables calculation.\n");
00509
00510 if ( (aid < 0) || (aid >= Node::Object()->molecule->numAtoms) )
00511 cvm::fatal_error ("Error: invalid atom number specified, "+
00512 cvm::to_str (atom_number)+"\n");
00513 this->index = ((colvarproxy_namd *) cvm::proxy)->init_namd_atom (aid);
00514 if (cvm::debug())
00515 cvm::log ("The index of this atom in the colvarproxy_namd arrays is "+
00516 cvm::to_str (this->index)+".\n");
00517 this->id = aid;
00518 this->mass = Node::Object()->molecule->atommass (aid);
00519 this->reset_data();
00520 }
00521
00522
00526 cvm::atom::atom (cvm::residue_id const &residue,
00527 std::string const &atom_name,
00528 std::string const &segment_id)
00529 {
00530 AtomID const aid =
00531 (segment_id.size() ?
00532 Node::Object()->molecule->get_atom_from_name (segment_id.c_str(),
00533 residue,
00534 atom_name.c_str()) :
00535 Node::Object()->molecule->get_atom_from_name ("MAIN",
00536 residue,
00537 atom_name.c_str()));
00538
00539
00540 if (cvm::debug())
00541 cvm::log ("Adding atom \""+
00542 atom_name+"\" in residue "+
00543 cvm::to_str (residue)+
00544 " (index "+cvm::to_str (aid)+
00545 ") for collective variables calculation.\n");
00546
00547 if (aid < 0) {
00548
00549 cvm::fatal_error ("Error: could not find atom \""+
00550 atom_name+"\" in residue "+
00551 cvm::to_str (residue)+
00552 ( (segment_id != "MAIN") ?
00553 (", segment \""+segment_id+"\"") :
00554 ("") )+
00555 "\n");
00556 }
00557
00558 this->index = ((colvarproxy_namd *) cvm::proxy)->init_namd_atom (aid);
00559 if (cvm::debug())
00560 cvm::log ("The index of this atom in the colvarproxy_namd arrays is "+
00561 cvm::to_str (this->index)+".\n");
00562 this->id = aid;
00563 this->mass = Node::Object()->molecule->atommass (aid);
00564 this->reset_data();
00565 }
00566
00567
00568
00569 cvm::atom::atom (cvm::atom const &a)
00570 : index (a.index), id (a.id), mass (a.mass)
00571 {
00572
00573
00574
00575
00576 colvarproxy_namd *gm = (colvarproxy_namd *) cvm::proxy;
00577 gm->colvars_atoms_ncopies[this->index] += 1;
00578 }
00579
00580
00581 cvm::atom::~atom()
00582 {
00583 colvarproxy_namd *gm = (colvarproxy_namd *) cvm::proxy;
00584 if (gm->colvars_atoms_ncopies[this->index] > 0)
00585 gm->colvars_atoms_ncopies[this->index] -= 1;
00586 }
00587
00588
00589 void cvm::atom::read_position()
00590 {
00591 colvarproxy_namd const * const gm = (colvarproxy_namd *) cvm::proxy;
00592 this->pos = gm->positions[this->index];
00593 }
00594
00595
00596 void cvm::atom::read_velocity()
00597 {
00598 cvm::fatal_error ("Error: NAMD does not have yet a way to communicate "
00599 "atom velocities to the colvars.\n");
00600 }
00601
00602
00603 void cvm::atom::read_system_force()
00604 {
00605 colvarproxy_namd const * const gm = (colvarproxy_namd *) cvm::proxy;
00606 this->system_force = gm->total_forces[this->index] - gm->applied_forces[this->index];
00607 }
00608
00609
00610 void cvm::atom::apply_force (cvm::rvector const &new_force)
00611 {
00612 colvarproxy_namd *gm = (colvarproxy_namd *) cvm::proxy;
00613 gm->modifyForcedAtoms().add (this->id);
00614 gm->modifyAppliedForces().add (Vector (new_force.x, new_force.y, new_force.z));
00615 }
00616