colvarproxy_namd.C

Go to the documentation of this file.
00001 // -*- c++ -*-
00002 
00003 // This file is part of the Collective Variables module (Colvars).
00004 // The original version of Colvars and its updates are located at:
00005 // https://github.com/colvars/colvars
00006 // Please update all Colvars source files before making any changes.
00007 // If you wish to distribute your changes, please submit them to the
00008 // Colvars repository at GitHub.
00009 
00010 #include <errno.h>
00011 
00012 #include "common.h"
00013 #include "fstream_namd.h"
00014 #include "BackEnd.h"
00015 #include "InfoStream.h"
00016 #include "Node.h"
00017 #include "Molecule.h"
00018 #include "PDB.h"
00019 #include "PDBData.h"
00020 #include "ReductionMgr.h"
00021 #include "ScriptTcl.h"
00022 
00023 #ifdef NAMD_TCL
00024 #include <tcl.h>
00025 #endif
00026 
00027 #include "colvarmodule.h"
00028 #include "colvaratoms.h"
00029 #include "colvarproxy.h"
00030 #include "colvarproxy_namd.h"
00031 #include "colvarscript.h"
00032 
00033 
00034 colvarproxy_namd::colvarproxy_namd()
00035 {
00036   first_timestep = true;
00037   total_force_requested = false;
00038   requestTotalForce(total_force_requested);
00039 
00040   // initialize pointers to NAMD configuration data
00041   simparams = Node::Object()->simParameters;
00042 
00043   if (cvm::debug())
00044     iout << "Info: initializing the colvars proxy object.\n" << endi;
00045 
00046   // find the configuration file, if provided
00047   StringList *config = Node::Object()->configList->find("colvarsConfig");
00048 
00049   // find the input state file
00050   StringList *input_restart = Node::Object()->configList->find("colvarsInput");
00051   input_prefix_str = std::string(input_restart ? input_restart->data : "");
00052   if (input_prefix_str.rfind(".colvars.state") != std::string::npos) {
00053     // strip the extension, if present
00054     input_prefix_str.erase(input_prefix_str.rfind(".colvars.state"),
00055                            std::string(".colvars.state").size());
00056   }
00057 
00058   // get the thermostat temperature
00059   if (simparams->rescaleFreq > 0)
00060     thermostat_temperature = simparams->rescaleTemp;
00061   else if (simparams->reassignFreq > 0)
00062     thermostat_temperature = simparams->reassignTemp;
00063   else if (simparams->langevinOn)
00064     thermostat_temperature = simparams->langevinTemp;
00065   else if (simparams->tCoupleOn)
00066     thermostat_temperature = simparams->tCoupleTemp;
00067   //else if (simparams->loweAndersenOn)
00068   //  thermostat_temperature = simparams->loweAndersenTemp;
00069   else
00070     thermostat_temperature = 0.0;
00071 
00072   random = Random(simparams->randomSeed);
00073 
00074   // take the output prefixes from the namd input
00075   output_prefix_str = std::string(simparams->outputFilename);
00076   restart_output_prefix_str = std::string(simparams->restartFilename);
00077   restart_frequency_s = simparams->restartFrequency;
00078 
00079   // check if it is possible to save output configuration
00080   if ((!output_prefix_str.size()) && (!restart_output_prefix_str.size())) {
00081     fatal_error("Error: neither the final output state file or "
00082                 "the output restart file could be defined, exiting.\n");
00083   }
00084 
00085 
00086 #ifdef NAMD_TCL
00087   have_scripts = true;
00088   // Store pointer to NAMD's Tcl interpreter
00089   interp = Node::Object()->getScript()->interp;
00090 
00091   // See is user-scripted forces are defined
00092   if (Tcl_FindCommand(interp, "calc_colvar_forces", NULL, 0) == NULL) {
00093     force_script_defined = false;
00094   } else {
00095     force_script_defined = true;
00096   }
00097 #else
00098   force_script_defined = false;
00099   have_scripts = false;
00100 #endif
00101 
00102 
00103   // initiate module: this object will be the communication proxy
00104   colvars = new colvarmodule(this);
00105   log("Using NAMD interface, version "+
00106       cvm::to_str(COLVARPROXY_VERSION)+".\n");
00107 
00108   if (config) {
00109     colvars->read_config_file(config->data);
00110   }
00111 
00112   setup();
00113   colvars->setup();
00114   colvars->setup_input();
00115   colvars->setup_output();
00116 
00117   // save to Node for Tcl script access
00118   Node::Object()->colvars = colvars;
00119 
00120 #ifdef NAMD_TCL
00121   // Construct instance of colvars scripting interface
00122   script = new colvarscript(this);
00123 #endif
00124 
00125   if (simparams->firstTimestep != 0) {
00126     log("Initializing step number as firstTimestep.\n");
00127     colvars->it = colvars->it_restart = simparams->firstTimestep;
00128   }
00129 
00130   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00131 
00132   if (cvm::debug())
00133     iout << "Info: done initializing the colvars proxy object.\n" << endi;
00134 }
00135 
00136 
00137 colvarproxy_namd::~colvarproxy_namd()
00138 {
00139   delete reduction;
00140   if (script != NULL) {
00141     delete script;
00142     script = NULL;
00143   }
00144   if (colvars != NULL) {
00145     delete colvars;
00146     colvars = NULL;
00147   }
00148 }
00149 
00150 
00151 int colvarproxy_namd::update_atoms_map(AtomIDList::const_iterator begin,
00152                                        AtomIDList::const_iterator end)
00153 {
00154   for (AtomIDList::const_iterator a_i = begin; a_i != end; a_i++) {
00155 
00156     if (cvm::debug()) {
00157       cvm::log("Updating atoms_map for atom ID "+cvm::to_str(*a_i)+"\n");
00158     }
00159 
00160     if (atoms_map[*a_i] >= 0) continue;
00161 
00162     for (size_t i = 0; i < atoms_ids.size(); i++) {
00163       if (atoms_ids[i] == *a_i) {
00164         atoms_map[*a_i] = i;
00165         break;
00166       }
00167     }
00168 
00169     if (atoms_map[*a_i] < 0) {
00170       // this atom is probably managed by another GlobalMaster:
00171       // add it here anyway to avoid having to test for array boundaries at each step
00172       int const index = add_atom_slot(*a_i);
00173       atoms_map[*a_i] = index;
00174       update_atom_properties(index);
00175     }
00176   }
00177 
00178   if (cvm::debug()) {
00179     log("atoms_map = "+cvm::to_str(atoms_map)+".\n");
00180   }
00181 
00182   return COLVARS_OK;
00183 }
00184 
00185 
00186 int colvarproxy_namd::setup()
00187 {
00188   if (colvars->size() == 0) return COLVARS_OK;
00189 
00190   log("Updating NAMD interface:\n");
00191 
00192   if (simparams->wrapAll) {
00193     cvm::log("Warning: enabling wrapAll can lead to inconsistent results "
00194              "for Colvars calculations: please disable wrapAll, "
00195              "as is the default option in NAMD.\n");
00196   }
00197 
00198   log("updating atomic data ("+cvm::to_str(atoms_ids.size())+" atoms).\n");
00199 
00200   size_t i;
00201   for (i = 0; i < atoms_ids.size(); i++) {
00202     update_atom_properties(i);
00203 
00204     // zero out mutable arrays
00205     atoms_positions[i] = cvm::rvector(0.0, 0.0, 0.0);
00206     atoms_total_forces[i] = cvm::rvector(0.0, 0.0, 0.0);
00207     atoms_new_colvar_forces[i] = cvm::rvector(0.0, 0.0, 0.0);
00208   }
00209 
00210   size_t n_group_atoms = 0;
00211   for (int ig = 0; ig < modifyRequestedGroups().size(); ig++) {
00212     n_group_atoms += modifyRequestedGroups()[ig].size();
00213   }
00214 
00215   log("updating group data ("+cvm::to_str(atom_groups_ids.size())+" scalable groups, "+
00216       cvm::to_str(n_group_atoms)+" atoms in total).\n");
00217 
00218   // Note: groupMassBegin, groupMassEnd may be used here, but they won't work for charges
00219   for (int ig = 0; ig < modifyRequestedGroups().size(); ig++) {
00220 
00221     // update mass and charge
00222     update_group_properties(ig);
00223 
00224     atom_groups_coms[ig] = cvm::rvector(0.0, 0.0, 0.0);
00225     atom_groups_total_forces[ig] = cvm::rvector(0.0, 0.0, 0.0);
00226     atom_groups_new_colvar_forces[ig] = cvm::rvector(0.0, 0.0, 0.0);
00227   }
00228 
00229   return COLVARS_OK;
00230 }
00231 
00232 
00233 int colvarproxy_namd::reset()
00234 {
00235   int error_code = COLVARS_OK;
00236 
00237   // Unrequest all atoms and group from NAMD
00238   modifyRequestedAtoms().clear();
00239   modifyRequestedGroups().clear();
00240 
00241   atoms_map.clear();
00242 
00243   // Clear internal Proxy records
00244   error_code |= colvarproxy::reset();
00245 
00246   return error_code;
00247 }
00248 
00249 
00250 void colvarproxy_namd::calculate()
00251 {
00252   if (first_timestep) {
00253 
00254     this->setup();
00255     colvars->setup();
00256     colvars->setup_input();
00257     colvars->setup_output();
00258 
00259     first_timestep = false;
00260 
00261   } else {
00262     // Use the time step number inherited from GlobalMaster
00263     if ( step - previous_NAMD_step == 1 ) {
00264       colvars->it++;
00265     }
00266     // Other cases could mean:
00267     // - run 0
00268     // - beginning of a new run statement
00269     // then the internal counter should not be incremented
00270   }
00271 
00272   previous_NAMD_step = step;
00273 
00274   if (cvm::debug()) {
00275     log(std::string(cvm::line_marker)+
00276         "colvarproxy_namd, step no. "+cvm::to_str(colvars->it)+"\n"+
00277         "Updating atomic data arrays.\n");
00278   }
00279 
00280   // must delete the forces applied at the previous step: we can do
00281   // that because they have already been used and copied to other
00282   // memory locations
00283   modifyForcedAtoms().clear();
00284   modifyAppliedForces().clear();
00285 
00286   // prepare local arrays
00287   for (size_t i = 0; i < atoms_ids.size(); i++) {
00288     atoms_positions[i] = cvm::rvector(0.0, 0.0, 0.0);
00289     atoms_total_forces[i] = cvm::rvector(0.0, 0.0, 0.0);
00290     atoms_new_colvar_forces[i] = cvm::rvector(0.0, 0.0, 0.0);
00291   }
00292 
00293   for (size_t i = 0; i < atom_groups_ids.size(); i++) {
00294     atom_groups_total_forces[i] = cvm::rvector(0.0, 0.0, 0.0);
00295     atom_groups_new_colvar_forces[i] = cvm::rvector(0.0, 0.0, 0.0);
00296   }
00297 
00298   // create the atom map if needed
00299   size_t const n_all_atoms = Node::Object()->molecule->numAtoms;
00300   if (atoms_map.size() != n_all_atoms) {
00301     atoms_map.resize(n_all_atoms);
00302     atoms_map.assign(n_all_atoms, -1);
00303     update_atoms_map(getAtomIdBegin(), getAtomIdEnd());
00304   }
00305 
00306   // if new atomic positions or forces have been communicated by other GlobalMasters, add them to the atom map
00307   if ((int(atoms_ids.size()) < (getAtomIdEnd() - getAtomIdBegin())) ||
00308       (int(atoms_ids.size()) < (getForceIdEnd() - getForceIdBegin()))) {
00309     update_atoms_map(getAtomIdBegin(), getAtomIdEnd());
00310     update_atoms_map(getForceIdBegin(), getForceIdEnd());
00311   }
00312 
00313   {
00314     if (cvm::debug()) {
00315       log("Updating positions arrays.\n");
00316     }
00317     size_t n_positions = 0;
00318     AtomIDList::const_iterator a_i = getAtomIdBegin();
00319     AtomIDList::const_iterator a_e = getAtomIdEnd();
00320     PositionList::const_iterator p_i = getAtomPositionBegin();
00321 
00322     for ( ; a_i != a_e; ++a_i, ++p_i ) {
00323       atoms_positions[atoms_map[*a_i]] = cvm::rvector((*p_i).x, (*p_i).y, (*p_i).z);
00324       n_positions++;
00325     }
00326 
00327     // The following had to be relaxed because some atoms may be forced without their position being requested
00328     // if (n_positions < atoms_ids.size()) {
00329     //   cvm::error("Error: did not receive the positions of all atoms.\n", BUG_ERROR);
00330     // }
00331   }
00332 
00333   if (total_force_requested && cvm::step_relative() > 0) {
00334 
00335     // sort the force arrays the previous step
00336     // (but only do so if there *is* a previous step!)
00337 
00338     {
00339       if (cvm::debug()) {
00340         log("Updating total forces arrays.\n");
00341       }
00342       size_t n_total_forces = 0;
00343       AtomIDList::const_iterator a_i = getForceIdBegin();
00344       AtomIDList::const_iterator a_e = getForceIdEnd();
00345       ForceList::const_iterator f_i = getTotalForce();
00346 
00347       for ( ; a_i != a_e; ++a_i, ++f_i ) {
00348         atoms_total_forces[atoms_map[*a_i]] = cvm::rvector((*f_i).x, (*f_i).y, (*f_i).z);
00349         n_total_forces++;
00350       }
00351 
00352       if (n_total_forces < atoms_ids.size()) {
00353         cvm::error("Error: total forces were requested, but total forces "
00354                    "were not received for all atoms.\n"
00355                    "The most probable cause is combination of energy "
00356                    "minimization with a biasing method that requires MD (e.g. ABF).\n"
00357                    "Always run minimization and ABF separately.", INPUT_ERROR);
00358       }
00359     }
00360 
00361     {
00362       if (cvm::debug()) {
00363         log("Updating group total forces arrays.\n");
00364       }
00365       ForceList::const_iterator f_i = getGroupTotalForceBegin();
00366       ForceList::const_iterator f_e = getGroupTotalForceEnd();
00367       size_t i = 0;
00368       if ((f_e - f_i) != ((int) atom_groups_ids.size())) {
00369         cvm::error("Error: total forces were requested for scalable groups, "
00370                    "but they are not in the same number from the number of groups.\n"
00371                    "The most probable cause is combination of energy "
00372                    "minimization with a biasing method that requires MD (e.g. ABF).\n"
00373                    "Always run minimization and ABF separately.", INPUT_ERROR);
00374       }
00375       for ( ; f_i != f_e; f_i++, i++) {
00376         atom_groups_total_forces[i] = cvm::rvector((*f_i).x, (*f_i).y, (*f_i).z);
00377       }
00378     }
00379   }
00380 
00381   {
00382     if (cvm::debug()) {
00383       log("Updating group positions arrays.\n");
00384     }
00385     // update group data (only coms available so far)
00386     size_t ig;
00387     // note: getGroupMassBegin() could be used here, but masses and charges
00388     // have already been calculated from the last call to setup()
00389     PositionList::const_iterator gp_i = getGroupPositionBegin();
00390     for (ig = 0; gp_i != getGroupPositionEnd(); gp_i++, ig++) {
00391       atom_groups_coms[ig] = cvm::rvector(gp_i->x, gp_i->y, gp_i->z);
00392     }
00393   }
00394 
00395   if (cvm::debug()) {
00396     log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n");
00397     log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n");
00398     log("atoms_masses = "+cvm::to_str(atoms_masses)+"\n");
00399     log("atoms_charges = "+cvm::to_str(atoms_charges)+"\n");
00400     log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n");
00401     log("atoms_total_forces = "+cvm::to_str(atoms_total_forces)+"\n");
00402     log(cvm::line_marker);
00403 
00404     log("atom_groups_ids = "+cvm::to_str(atom_groups_ids)+"\n");
00405     log("atom_groups_ncopies = "+cvm::to_str(atom_groups_ncopies)+"\n");
00406     log("atom_groups_masses = "+cvm::to_str(atom_groups_masses)+"\n");
00407     log("atom_groups_charges = "+cvm::to_str(atom_groups_charges)+"\n");
00408     log("atom_groups_coms = "+cvm::to_str(atom_groups_coms)+"\n");
00409     log("atom_groups_total_forces = "+cvm::to_str(atom_groups_total_forces)+"\n");
00410     log(cvm::line_marker);
00411   }
00412 
00413   // call the collective variable module
00414   if (colvars->calc() != COLVARS_OK) {
00415     cvm::error("Error in the collective variables module.\n", COLVARS_ERROR);
00416   }
00417 
00418   if (cvm::debug()) {
00419     log(cvm::line_marker);
00420     log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n");
00421     log(cvm::line_marker);
00422     log("atom_groups_new_colvar_forces = "+cvm::to_str(atom_groups_new_colvar_forces)+"\n");
00423     log(cvm::line_marker);
00424   }
00425 
00426   // communicate all forces to the MD integrator
00427   for (size_t i = 0; i < atoms_ids.size(); i++) {
00428     cvm::rvector const &f = atoms_new_colvar_forces[i];
00429     modifyForcedAtoms().add(atoms_ids[i]);
00430     modifyAppliedForces().add(Vector(f.x, f.y, f.z));
00431   }
00432 
00433   {
00434     // zero out the applied forces on each group
00435     modifyGroupForces().resize(modifyRequestedGroups().size());
00436     ForceList::iterator gf_i = modifyGroupForces().begin();
00437     for (int ig = 0; gf_i != modifyGroupForces().end(); gf_i++, ig++) {
00438       cvm::rvector const &f = atom_groups_new_colvar_forces[ig];
00439       *gf_i = Vector(f.x, f.y, f.z);
00440     }
00441   }
00442 
00443   // send MISC energy
00444   reduction->submit();
00445 
00446   // NAMD does not destruct GlobalMaster objects, so we must remember
00447   // to write all output files at the end of a run
00448   if (step == simparams->N) {
00449     colvars->write_output_files();
00450   }
00451 }
00452 
00453 
00454 // Callback functions
00455 
00456 int colvarproxy_namd::run_force_callback() {
00457 #ifdef NAMD_TCL
00458   std::string cmd = std::string("calc_colvar_forces ")
00459     + cvm::to_str(cvm::step_absolute());
00460   int err = Tcl_Eval(interp, cmd.c_str());
00461   if (err != TCL_OK) {
00462     log(std::string("Error while executing calc_colvar_forces:\n"));
00463     cvm::error(Tcl_GetStringResult(interp));
00464     return COLVARS_ERROR;
00465   }
00466   return COLVARS_OK;
00467 #else
00468   return COLVARS_NOT_IMPLEMENTED;
00469 #endif
00470 }
00471 
00472 int colvarproxy_namd::run_colvar_callback(std::string const &name,
00473                                           std::vector<const colvarvalue *> const &cvc_values,
00474                                           colvarvalue &value)
00475 {
00476 #ifdef NAMD_TCL
00477   size_t i;
00478   std::string cmd = std::string("calc_") + name;
00479   for (i = 0; i < cvc_values.size(); i++) {
00480     cmd += std::string(" {") +  (*(cvc_values[i])).to_simple_string() + std::string("}");
00481   }
00482   int err = Tcl_Eval(interp, cmd.c_str());
00483   const char *result = Tcl_GetStringResult(interp);
00484   if (err != TCL_OK) {
00485     log(std::string("Error while executing ")
00486         + cmd + std::string(":\n"));
00487     cvm::error(result);
00488     return COLVARS_ERROR;
00489   }
00490   std::istringstream is(result);
00491   if (value.from_simple_string(is.str()) != COLVARS_OK) {
00492     log("Error parsing colvar value from script:");
00493     cvm::error(result);
00494     return COLVARS_ERROR;
00495   }
00496   return COLVARS_OK;
00497 #else
00498   return COLVARS_NOT_IMPLEMENTED;
00499 #endif
00500 }
00501 
00502 int colvarproxy_namd::run_colvar_gradient_callback(std::string const &name,
00503                                                    std::vector<const colvarvalue *> const &cvc_values,
00504                                                    std::vector<cvm::matrix2d<cvm::real> > &gradient)
00505 {
00506 #ifdef NAMD_TCL
00507   size_t i;
00508   std::string cmd = std::string("calc_") + name + "_gradient";
00509   for (i = 0; i < cvc_values.size(); i++) {
00510     cmd += std::string(" {") +  (*(cvc_values[i])).to_simple_string() + std::string("}");
00511   }
00512   int err = Tcl_Eval(interp, cmd.c_str());
00513   if (err != TCL_OK) {
00514     log(std::string("Error while executing ")
00515         + cmd + std::string(":\n"));
00516     cvm::error(Tcl_GetStringResult(interp));
00517     return COLVARS_ERROR;
00518   }
00519   Tcl_Obj **list;
00520   int n;
00521   Tcl_ListObjGetElements(interp, Tcl_GetObjResult(interp),
00522                          &n, &list);
00523   if (n != int(gradient.size())) {
00524     cvm::error("Error parsing list of gradient values from script: found "
00525                + cvm::to_str(n) + " values instead of " + cvm::to_str(gradient.size()));
00526     return COLVARS_ERROR;
00527   }
00528   for (i = 0; i < gradient.size(); i++) {
00529     std::istringstream is(Tcl_GetString(list[i]));
00530     if (gradient[i].from_simple_string(is.str()) != COLVARS_OK) {
00531       log("Gradient matrix size: " + cvm::to_str(gradient[i].size()));
00532       log("Gradient string: " + cvm::to_str(Tcl_GetString(list[i])));
00533       cvm::error("Error parsing gradient value from script");
00534       return COLVARS_ERROR;
00535     }
00536   }
00537   return (err == TCL_OK) ? COLVARS_OK : COLVARS_ERROR;
00538 #else
00539   return COLVARS_NOT_IMPLEMENTED;
00540 #endif
00541 }
00542 
00543 
00544 void colvarproxy_namd::add_energy(cvm::real energy)
00545 {
00546   reduction->item(REDUCTION_MISC_ENERGY) += energy;
00547 }
00548 
00549 void colvarproxy_namd::request_total_force(bool yesno)
00550 {
00551   if (cvm::debug()) {
00552     cvm::log("colvarproxy_namd::request_total_force()\n");
00553   }
00554   total_force_requested = yesno;
00555   requestTotalForce(total_force_requested);
00556   if (cvm::debug()) {
00557     cvm::log("colvarproxy_namd::request_total_force() end\n");
00558   }
00559 }
00560 
00561 void colvarproxy_namd::log(std::string const &message)
00562 {
00563   std::istringstream is(message);
00564   std::string line;
00565   while (std::getline(is, line))
00566     iout << "colvars: " << line << "\n";
00567   iout << endi;
00568 }
00569 
00570 void colvarproxy_namd::error(std::string const &message)
00571 {
00572   // In NAMD, all errors are fatal
00573   fatal_error(message);
00574 }
00575 
00576 
00577 void colvarproxy_namd::fatal_error(std::string const &message)
00578 {
00579   log(message);
00580   if (errno) {
00581     log(strerror(errno));
00582     NAMD_err("Error in the collective variables module");
00583   } else {
00584     NAMD_die("Error in the collective variables module: exiting.\n");
00585   }
00586 }
00587 
00588 
00589 void colvarproxy_namd::exit(std::string const &message)
00590 {
00591   log(message);
00592   BackEnd::exit();
00593 }
00594 
00595 
00596 int colvarproxy_namd::check_atom_id(int atom_number)
00597 {
00598   // NAMD's internal numbering starts from zero
00599   int const aid = (atom_number-1);
00600 
00601   if (cvm::debug())
00602     log("Adding atom "+cvm::to_str(atom_number)+
00603         " for collective variables calculation.\n");
00604 
00605   if ( (aid < 0) || (aid >= Node::Object()->molecule->numAtoms) ) {
00606     cvm::error("Error: invalid atom number specified, "+
00607                cvm::to_str(atom_number)+"\n", INPUT_ERROR);
00608     return INPUT_ERROR;
00609   }
00610 
00611   return aid;
00612 }
00613 
00614 
00615 int colvarproxy_namd::init_atom(int atom_number)
00616 {
00617   // save time by checking first whether this atom has been requested before
00618   // (this is more common than a non-valid atom number)
00619   int aid = (atom_number-1);
00620 
00621   for (size_t i = 0; i < atoms_ids.size(); i++) {
00622     if (atoms_ids[i] == aid) {
00623       // this atom id was already recorded
00624       atoms_ncopies[i] += 1;
00625       return i;
00626     }
00627   }
00628 
00629   aid = check_atom_id(atom_number);
00630 
00631   if (aid < 0) {
00632     return INPUT_ERROR;
00633   }
00634 
00635   int const index = add_atom_slot(aid);
00636   modifyRequestedAtoms().add(aid);
00637   update_atom_properties(index);
00638   return index;
00639 }
00640 
00641 
00642 int colvarproxy_namd::check_atom_id(cvm::residue_id const &residue,
00643                                     std::string const     &atom_name,
00644                                     std::string const     &segment_id)
00645 {
00646   int const aid =
00647     (segment_id.size() ?
00648      Node::Object()->molecule->get_atom_from_name(segment_id.c_str(),
00649                                                   residue,
00650                                                   atom_name.c_str()) :
00651      Node::Object()->molecule->get_atom_from_name("MAIN",
00652                                                   residue,
00653                                                   atom_name.c_str()));
00654 
00655   if (aid < 0) {
00656     // get_atom_from_name() has returned an error value
00657     cvm::error("Error: could not find atom \""+
00658                atom_name+"\" in residue "+
00659                cvm::to_str(residue)+
00660                ( (segment_id != "MAIN") ?
00661                  (", segment \""+segment_id+"\"") :
00662                  ("") )+
00663                "\n", INPUT_ERROR);
00664     return INPUT_ERROR;
00665   }
00666 
00667   return aid;
00668 }
00669 
00670 
00671 
00675 int colvarproxy_namd::init_atom(cvm::residue_id const &residue,
00676                                 std::string const     &atom_name,
00677                                 std::string const     &segment_id)
00678 {
00679   int const aid = check_atom_id(residue, atom_name, segment_id);
00680 
00681   for (size_t i = 0; i < atoms_ids.size(); i++) {
00682     if (atoms_ids[i] == aid) {
00683       // this atom id was already recorded
00684       atoms_ncopies[i] += 1;
00685       return i;
00686     }
00687   }
00688 
00689   if (cvm::debug())
00690     log("Adding atom \""+
00691         atom_name+"\" in residue "+
00692         cvm::to_str(residue)+
00693         " (index "+cvm::to_str(aid)+
00694         ") for collective variables calculation.\n");
00695 
00696   int const index = add_atom_slot(aid);
00697   modifyRequestedAtoms().add(aid);
00698   update_atom_properties(index);
00699   return index;
00700 }
00701 
00702 
00703 void colvarproxy_namd::clear_atom(int index)
00704 {
00705   colvarproxy::clear_atom(index);
00706   // TODO remove it from GlobalMaster arrays?
00707 }
00708 
00709 
00710 enum e_pdb_field {
00711   e_pdb_none,
00712   e_pdb_occ,
00713   e_pdb_beta,
00714   e_pdb_x,
00715   e_pdb_y,
00716   e_pdb_z,
00717   e_pdb_ntot
00718 };
00719 
00720 
00721 e_pdb_field pdb_field_str2enum(std::string const &pdb_field_str)
00722 {
00723   e_pdb_field pdb_field = e_pdb_none;
00724 
00725   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00726       colvarparse::to_lower_cppstr("O")) {
00727     pdb_field = e_pdb_occ;
00728   }
00729 
00730   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00731       colvarparse::to_lower_cppstr("B")) {
00732     pdb_field = e_pdb_beta;
00733   }
00734 
00735   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00736       colvarparse::to_lower_cppstr("X")) {
00737     pdb_field = e_pdb_x;
00738   }
00739 
00740   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00741       colvarparse::to_lower_cppstr("Y")) {
00742     pdb_field = e_pdb_y;
00743   }
00744 
00745   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00746       colvarparse::to_lower_cppstr("Z")) {
00747     pdb_field = e_pdb_z;
00748   }
00749 
00750   if (pdb_field == e_pdb_none) {
00751     cvm::error("Error: unsupported PDB field, \""+
00752                pdb_field_str+"\".\n", INPUT_ERROR);
00753   }
00754 
00755   return pdb_field;
00756 }
00757 
00758 
00759 int colvarproxy_namd::load_coords(char const *pdb_filename,
00760                                   std::vector<cvm::atom_pos> &pos,
00761                                   const std::vector<int> &indices,
00762                                   std::string const &pdb_field_str,
00763                                   double const pdb_field_value)
00764 {
00765   if (pdb_field_str.size() == 0 && indices.size() == 0) {
00766     cvm::error("Bug alert: either PDB field should be defined or list of "
00767                "atom IDs should be available when loading atom coordinates!\n", BUG_ERROR);
00768   }
00769 
00770   e_pdb_field pdb_field_index;
00771   bool const use_pdb_field = (pdb_field_str.size() > 0);
00772   if (use_pdb_field) {
00773     pdb_field_index = pdb_field_str2enum(pdb_field_str);
00774   }
00775 
00776   // next index to be looked up in PDB file (if list is supplied)
00777   std::vector<int>::const_iterator current_index = indices.begin();
00778 
00779   PDB *pdb = new PDB(pdb_filename);
00780   size_t const pdb_natoms = pdb->num_atoms();
00781 
00782   if (pos.size() != pdb_natoms) {
00783 
00784     bool const pos_allocated = (pos.size() > 0);
00785 
00786     size_t ipos = 0, ipdb = 0;
00787     for ( ; ipdb < pdb_natoms; ipdb++) {
00788 
00789       if (use_pdb_field) {
00790         // PDB field mode: skip atoms with wrong value in PDB field
00791         double atom_pdb_field_value = 0.0;
00792 
00793         switch (pdb_field_index) {
00794         case e_pdb_occ:
00795           atom_pdb_field_value = (pdb->atom(ipdb))->occupancy();
00796           break;
00797         case e_pdb_beta:
00798           atom_pdb_field_value = (pdb->atom(ipdb))->temperaturefactor();
00799           break;
00800         case e_pdb_x:
00801           atom_pdb_field_value = (pdb->atom(ipdb))->xcoor();
00802           break;
00803         case e_pdb_y:
00804           atom_pdb_field_value = (pdb->atom(ipdb))->ycoor();
00805           break;
00806         case e_pdb_z:
00807           atom_pdb_field_value = (pdb->atom(ipdb))->zcoor();
00808           break;
00809         default:
00810           break;
00811         }
00812 
00813         if ( (pdb_field_value) &&
00814              (atom_pdb_field_value != pdb_field_value) ) {
00815           continue;
00816         } else if (atom_pdb_field_value == 0.0) {
00817           continue;
00818         }
00819 
00820       } else {
00821         // Atom ID mode: use predefined atom IDs from the atom group
00822         if (((int) ipdb) != *current_index) {
00823           // Skip atoms not in the list
00824           continue;
00825         } else {
00826           current_index++;
00827         }
00828       }
00829 
00830       if (!pos_allocated) {
00831         pos.push_back(cvm::atom_pos(0.0, 0.0, 0.0));
00832       } else if (ipos >= pos.size()) {
00833         cvm::error("Error: the PDB file \""+
00834                    std::string(pdb_filename)+
00835                    "\" contains coordinates for "
00836                    "more atoms than needed.\n", BUG_ERROR);
00837       }
00838 
00839       pos[ipos] = cvm::atom_pos((pdb->atom(ipdb))->xcoor(),
00840                                 (pdb->atom(ipdb))->ycoor(),
00841                                 (pdb->atom(ipdb))->zcoor());
00842       ipos++;
00843       if (!use_pdb_field && current_index == indices.end())
00844         break;
00845     }
00846 
00847     if ((ipos < pos.size()) || (current_index != indices.end()))
00848       cvm::error("Error: the number of records in the PDB file \""+
00849                  std::string(pdb_filename)+
00850                  "\" does not appear to match either the total number of atoms,"+
00851                  " or the number of coordinates requested at this point("+
00852                  cvm::to_str(pos.size())+").\n", BUG_ERROR);
00853 
00854   } else {
00855 
00856     // when the PDB contains exactly the number of atoms of the array,
00857     // ignore the fields and just read coordinates
00858     for (size_t ia = 0; ia < pos.size(); ia++) {
00859       pos[ia] = cvm::atom_pos((pdb->atom(ia))->xcoor(),
00860                               (pdb->atom(ia))->ycoor(),
00861                               (pdb->atom(ia))->zcoor());
00862     }
00863   }
00864 
00865   delete pdb;
00866   return COLVARS_OK;
00867 }
00868 
00869 
00870 int colvarproxy_namd::load_atoms(char const *pdb_filename,
00871                                  cvm::atom_group &atoms,
00872                                  std::string const &pdb_field_str,
00873                                  double const pdb_field_value)
00874 {
00875   if (pdb_field_str.size() == 0)
00876     cvm::error("Error: must define which PDB field to use "
00877                "in order to define atoms from a PDB file.\n", INPUT_ERROR);
00878 
00879   PDB *pdb = new PDB(pdb_filename);
00880   size_t const pdb_natoms = pdb->num_atoms();
00881 
00882   e_pdb_field pdb_field_index = pdb_field_str2enum(pdb_field_str);
00883 
00884   for (size_t ipdb = 0; ipdb < pdb_natoms; ipdb++) {
00885 
00886     double atom_pdb_field_value = 0.0;
00887 
00888     switch (pdb_field_index) {
00889     case e_pdb_occ:
00890       atom_pdb_field_value = (pdb->atom(ipdb))->occupancy();
00891       break;
00892     case e_pdb_beta:
00893       atom_pdb_field_value = (pdb->atom(ipdb))->temperaturefactor();
00894       break;
00895     case e_pdb_x:
00896       atom_pdb_field_value = (pdb->atom(ipdb))->xcoor();
00897       break;
00898     case e_pdb_y:
00899       atom_pdb_field_value = (pdb->atom(ipdb))->ycoor();
00900       break;
00901     case e_pdb_z:
00902       atom_pdb_field_value = (pdb->atom(ipdb))->zcoor();
00903       break;
00904     default:
00905       break;
00906     }
00907 
00908     if ( (pdb_field_value) &&
00909          (atom_pdb_field_value != pdb_field_value) ) {
00910       continue;
00911     } else if (atom_pdb_field_value == 0.0) {
00912       continue;
00913     }
00914 
00915     if (atoms.is_enabled(colvardeps::f_ag_scalable)) {
00916       atoms.add_atom_id(ipdb);
00917     } else {
00918       atoms.add_atom(cvm::atom(ipdb+1));
00919     }
00920   }
00921 
00922   delete pdb;
00923   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
00924 }
00925 
00926 
00927 std::ostream * colvarproxy_namd::output_stream(std::string const &output_name,
00928                                                std::ios_base::openmode mode)
00929 {
00930   if (cvm::debug()) {
00931     cvm::log("Using colvarproxy_namd::output_stream()\n");
00932   }
00933   std::list<std::ostream *>::iterator osi  = output_files.begin();
00934   std::list<std::string>::iterator    osni = output_stream_names.begin();
00935   for ( ; osi != output_files.end(); osi++, osni++) {
00936     if (*osni == output_name) {
00937       return *osi;
00938     }
00939   }
00940   if (!(mode & (std::ios_base::app | std::ios_base::ate))) {
00941     colvarproxy::backup_file(output_name);
00942   }
00943   ofstream_namd *os = new ofstream_namd(output_name.c_str(), mode);
00944   if (!os->is_open()) {
00945     cvm::error("Error: cannot write to file \""+output_name+"\".\n",
00946                FILE_ERROR);
00947     return NULL;
00948   }
00949   output_stream_names.push_back(output_name);
00950   output_files.push_back(os);
00951   return os;
00952 }
00953 
00954 
00955 int colvarproxy_namd::flush_output_stream(std::ostream *os)
00956 {
00957   std::list<std::ostream *>::iterator osi  = output_files.begin();
00958   std::list<std::string>::iterator    osni = output_stream_names.begin();
00959   for ( ; osi != output_files.end(); osi++, osni++) {
00960     if (*osi == os) {
00961       ((ofstream_namd *) *osi)->flush();
00962       return COLVARS_OK;
00963     }
00964   }
00965   return COLVARS_ERROR;
00966 }
00967 
00968 
00969 int colvarproxy_namd::close_output_stream(std::string const &output_name)
00970 {
00971   std::list<std::ostream *>::iterator osi  = output_files.begin();
00972   std::list<std::string>::iterator    osni = output_stream_names.begin();
00973   for ( ; osi != output_files.end(); osi++, osni++) {
00974     if (*osni == output_name) {
00975       if (((ofstream_namd *) *osi)->is_open()) {
00976         ((ofstream_namd *) *osi)->close();
00977       }
00978       output_files.erase(osi);
00979       output_stream_names.erase(osni);
00980       return COLVARS_OK;
00981     }
00982   }
00983   return COLVARS_ERROR;
00984 }
00985 
00986 
00987 int colvarproxy_namd::backup_file(char const *filename)
00988 {
00989   if (std::string(filename).rfind(std::string(".colvars.state")) != std::string::npos) {
00990     NAMD_backup_file(filename, ".old");
00991   } else {
00992     NAMD_backup_file(filename, ".BAK");
00993   }
00994   return COLVARS_OK;
00995 }
00996 
00997 
00998 char *colvarproxy_namd::script_obj_to_str(unsigned char *obj)
00999 {
01000 #ifdef NAMD_TCL
01001   return Tcl_GetString(reinterpret_cast<Tcl_Obj *>(obj));
01002 #else
01003   // This is most likely not going to be executed
01004   return colvarproxy::script_obj_to_str(obj);
01005 #endif
01006 }
01007 
01008 int colvarproxy_namd::init_atom_group(std::vector<int> const &atoms_ids)
01009 {
01010   if (cvm::debug())
01011     log("Reguesting from NAMD a group of size "+cvm::to_str(atoms_ids.size())+
01012         " for collective variables calculation.\n");
01013 
01014   // Note: modifyRequestedGroups is supposed to be in sync with the colvarproxy arrays,
01015   // and to stay that way during a simulation
01016 
01017   // compare this new group to those already allocated inside GlobalMaster
01018   int ig;
01019   for (ig = 0; ig < modifyRequestedGroups().size(); ig++) {
01020     AtomIDList const &namd_group = modifyRequestedGroups()[ig];
01021     bool b_match = true;
01022 
01023     if (namd_group.size() != ((int) atoms_ids.size())) {
01024       b_match = false;
01025     } else {
01026       int ia;
01027       for (ia = 0; ia < namd_group.size(); ia++) {
01028         int const aid = atoms_ids[ia];
01029         if (namd_group[ia] != aid) {
01030           b_match = false;
01031           break;
01032         }
01033       }
01034     }
01035 
01036     if (b_match) {
01037       if (cvm::debug())
01038         log("Group was already added.\n");
01039       // this group already exists
01040       atom_groups_ncopies[ig] += 1;
01041       return ig;
01042     }
01043   }
01044 
01045   // add this group (note: the argument of add_atom_group_slot() is redundant for NAMD, and provided only for consistency)
01046   size_t const index = add_atom_group_slot(atom_groups_ids.size());
01047   modifyRequestedGroups().resize(atom_groups_ids.size());
01048   // the following is done in calculate()
01049   // modifyGroupForces().resize(atom_groups_ids.size());
01050   AtomIDList &namd_group = modifyRequestedGroups()[index];
01051   namd_group.resize(atoms_ids.size());
01052   int const n_all_atoms = Node::Object()->molecule->numAtoms;
01053   for (size_t ia = 0; ia < atoms_ids.size(); ia++) {
01054     int const aid = atoms_ids[ia];
01055     if (cvm::debug())
01056       log("Adding atom "+cvm::to_str(aid+1)+
01057           " for collective variables calculation.\n");
01058     if ( (aid < 0) || (aid >= n_all_atoms) ) {
01059       cvm::error("Error: invalid atom number specified, "+
01060                  cvm::to_str(aid+1)+"\n", INPUT_ERROR);
01061       return -1;
01062     }
01063     namd_group[ia] = aid;
01064   }
01065 
01066   update_group_properties(index);
01067 
01068   if (cvm::debug()) {
01069     log("Group has index "+cvm::to_str(index)+"\n");
01070     log("modifyRequestedGroups length = "+cvm::to_str(modifyRequestedGroups().size())+
01071         ", modifyGroupForces length = "+cvm::to_str(modifyGroupForces().size())+"\n");
01072   }
01073 
01074   return index;
01075 }
01076 
01077 
01078 void colvarproxy_namd::clear_atom_group(int index)
01079 {
01080   // do nothing, keep the NAMD arrays in sync with the colvarproxy ones
01081   colvarproxy::clear_atom_group(index);
01082 }
01083 
01084 
01085 int colvarproxy_namd::update_group_properties(int index)
01086 {
01087   AtomIDList const &namd_group = modifyRequestedGroups()[index];
01088   if (cvm::debug()) {
01089     log("Re-calculating total mass and charge for scalable group no. "+cvm::to_str(index+1)+" ("+
01090         cvm::to_str(namd_group.size())+" atoms).\n");
01091   }
01092 
01093   cvm::real total_mass = 0.0;
01094   cvm::real total_charge = 0.0;
01095   for (int i = 0; i < namd_group.size(); i++) {
01096     total_mass += Node::Object()->molecule->atommass(namd_group[i]);
01097     total_charge += Node::Object()->molecule->atomcharge(namd_group[i]);
01098   }
01099   atom_groups_masses[index] = total_mass;
01100   atom_groups_charges[index] = total_charge;
01101 
01102   if (cvm::debug()) {
01103     log("total mass = "+cvm::to_str(total_mass)+
01104         ", total charge = "+cvm::to_str(total_charge)+"\n");
01105   }
01106 
01107   return COLVARS_OK;
01108 }
01109 
01110 #if CMK_SMP && USE_CKLOOP // SMP only
01111 
01112 void calc_colvars_items_smp(int first, int last, void *result, int paramNum, void *param)
01113 {
01114   colvarproxy_namd *proxy = (colvarproxy_namd *) param;
01115   colvarmodule *cv = proxy->colvars;
01116 
01117   cvm::increase_depth();
01118   for (int i = first; i <= last; i++) {
01119     colvar *x = (*(cv->variables_active_smp()))[i];
01120     int x_item = (*(cv->variables_active_smp_items()))[i];
01121     if (cvm::debug()) {
01122       cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+cvm::to_str(proxy->smp_num_threads())+
01123                "]: calc_colvars_items_smp(), first = "+cvm::to_str(first)+
01124                ", last = "+cvm::to_str(last)+", cv = "+
01125                x->name+", cvc = "+cvm::to_str(x_item)+"\n");
01126     }
01127     x->calc_cvcs(x_item, 1);
01128   }
01129   cvm::decrease_depth();
01130 }
01131 
01132 
01133 int colvarproxy_namd::smp_colvars_loop()
01134 {
01135   colvarmodule *cv = this->colvars;
01136   CkLoop_Parallelize(calc_colvars_items_smp, 1, this,
01137                      cv->variables_active_smp()->size(),
01138                      0, cv->variables_active_smp()->size()-1);
01139   return cvm::get_error();
01140 }
01141 
01142 
01143 void calc_cv_biases_smp(int first, int last, void *result, int paramNum, void *param)
01144 {
01145   colvarproxy_namd *proxy = (colvarproxy_namd *) param;
01146   colvarmodule *cv = proxy->colvars;
01147 
01148   cvm::increase_depth();
01149   for (int i = first; i <= last; i++) {
01150     colvarbias *b = (*(cv->biases_active()))[i];
01151     if (cvm::debug()) {
01152       cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+cvm::to_str(proxy->smp_num_threads())+
01153                "]: calc_cv_biases_smp(), first = "+cvm::to_str(first)+
01154                ", last = "+cvm::to_str(last)+", bias = "+
01155                b->name+"\n");
01156     }
01157     b->update();
01158   }
01159   cvm::decrease_depth();
01160 }
01161 
01162 
01163 int colvarproxy_namd::smp_biases_loop()
01164 {
01165   colvarmodule *cv = this->colvars;
01166   CkLoop_Parallelize(calc_cv_biases_smp, 1, this,
01167                      cv->biases_active()->size(), 0, cv->biases_active()->size()-1);
01168   return cvm::get_error();
01169 }
01170 
01171 
01172 void calc_cv_scripted_forces(int paramNum, void *param)
01173 {
01174   colvarproxy_namd *proxy = (colvarproxy_namd *) param;
01175   colvarmodule *cv = proxy->colvars;
01176   if (cvm::debug()) {
01177     cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+cvm::to_str(proxy->smp_num_threads())+
01178              "]: calc_cv_scripted_forces()\n");
01179   }
01180   cv->calc_scripted_forces();
01181 }
01182 
01183 
01184 int colvarproxy_namd::smp_biases_script_loop()
01185 {
01186   colvarmodule *cv = this->colvars;
01187   CkLoop_Parallelize(calc_cv_biases_smp, 1, this,
01188                      cv->biases_active()->size(), 0, cv->biases_active()->size()-1,
01189                      1, NULL, CKLOOP_NONE,
01190                      calc_cv_scripted_forces, 1, this);
01191   return cvm::get_error();
01192 }
01193 
01194 #endif  // SMP section

Generated on Sat Sep 23 01:17:11 2017 for NAMD by  doxygen 1.4.7