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 cvm::rvector colvarproxy_namd::position_distance(cvm::atom_pos const &pos1,
00681                                                  cvm::atom_pos const &pos2)
00682   const
00683 {
00684   Position const p1(pos1.x, pos1.y, pos1.z);
00685   Position const p2(pos2.x, pos2.y, pos2.z);
00686   // return p2 - p1
00687   Vector const d = this->lattice->delta(p2, p1);
00688   return cvm::rvector(d.x, d.y, d.z);
00689 }
00690 
00691 
00692 
00693 enum e_pdb_field {
00694   e_pdb_none,
00695   e_pdb_occ,
00696   e_pdb_beta,
00697   e_pdb_x,
00698   e_pdb_y,
00699   e_pdb_z,
00700   e_pdb_ntot
00701 };
00702 
00703 
00704 e_pdb_field pdb_field_str2enum(std::string const &pdb_field_str)
00705 {
00706   e_pdb_field pdb_field = e_pdb_none;
00707 
00708   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00709       colvarparse::to_lower_cppstr("O")) {
00710     pdb_field = e_pdb_occ;
00711   }
00712 
00713   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00714       colvarparse::to_lower_cppstr("B")) {
00715     pdb_field = e_pdb_beta;
00716   }
00717 
00718   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00719       colvarparse::to_lower_cppstr("X")) {
00720     pdb_field = e_pdb_x;
00721   }
00722 
00723   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00724       colvarparse::to_lower_cppstr("Y")) {
00725     pdb_field = e_pdb_y;
00726   }
00727 
00728   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00729       colvarparse::to_lower_cppstr("Z")) {
00730     pdb_field = e_pdb_z;
00731   }
00732 
00733   if (pdb_field == e_pdb_none) {
00734     cvm::error("Error: unsupported PDB field, \""+
00735                pdb_field_str+"\".\n", INPUT_ERROR);
00736   }
00737 
00738   return pdb_field;
00739 }
00740 
00741 
00742 int colvarproxy_namd::load_coords(char const *pdb_filename,
00743                                   std::vector<cvm::atom_pos> &pos,
00744                                   const std::vector<int> &indices,
00745                                   std::string const &pdb_field_str,
00746                                   double const pdb_field_value)
00747 {
00748   if (pdb_field_str.size() == 0 && indices.size() == 0) {
00749     cvm::error("Bug alert: either PDB field should be defined or list of "
00750                "atom IDs should be available when loading atom coordinates!\n", BUG_ERROR);
00751   }
00752 
00753   e_pdb_field pdb_field_index;
00754   bool const use_pdb_field = (pdb_field_str.size() > 0);
00755   if (use_pdb_field) {
00756     pdb_field_index = pdb_field_str2enum(pdb_field_str);
00757   }
00758 
00759   // next index to be looked up in PDB file (if list is supplied)
00760   std::vector<int>::const_iterator current_index = indices.begin();
00761 
00762   PDB *pdb = new PDB(pdb_filename);
00763   size_t const pdb_natoms = pdb->num_atoms();
00764 
00765   if (pos.size() != pdb_natoms) {
00766 
00767     bool const pos_allocated = (pos.size() > 0);
00768 
00769     size_t ipos = 0, ipdb = 0;
00770     for ( ; ipdb < pdb_natoms; ipdb++) {
00771 
00772       if (use_pdb_field) {
00773         // PDB field mode: skip atoms with wrong value in PDB field
00774         double atom_pdb_field_value = 0.0;
00775 
00776         switch (pdb_field_index) {
00777         case e_pdb_occ:
00778           atom_pdb_field_value = (pdb->atom(ipdb))->occupancy();
00779           break;
00780         case e_pdb_beta:
00781           atom_pdb_field_value = (pdb->atom(ipdb))->temperaturefactor();
00782           break;
00783         case e_pdb_x:
00784           atom_pdb_field_value = (pdb->atom(ipdb))->xcoor();
00785           break;
00786         case e_pdb_y:
00787           atom_pdb_field_value = (pdb->atom(ipdb))->ycoor();
00788           break;
00789         case e_pdb_z:
00790           atom_pdb_field_value = (pdb->atom(ipdb))->zcoor();
00791           break;
00792         default:
00793           break;
00794         }
00795 
00796         if ( (pdb_field_value) &&
00797              (atom_pdb_field_value != pdb_field_value) ) {
00798           continue;
00799         } else if (atom_pdb_field_value == 0.0) {
00800           continue;
00801         }
00802 
00803       } else {
00804         // Atom ID mode: use predefined atom IDs from the atom group
00805         if (((int) ipdb) != *current_index) {
00806           // Skip atoms not in the list
00807           continue;
00808         } else {
00809           current_index++;
00810         }
00811       }
00812 
00813       if (!pos_allocated) {
00814         pos.push_back(cvm::atom_pos(0.0, 0.0, 0.0));
00815       } else if (ipos >= pos.size()) {
00816         cvm::error("Error: the PDB file \""+
00817                    std::string(pdb_filename)+
00818                    "\" contains coordinates for "
00819                    "more atoms than needed.\n", BUG_ERROR);
00820       }
00821 
00822       pos[ipos] = cvm::atom_pos((pdb->atom(ipdb))->xcoor(),
00823                                 (pdb->atom(ipdb))->ycoor(),
00824                                 (pdb->atom(ipdb))->zcoor());
00825       ipos++;
00826       if (!use_pdb_field && current_index == indices.end())
00827         break;
00828     }
00829 
00830     if ((ipos < pos.size()) || (current_index != indices.end()))
00831       cvm::error("Error: the number of records in the PDB file \""+
00832                  std::string(pdb_filename)+
00833                  "\" does not appear to match either the total number of atoms,"+
00834                  " or the number of coordinates requested at this point("+
00835                  cvm::to_str(pos.size())+").\n", BUG_ERROR);
00836 
00837   } else {
00838 
00839     // when the PDB contains exactly the number of atoms of the array,
00840     // ignore the fields and just read coordinates
00841     for (size_t ia = 0; ia < pos.size(); ia++) {
00842       pos[ia] = cvm::atom_pos((pdb->atom(ia))->xcoor(),
00843                               (pdb->atom(ia))->ycoor(),
00844                               (pdb->atom(ia))->zcoor());
00845     }
00846   }
00847 
00848   delete pdb;
00849   return COLVARS_OK;
00850 }
00851 
00852 
00853 int colvarproxy_namd::load_atoms(char const *pdb_filename,
00854                                  cvm::atom_group &atoms,
00855                                  std::string const &pdb_field_str,
00856                                  double const pdb_field_value)
00857 {
00858   if (pdb_field_str.size() == 0)
00859     cvm::error("Error: must define which PDB field to use "
00860                "in order to define atoms from a PDB file.\n", INPUT_ERROR);
00861 
00862   PDB *pdb = new PDB(pdb_filename);
00863   size_t const pdb_natoms = pdb->num_atoms();
00864 
00865   e_pdb_field pdb_field_index = pdb_field_str2enum(pdb_field_str);
00866 
00867   for (size_t ipdb = 0; ipdb < pdb_natoms; ipdb++) {
00868 
00869     double atom_pdb_field_value = 0.0;
00870 
00871     switch (pdb_field_index) {
00872     case e_pdb_occ:
00873       atom_pdb_field_value = (pdb->atom(ipdb))->occupancy();
00874       break;
00875     case e_pdb_beta:
00876       atom_pdb_field_value = (pdb->atom(ipdb))->temperaturefactor();
00877       break;
00878     case e_pdb_x:
00879       atom_pdb_field_value = (pdb->atom(ipdb))->xcoor();
00880       break;
00881     case e_pdb_y:
00882       atom_pdb_field_value = (pdb->atom(ipdb))->ycoor();
00883       break;
00884     case e_pdb_z:
00885       atom_pdb_field_value = (pdb->atom(ipdb))->zcoor();
00886       break;
00887     default:
00888       break;
00889     }
00890 
00891     if ( (pdb_field_value) &&
00892          (atom_pdb_field_value != pdb_field_value) ) {
00893       continue;
00894     } else if (atom_pdb_field_value == 0.0) {
00895       continue;
00896     }
00897 
00898     if (atoms.is_enabled(colvardeps::f_ag_scalable)) {
00899       atoms.add_atom_id(ipdb);
00900     } else {
00901       atoms.add_atom(cvm::atom(ipdb+1));
00902     }
00903   }
00904 
00905   delete pdb;
00906   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
00907 }
00908 
00909 
00910 std::ostream * colvarproxy_namd::output_stream(std::string const &output_name,
00911                                                std::ios_base::openmode mode)
00912 {
00913   if (cvm::debug()) {
00914     cvm::log("Using colvarproxy_namd::output_stream()\n");
00915   }
00916   std::list<std::ostream *>::iterator osi  = output_files.begin();
00917   std::list<std::string>::iterator    osni = output_stream_names.begin();
00918   for ( ; osi != output_files.end(); osi++, osni++) {
00919     if (*osni == output_name) {
00920       return *osi;
00921     }
00922   }
00923   if (!(mode & (std::ios_base::app | std::ios_base::ate))) {
00924     colvarproxy::backup_file(output_name);
00925   }
00926   ofstream_namd *os = new ofstream_namd(output_name.c_str(), mode);
00927   if (!os->is_open()) {
00928     cvm::error("Error: cannot write to file \""+output_name+"\".\n",
00929                FILE_ERROR);
00930     return NULL;
00931   }
00932   output_stream_names.push_back(output_name);
00933   output_files.push_back(os);
00934   return os;
00935 }
00936 
00937 
00938 int colvarproxy_namd::flush_output_stream(std::ostream *os)
00939 {
00940   std::list<std::ostream *>::iterator osi  = output_files.begin();
00941   std::list<std::string>::iterator    osni = output_stream_names.begin();
00942   for ( ; osi != output_files.end(); osi++, osni++) {
00943     if (*osi == os) {
00944       ((ofstream_namd *) *osi)->flush();
00945       return COLVARS_OK;
00946     }
00947   }
00948   return COLVARS_ERROR;
00949 }
00950 
00951 
00952 int colvarproxy_namd::close_output_stream(std::string const &output_name)
00953 {
00954   std::list<std::ostream *>::iterator osi  = output_files.begin();
00955   std::list<std::string>::iterator    osni = output_stream_names.begin();
00956   for ( ; osi != output_files.end(); osi++, osni++) {
00957     if (*osni == output_name) {
00958       if (((ofstream_namd *) *osi)->is_open()) {
00959         ((ofstream_namd *) *osi)->close();
00960       }
00961       delete *osi;
00962       output_files.erase(osi);
00963       output_stream_names.erase(osni);
00964       return COLVARS_OK;
00965     }
00966   }
00967   return COLVARS_ERROR;
00968 }
00969 
00970 
00971 int colvarproxy_namd::backup_file(char const *filename)
00972 {
00973   if (std::string(filename).rfind(std::string(".colvars.state")) != std::string::npos) {
00974     NAMD_backup_file(filename, ".old");
00975   } else {
00976     NAMD_backup_file(filename, ".BAK");
00977   }
00978   return COLVARS_OK;
00979 }
00980 
00981 
00982 char const *colvarproxy_namd::script_obj_to_str(unsigned char *obj)
00983 {
00984   if (cvm::debug()) {
00985     cvm::log("Called colvarproxy_namd::script_obj_to_str().\n");
00986   }
00987 #ifdef NAMD_TCL
00988   return Tcl_GetString(reinterpret_cast<Tcl_Obj *>(obj));
00989 #else
00990   // This is most likely not going to be executed
00991   return colvarproxy::script_obj_to_str(obj);
00992 #endif
00993 }
00994 
00995 int colvarproxy_namd::init_atom_group(std::vector<int> const &atoms_ids)
00996 {
00997   if (cvm::debug())
00998     log("Reguesting from NAMD a group of size "+cvm::to_str(atoms_ids.size())+
00999         " for collective variables calculation.\n");
01000 
01001   // Note: modifyRequestedGroups is supposed to be in sync with the colvarproxy arrays,
01002   // and to stay that way during a simulation
01003 
01004   // compare this new group to those already allocated inside GlobalMaster
01005   int ig;
01006   for (ig = 0; ig < modifyRequestedGroups().size(); ig++) {
01007     AtomIDList const &namd_group = modifyRequestedGroups()[ig];
01008     bool b_match = true;
01009 
01010     if (namd_group.size() != ((int) atoms_ids.size())) {
01011       b_match = false;
01012     } else {
01013       int ia;
01014       for (ia = 0; ia < namd_group.size(); ia++) {
01015         int const aid = atoms_ids[ia];
01016         if (namd_group[ia] != aid) {
01017           b_match = false;
01018           break;
01019         }
01020       }
01021     }
01022 
01023     if (b_match) {
01024       if (cvm::debug())
01025         log("Group was already added.\n");
01026       // this group already exists
01027       atom_groups_ncopies[ig] += 1;
01028       return ig;
01029     }
01030   }
01031 
01032   // add this group (note: the argument of add_atom_group_slot() is redundant for NAMD, and provided only for consistency)
01033   size_t const index = add_atom_group_slot(atom_groups_ids.size());
01034   modifyRequestedGroups().resize(atom_groups_ids.size());
01035   // the following is done in calculate()
01036   // modifyGroupForces().resize(atom_groups_ids.size());
01037   AtomIDList &namd_group = modifyRequestedGroups()[index];
01038   namd_group.resize(atoms_ids.size());
01039   int const n_all_atoms = Node::Object()->molecule->numAtoms;
01040   for (size_t ia = 0; ia < atoms_ids.size(); ia++) {
01041     int const aid = atoms_ids[ia];
01042     if (cvm::debug())
01043       log("Adding atom "+cvm::to_str(aid+1)+
01044           " for collective variables calculation.\n");
01045     if ( (aid < 0) || (aid >= n_all_atoms) ) {
01046       cvm::error("Error: invalid atom number specified, "+
01047                  cvm::to_str(aid+1)+"\n", INPUT_ERROR);
01048       return -1;
01049     }
01050     namd_group[ia] = aid;
01051   }
01052 
01053   update_group_properties(index);
01054 
01055   if (cvm::debug()) {
01056     log("Group has index "+cvm::to_str(index)+"\n");
01057     log("modifyRequestedGroups length = "+cvm::to_str(modifyRequestedGroups().size())+
01058         ", modifyGroupForces length = "+cvm::to_str(modifyGroupForces().size())+"\n");
01059   }
01060 
01061   return index;
01062 }
01063 
01064 
01065 void colvarproxy_namd::clear_atom_group(int index)
01066 {
01067   // do nothing, keep the NAMD arrays in sync with the colvarproxy ones
01068   colvarproxy::clear_atom_group(index);
01069 }
01070 
01071 
01072 int colvarproxy_namd::update_group_properties(int index)
01073 {
01074   AtomIDList const &namd_group = modifyRequestedGroups()[index];
01075   if (cvm::debug()) {
01076     log("Re-calculating total mass and charge for scalable group no. "+cvm::to_str(index+1)+" ("+
01077         cvm::to_str(namd_group.size())+" atoms).\n");
01078   }
01079 
01080   cvm::real total_mass = 0.0;
01081   cvm::real total_charge = 0.0;
01082   for (int i = 0; i < namd_group.size(); i++) {
01083     total_mass += Node::Object()->molecule->atommass(namd_group[i]);
01084     total_charge += Node::Object()->molecule->atomcharge(namd_group[i]);
01085   }
01086   atom_groups_masses[index] = total_mass;
01087   atom_groups_charges[index] = total_charge;
01088 
01089   if (cvm::debug()) {
01090     log("total mass = "+cvm::to_str(total_mass)+
01091         ", total charge = "+cvm::to_str(total_charge)+"\n");
01092   }
01093 
01094   return COLVARS_OK;
01095 }
01096 
01097 
01098 #if CMK_SMP && USE_CKLOOP // SMP only
01099 
01100 void calc_colvars_items_smp(int first, int last, void *result, int paramNum, void *param)
01101 {
01102   colvarproxy_namd *proxy = (colvarproxy_namd *) param;
01103   colvarmodule *cv = proxy->colvars;
01104 
01105   cvm::increase_depth();
01106   for (int i = first; i <= last; i++) {
01107     colvar *x = (*(cv->variables_active_smp()))[i];
01108     int x_item = (*(cv->variables_active_smp_items()))[i];
01109     if (cvm::debug()) {
01110       cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+cvm::to_str(proxy->smp_num_threads())+
01111                "]: calc_colvars_items_smp(), first = "+cvm::to_str(first)+
01112                ", last = "+cvm::to_str(last)+", cv = "+
01113                x->name+", cvc = "+cvm::to_str(x_item)+"\n");
01114     }
01115     x->calc_cvcs(x_item, 1);
01116   }
01117   cvm::decrease_depth();
01118 }
01119 
01120 
01121 int colvarproxy_namd::smp_colvars_loop()
01122 {
01123   colvarmodule *cv = this->colvars;
01124   CkLoop_Parallelize(calc_colvars_items_smp, 1, this,
01125                      cv->variables_active_smp()->size(),
01126                      0, cv->variables_active_smp()->size()-1);
01127   return cvm::get_error();
01128 }
01129 
01130 
01131 void calc_cv_biases_smp(int first, int last, void *result, int paramNum, void *param)
01132 {
01133   colvarproxy_namd *proxy = (colvarproxy_namd *) param;
01134   colvarmodule *cv = proxy->colvars;
01135 
01136   cvm::increase_depth();
01137   for (int i = first; i <= last; i++) {
01138     colvarbias *b = (*(cv->biases_active()))[i];
01139     if (cvm::debug()) {
01140       cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+cvm::to_str(proxy->smp_num_threads())+
01141                "]: calc_cv_biases_smp(), first = "+cvm::to_str(first)+
01142                ", last = "+cvm::to_str(last)+", bias = "+
01143                b->name+"\n");
01144     }
01145     b->update();
01146   }
01147   cvm::decrease_depth();
01148 }
01149 
01150 
01151 int colvarproxy_namd::smp_biases_loop()
01152 {
01153   colvarmodule *cv = this->colvars;
01154   CkLoop_Parallelize(calc_cv_biases_smp, 1, this,
01155                      cv->biases_active()->size(), 0, cv->biases_active()->size()-1);
01156   return cvm::get_error();
01157 }
01158 
01159 
01160 void calc_cv_scripted_forces(int paramNum, void *param)
01161 {
01162   colvarproxy_namd *proxy = (colvarproxy_namd *) param;
01163   colvarmodule *cv = proxy->colvars;
01164   if (cvm::debug()) {
01165     cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+cvm::to_str(proxy->smp_num_threads())+
01166              "]: calc_cv_scripted_forces()\n");
01167   }
01168   cv->calc_scripted_forces();
01169 }
01170 
01171 
01172 int colvarproxy_namd::smp_biases_script_loop()
01173 {
01174   colvarmodule *cv = this->colvars;
01175   CkLoop_Parallelize(calc_cv_biases_smp, 1, this,
01176                      cv->biases_active()->size(), 0, cv->biases_active()->size()-1,
01177                      1, NULL, CKLOOP_NONE,
01178                      calc_cv_scripted_forces, 1, this);
01179   return cvm::get_error();
01180 }
01181 
01182 #endif  // #if CMK_SMP && USE_CKLOOP

Generated on Sat Jul 21 01:17:11 2018 for NAMD by  doxygen 1.4.7