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

Generated on Tue Sep 25 01:17:11 2018 for NAMD by  doxygen 1.4.7