Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

colvarproxy_namd.C

Go to the documentation of this file.
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 
00009 #include "colvarmodule.h"
00010 #include "colvaratoms.h"
00011 #include "colvarproxy.h"
00012 #include "colvarproxy_namd.h"
00013 
00014 
00015 colvarproxy_namd::colvarproxy_namd()
00016 {
00017   // initialize pointers to NAMD configuration data
00018   simparams = Node::Object()->simParameters;
00019   lattice = &(simparams->lattice);
00020 
00021   if (cvm::debug())
00022     iout << "Info: initializing the colvars proxy object.\n" << endi;
00023 
00024   // find the configuration file
00025   StringList *config = Node::Object()->configList->find ("colvarsConfig");
00026   if (!config)
00027     NAMD_die ("No configuration file for collective variables: exiting.\n");
00028 
00029   // find the input state file
00030   StringList *input_restart = Node::Object()->configList->find ("colvarsInput");
00031   input_prefix_str = std::string (input_restart ? input_restart->data : "");
00032   if (input_prefix_str.rfind (".colvars.state") != std::string::npos) {
00033     // strip the extension, if present
00034     input_prefix_str.erase (input_prefix_str.rfind (".colvars.state"),
00035                             std::string (".colvars.state").size());
00036   }
00037 
00038   // get the thermostat temperature
00039   if (simparams->rescaleFreq > 0)
00040     thermostat_temperature = simparams->rescaleTemp;
00041   else if (simparams->reassignFreq > 0)
00042     thermostat_temperature = simparams->reassignTemp;
00043   else if (simparams->langevinOn)
00044     thermostat_temperature = simparams->langevinTemp;
00045   else if (simparams->tCoupleOn)
00046     thermostat_temperature = simparams->tCoupleTemp;
00047   else 
00048     thermostat_temperature = 0.0;
00049 
00050   // take the output prefixes from the namd input
00051   output_prefix_str = std::string (simparams->outputFilename);
00052   restart_output_prefix_str = std::string (simparams->restartFilename);
00053   restart_frequency_s = simparams->restartFrequency;
00054 
00055   // initiate the colvarmodule, this object will be the communication
00056   // proxy
00057   colvars = new colvarmodule (config->data, this);
00058   colvars->dt = simparams->dt;
00059 
00060   if (cvm::debug()) {
00061     cvm::log ("colvars_atoms = "+cvm::to_str (colvars_atoms)+"\n");
00062     cvm::log ("colvars_atoms_ncopies = "+cvm::to_str (colvars_atoms_ncopies)+"\n");
00063     cvm::log ("positions = "+cvm::to_str (positions)+"\n");
00064     cvm::log ("total_forces = "+cvm::to_str (total_forces)+"\n");
00065     cvm::log ("applied_forces = "+cvm::to_str (applied_forces)+"\n");
00066     cvm::log (cvm::line_marker);
00067   }
00068 
00069   if (cvm::debug())
00070     iout << "Info: done initializing the colvars proxy object.\n" << endi;
00071 }
00072 
00073 
00074 colvarproxy_namd::~colvarproxy_namd()
00075 {
00076   if (colvars != NULL) {
00077     delete colvars;
00078     colvars = NULL;
00079   }
00080 }
00081 
00082 // Reimplemented function from GlobalMaster
00083 void colvarproxy_namd::calculate()
00084 {
00085   if (cvm::debug()) {
00086     cvm::log (cvm::line_marker+
00087               "colvarproxy_namd, step no. "+cvm::to_str (colvars->it)+"\n"+
00088               "Updating internal data.\n");
00089   }
00090 
00091   // must delete the forces applied at the previous step: they have
00092   // already been used and copied to other memory locations
00093   modifyForcedAtoms().resize (0);
00094   modifyAppliedForces().resize (0);
00095 
00096   // prepare the local arrays to contain the sorted copies of the NAMD
00097   // arrays
00098   for (size_t i = 0; i < colvars_atoms.size(); i++) {
00099     positions[i] = cvm::rvector (0.0, 0.0, 0.0);
00100     total_forces[i] = cvm::rvector (0.0, 0.0, 0.0);
00101     applied_forces[i] = cvm::rvector (0.0, 0.0, 0.0);
00102   }
00103 
00104   // sort the positions array
00105   for (size_t i = 0; i < colvars_atoms.size(); i++) {
00106     bool found_position = false;
00107     AtomIDList::const_iterator a_i = this->getAtomIdBegin();
00108     AtomIDList::const_iterator a_e = this->getAtomIdEnd();
00109     PositionList::const_iterator p_i = this->getAtomPositionBegin();
00110     for ( ; a_i != a_e; ++a_i, ++p_i ) {
00111       if ( *a_i == colvars_atoms[i] ) {
00112         found_position = true;
00113         Position const &namd_pos = *p_i;
00114         positions[i] = cvm::rvector (namd_pos.x, namd_pos.y, namd_pos.z);
00115         break;
00116       }
00117     }
00118     if (!found_position)
00119       cvm::fatal_error ("Error: cannot find the position of atom "+
00120                         cvm::to_str (colvars_atoms[i]+1)+"\n");
00121   }
00122 
00123 
00124   if (cvm::step_relative() > 0) {
00125 
00126     // sort the array of total forces from the previous step (but only
00127     // do it if there *is* a previous step!)
00128     for (size_t i = 0; i < colvars_atoms.size(); i++) {
00129       bool found_total_force = false;
00130       AtomIDList::const_iterator a_i = this->getForceIdBegin();
00131       AtomIDList::const_iterator a_e = this->getForceIdEnd();
00132       PositionList::const_iterator f_i = this->getTotalForce();
00133       for ( ; a_i != a_e; ++a_i, ++f_i ) {
00134         if ( *a_i == colvars_atoms[i] ) {
00135           found_total_force = true;
00136           Vector const &namd_force = *f_i;
00137           total_forces[i] = cvm::rvector (namd_force.x, namd_force.y, namd_force.z);
00138           //           if (cvm::debug()) 
00139           //             cvm::log ("Found the total force of atom "+
00140           //                       cvm::to_str (colvars_atoms[i]+1)+", which is "+
00141           //                       cvm::to_str (total_forces[i])+".\n");
00142           break;
00143         }
00144       }
00145       if (!found_total_force)
00146         cvm::fatal_error ("Error: cannot find the total force of atom "+
00147                           cvm::to_str (colvars_atoms[i]+1)+". Is energy minimization"+
00148                           " (instead of MD) enabled?\n");
00149     }
00150 
00151     // do the same for applied forces
00152     for (size_t i = 0; i < colvars_atoms.size(); i++) {
00153       AtomIDList::const_iterator a_i = this->getLastAtomsForcedBegin();
00154       AtomIDList::const_iterator a_e = this->getLastAtomsForcedEnd();
00155       PositionList::const_iterator f_i = this->getLastForcesBegin();
00156       for ( ; a_i != a_e; ++a_i, ++f_i ) {
00157         if ( *a_i == colvars_atoms[i] ) {
00158           Vector const &namd_force = *f_i;
00159           if (cvm::debug())
00160             cvm::log ("Found a force applied to atom "+
00161                       cvm::to_str (colvars_atoms[i]+1)+": "+
00162                       cvm::to_str (cvm::rvector (namd_force.x, namd_force.y, namd_force.z))+
00163                       "; current total is "+
00164                       cvm::to_str (applied_forces[i])+".\n");
00165           applied_forces[i] += cvm::rvector (namd_force.x, namd_force.y, namd_force.z);
00166         }
00167       }
00168     }
00169   }
00170 
00171   // call the collective variable module
00172   colvars->calc();
00173 
00174   // NAMD does not properly destruct GlobalMaster objects at the end
00175   // of the job
00176   if (colvars->it >= (colvars->it_restart + simparams->N)) {
00177     if (cvm::debug())
00178       iout << "Info: Deallocating the colvar module object.\n" << endi;
00179     colvars->finalise();
00180   }
00181 
00182   // increment the step number: XX TODO: syncronise with NAMD's counter?
00183   (colvars->it)++;
00184 }
00185 
00186 
00187 void colvarproxy_namd::log (std::string const &message)
00188 {
00189   std::istringstream is (message);
00190   std::string line;
00191   while (std::getline (is, line))
00192     iout << "colvars: " << line << "\n";
00193   iout << endi;
00194 }
00195 
00196 
00197 void colvarproxy_namd::fatal_error (std::string const &message)
00198 {
00199   cvm::log (message);
00200   if (!cvm::debug())
00201     cvm::log ("If this error message is unclear, "
00202               "try recompile with -DCOLVARS_DEBUG.\n");
00203   NAMD_die ("Error in the collective variables module: exiting.\n");
00204 }
00205 
00206 
00207 void colvarproxy_namd::exit (std::string const &message)
00208 {
00209   cvm::log (message);
00210   BackEnd::exit();
00211 }
00212 
00213 
00214 enum e_pdb_field {
00215   e_pdb_none,
00216   e_pdb_occ,
00217   e_pdb_beta,
00218   e_pdb_x,
00219   e_pdb_y,
00220   e_pdb_z,
00221   e_pdb_ntot
00222 };
00223 
00224 
00225 e_pdb_field pdb_field_str2enum (std::string const &pdb_field_str)
00226 {
00227   e_pdb_field pdb_field = e_pdb_none;
00228 
00229   if (colvarparse::to_lower_cppstr (pdb_field_str) ==
00230       colvarparse::to_lower_cppstr ("O")) {
00231     pdb_field = e_pdb_occ;
00232   }
00233 
00234   if (colvarparse::to_lower_cppstr (pdb_field_str) ==
00235       colvarparse::to_lower_cppstr ("B")) {
00236     pdb_field = e_pdb_beta;
00237   }
00238 
00239   if (colvarparse::to_lower_cppstr (pdb_field_str) ==
00240       colvarparse::to_lower_cppstr ("X")) {
00241     pdb_field = e_pdb_x;
00242   }
00243   
00244   if (colvarparse::to_lower_cppstr (pdb_field_str) ==
00245       colvarparse::to_lower_cppstr ("Y")) {
00246     pdb_field = e_pdb_y;
00247   }
00248 
00249   if (colvarparse::to_lower_cppstr (pdb_field_str) ==
00250       colvarparse::to_lower_cppstr ("Z")) {
00251     pdb_field = e_pdb_z;
00252   }
00253 
00254   if (pdb_field == e_pdb_none) {
00255     cvm::fatal_error ("Error: unsupported PDB field, \""+
00256                       pdb_field_str+"\".\n");
00257   }
00258 
00259   return pdb_field;
00260 }
00261 
00262 
00263 void colvarproxy_namd::load_coords (char const *pdb_filename,
00264                                     std::vector<cvm::atom_pos> &pos,
00265                                     std::string const pdb_field_str,
00266                                     double const pdb_field_value)
00267 {
00268   if (pdb_field_str.size() == 0)
00269     cvm::fatal_error ("Error: must define which PDB field to use "
00270                       "in order to define atoms from a PDB file.\n");
00271 
00272   e_pdb_field pdb_field_index = pdb_field_str2enum (pdb_field_str);
00273 
00274   PDB *pdb = new PDB (pdb_filename);
00275   size_t const pdb_natoms = pdb->num_atoms();
00276   
00277   if (pos.size() != pdb_natoms) {
00278 
00279     bool const pos_allocated = (pos.size() > 0);
00280 
00281     size_t ipos = 0, ipdb = 0;
00282     for ( ; ipdb < pdb_natoms; ipdb++) {
00283 
00284       double atom_pdb_field_value = 0.0;
00285 
00286       switch (pdb_field_index) {
00287       case e_pdb_occ:
00288         atom_pdb_field_value = (pdb->atom (ipdb))->occupancy();
00289         break;
00290       case e_pdb_beta:
00291         atom_pdb_field_value = (pdb->atom (ipdb))->temperaturefactor();
00292         break;
00293       case e_pdb_x:
00294         atom_pdb_field_value = (pdb->atom (ipdb))->xcoor();
00295         break;
00296       case e_pdb_y:
00297         atom_pdb_field_value = (pdb->atom (ipdb))->ycoor();
00298         break;
00299       case e_pdb_z:
00300         atom_pdb_field_value = (pdb->atom (ipdb))->zcoor();
00301         break;
00302       default:
00303         break;
00304       }
00305 
00306       if ( (pdb_field_value) &&
00307            (atom_pdb_field_value != pdb_field_value) ) {
00308         continue;
00309       } else if (atom_pdb_field_value == 0.0) {
00310         continue;
00311       }
00312       
00313       if (!pos_allocated)
00314         pos.push_back (cvm::atom_pos (0.0, 0.0, 0.0));
00315 
00316       if (ipos < pos.size()) {
00317         pos[ipos] = cvm::atom_pos ((pdb->atom (ipdb))->xcoor(),
00318                                    (pdb->atom (ipdb))->ycoor(),
00319                                    (pdb->atom (ipdb))->zcoor());
00320         ipos++;
00321       } else {
00322         cvm::fatal_error ("Error: the PDB file \""+
00323                           std::string (pdb_filename)+
00324                           "\" contains coordinates for "
00325                           "more atoms than needed ("+
00326                           cvm::to_str (pos.size())+").\n");
00327       }
00328     }
00329 
00330     if (ipos < pos.size())
00331       cvm::fatal_error ("Error: the PDB file \""+
00332                         std::string (pdb_filename)+
00333                         "\" contains coordinates for only "+
00334                         cvm::to_str (ipdb)+
00335                         " atoms, but "+cvm::to_str (pos.size())+
00336                         " are needed.\n");
00337 
00338   } else {
00339 
00340     // when the PDB contains exactly the number of atoms of the array,
00341     // ignore the fields and just read coordinates
00342     for (size_t ia = 0; ia < pos.size(); ia++) {
00343       pos[ia] = cvm::atom_pos ((pdb->atom (ia))->xcoor(),
00344                                (pdb->atom (ia))->ycoor(),
00345                                (pdb->atom (ia))->zcoor());
00346     }
00347   }
00348 
00349   delete pdb;
00350 }
00351 
00352 
00353 void colvarproxy_namd::load_atoms (char const *pdb_filename,
00354                                    std::vector<cvm::atom> &atoms,
00355                                    std::string const pdb_field_str,
00356                                    double const pdb_field_value)
00357 {
00358   if (pdb_field_str.size() == 0)
00359     cvm::fatal_error ("Error: must define which PDB field to use "
00360                       "in order to define atoms from a PDB file.\n");
00361 
00362   PDB *pdb = new PDB (pdb_filename);
00363   size_t const pdb_natoms = pdb->num_atoms();
00364 
00365   e_pdb_field pdb_field_index = pdb_field_str2enum (pdb_field_str);
00366 
00367   for (size_t ipdb = 0; ipdb < pdb_natoms; ipdb++) {
00368 
00369     double atom_pdb_field_value = 0.0;
00370 
00371     switch (pdb_field_index) {
00372     case e_pdb_occ:
00373       atom_pdb_field_value = (pdb->atom (ipdb))->occupancy();
00374       break;
00375     case e_pdb_beta:
00376       atom_pdb_field_value = (pdb->atom (ipdb))->temperaturefactor();
00377       break;
00378     case e_pdb_x:
00379       atom_pdb_field_value = (pdb->atom (ipdb))->xcoor();
00380       break;
00381     case e_pdb_y:
00382       atom_pdb_field_value = (pdb->atom (ipdb))->ycoor();
00383       break;
00384     case e_pdb_z:
00385       atom_pdb_field_value = (pdb->atom (ipdb))->zcoor();
00386       break;
00387     default:
00388       break;
00389     }
00390 
00391     if ( (pdb_field_value) &&
00392          (atom_pdb_field_value != pdb_field_value) ) {
00393       continue;
00394     } else if (atom_pdb_field_value == 0.0) {
00395       continue;
00396     }
00397      
00398     atoms.push_back (cvm::atom (ipdb+1));
00399   }
00400 
00401   delete pdb;
00402 }
00403 
00404 
00405 void colvarproxy_namd::backup_file (char const *filename)
00406 {
00407   if (std::string (filename).rfind (std::string (".colvars.state")) != std::string::npos) {
00408     NAMD_backup_file (filename, ".old");
00409   } else {
00410     NAMD_backup_file (filename, ".BAK");
00411   }
00412 }
00413 
00414 
00415 size_t colvarproxy_namd::init_namd_atom (AtomID const &aid)
00416 {
00417   modifyRequestedAtoms().add (aid);
00418   for (size_t i = 0; i < colvars_atoms.size(); i++) {
00419     if (colvars_atoms[i] == aid) {
00420       // this atom id was already recorded
00421       colvars_atoms_ncopies[i] += 1;
00422       return i;
00423     }
00424   }
00425 
00426   // allocate a new slot for this atom
00427   colvars_atoms_ncopies.push_back (1);
00428   colvars_atoms.push_back (aid);
00429   positions.push_back (cvm::rvector());
00430   total_forces.push_back (cvm::rvector());
00431   applied_forces.push_back (cvm::rvector());
00432 
00433   return (colvars_atoms.size()-1);
00434 }
00435 
00436 
00437 // atom member functions, NAMD specific implementations
00438 
00439 cvm::atom::atom (int const &atom_number)
00440 {
00441   // NAMD internal numbering starts from zero
00442   AtomID const aid (atom_number-1);
00443 
00444   if (cvm::debug())
00445     cvm::log ("Adding atom "+cvm::to_str (aid+1)+
00446               " for collective variables calculation.\n");
00447 
00448   if ( (aid < 0) || (aid >= Node::Object()->molecule->numAtoms) ) 
00449     cvm::fatal_error ("Error: invalid atom number specified, "+
00450                       cvm::to_str (atom_number)+"\n");
00451   this->index = ((colvarproxy_namd *) cvm::proxy)->init_namd_atom (aid);
00452   if (cvm::debug())
00453     cvm::log ("The index of this atom in the colvarproxy_namd arrays is "+
00454               cvm::to_str (this->index)+".\n");
00455   this->id = aid;
00456   this->mass = Node::Object()->molecule->atommass (aid);
00457   this->reset_data();
00458 }
00459 
00460 
00464 cvm::atom::atom (cvm::residue_id const &residue,
00465                  std::string const     &atom_name,
00466                  std::string const     &segment_id = "MAIN")
00467 {
00468   AtomID const aid =
00469     Node::Object()->molecule->get_atom_from_name (segment_id.c_str(),
00470                                                   residue,
00471                                                   atom_name.c_str());
00472 
00473   if (cvm::debug())
00474     cvm::log ("Adding atom \""+
00475               atom_name+"\" in residue "+
00476               cvm::to_str (residue)+
00477               " (index "+cvm::to_str (aid)+
00478               ") for collective variables calculation.\n");
00479 
00480   if (aid < 0) {
00481     // get_atom_from_name() has returned an error value
00482     cvm::fatal_error ("Error: could not find atom \""+
00483                       atom_name+"\" in residue "+
00484                       cvm::to_str (residue)+
00485                       ( (segment_id != "MAIN") ?
00486                         (", segment \""+segment_id+"\"") :
00487                         ("") )+
00488                       "\n");
00489   }
00490 
00491   this->index = ((colvarproxy_namd *) cvm::proxy)->init_namd_atom (aid);
00492   if (cvm::debug())
00493     cvm::log ("The index of this atom in the colvarproxy_namd arrays is "+
00494               cvm::to_str (this->index)+".\n");
00495   this->id = aid;
00496   this->mass = Node::Object()->molecule->atommass (aid);
00497   this->reset_data();
00498 }
00499 
00500 
00501 // copy constructor
00502 cvm::atom::atom (cvm::atom const &a)
00503   : index (a.index), id (a.id), mass (a.mass)
00504 {
00505   // init_namd_atom() has already been called by a's constructor, no
00506   // need to call it again
00507 
00508   // need to increment the counter anyway
00509   colvarproxy_namd *gm = (colvarproxy_namd *) cvm::proxy;
00510   gm->colvars_atoms_ncopies[this->index] += 1;
00511 }
00512 
00513 
00514 cvm::atom::~atom() 
00515 {
00516   colvarproxy_namd *gm = (colvarproxy_namd *) cvm::proxy;
00517   if (gm->colvars_atoms_ncopies[this->index] > 0)
00518     gm->colvars_atoms_ncopies[this->index] -= 1;
00519 }
00520 
00521 
00522 void cvm::atom::read_position()
00523 {
00524   colvarproxy_namd const * const gm = (colvarproxy_namd *) cvm::proxy;
00525   this->pos = gm->positions[this->index];
00526 }
00527 
00528 
00529 void cvm::atom::read_velocity()
00530 {
00531   cvm::fatal_error ("Error: NAMD does not have yet a way to communicate "
00532                     "atom velocities to the colvars.\n");
00533 }
00534 
00535 
00536 void cvm::atom::read_system_force()
00537 {
00538   colvarproxy_namd const * const gm = (colvarproxy_namd *) cvm::proxy;
00539   this->system_force = gm->total_forces[this->index] - gm->applied_forces[this->index];
00540 }
00541 
00542 
00543 void cvm::atom::apply_force (cvm::rvector const &new_force)
00544 {
00545   colvarproxy_namd *gm = (colvarproxy_namd *) cvm::proxy;
00546   gm->modifyForcedAtoms().add (this->id);
00547   gm->modifyAppliedForces().add (Vector (new_force.x, new_force.y, new_force.z));
00548 }
00549 

Generated on Mon Nov 23 04:59:19 2009 for NAMD by  doxygen 1.3.9.1