colvarproxy_vmd.C

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 <tcl.h>
00011 
00012 #include "VMDApp.h"
00013 #include "DrawMolecule.h"
00014 #include "MoleculeList.h"
00015 #include "Timestep.h"
00016 #include "Residue.h"
00017 #include "Inform.h"
00018 #include "utilities.h"
00019 
00020 #include "colvarmodule.h"
00021 #include "colvarscript.h"
00022 #include "colvaratoms.h"
00023 #include "colvarscript.h"
00024 #include "colvarproxy.h"
00025 #include "colvarproxy_vmd.h"
00026 
00027 
00028 namespace {
00029 // Keep pointers to relevant runtime objects
00030 VMDApp *colvars_vmd_ptr = NULL;
00031 colvarproxy_vmd *colvarproxy_vmd_ptr = NULL;
00032 }
00033 
00034 
00035 // Copy of declaration from colvarscript.cpp (this alone doesn't warrant a
00036 // new header file)
00037 // COLVARS_TCL is defined when relevant in colvarproxy_tcl.h, included via colvarproxy.h
00038 #ifdef COLVARS_TCL
00039 extern "C" int tcl_run_colvarscript_command(ClientData clientData,
00040 Tcl_Interp *interp_in,
00041 int objc, Tcl_Obj *const objv[]);
00042 #endif
00043 
00044 
00045 int tcl_colvars(ClientData clientData, Tcl_Interp *interp,
00046 int objc, Tcl_Obj *const objv[])
00047 {
00048 // Get pointer to VMD object
00049 if (colvars_vmd_ptr == NULL) {
00050 colvars_vmd_ptr = (VMDApp *) clientData;
00051 }
00052 return tcl_run_colvarscript_command(clientData, interp, objc, objv);
00053 }
00054 
00055 
00056 int tcl_colvars_vmd_init(Tcl_Interp *interp, int molid_input)
00057 {
00058 VMDApp *const vmd = colvars_vmd_ptr;
00059 int molid = molid_input == -(1<<16) ? vmd->molecule_top() : molid_input;
00060 if (vmd->molecule_valid_id(molid)) {
00061 colvarproxy_vmd_ptr = new colvarproxy_vmd(interp, vmd, molid);
00062 return (cvm::get_error() == COLVARS_OK) ? TCL_OK : TCL_ERROR;
00063 } else {
00064 Tcl_SetResult(interp, (char *) "Error: molecule not found.",
00065 TCL_STATIC);
00066 return TCL_ERROR;
00067 }
00068 }
00069 
00070 
00071 colvarproxy_vmd::colvarproxy_vmd(Tcl_Interp *interp, VMDApp *v, int molid)
00072 : vmd(v),
00073 vmdmolid(molid),
00074 #if defined(VMDTKCON)
00075 msgColvars("colvars: ", VMDCON_INFO)
00076 #else
00077 msgColvars("colvars: ")
00078 #endif
00079 {
00080 version_int = get_version_from_string(COLVARPROXY_VERSION);
00081 b_simulation_running = false;
00082 
00083 // both fields are taken from data structures already available
00084 updated_masses_ = updated_charges_ = true;
00085 
00086 // Do we have scripts?
00087 // For now colvars depend on Tcl, but this may not always be the case
00088 // in the future
00089 #if defined(VMDTCL)
00090 have_scripts = true;
00091 // Need to set this before constructing colvarmodule, which creates colvarscript object
00092 set_tcl_interp(interp);
00093 #else
00094 have_scripts = false;
00095 #endif
00096 
00097 colvars = new colvarmodule(this);
00098 cvm::log("Using VMD interface, version "+
00099 cvm::to_str(COLVARPROXY_VERSION)+".\n");
00100 
00101 colvars->cite_feature("VMD engine");
00102 colvars->cite_feature("Colvars-VMD interface (command line)");
00103 
00104 colvars->cv_traj_freq = 0; // I/O will be handled explicitly
00105 colvars->restart_out_freq = 0;
00106 cvm::rotation::monitor_crossings = false; // Avoid unnecessary error messages
00107 
00108 // Default to VMD's native unit system, but do not set the units string
00109 // to preserve the native workflow of VMD / NAMD / LAMMPS-real
00110 angstrom_value_ = 1.;
00111 kcal_mol_value_ = 1.;
00112 
00113 colvars->setup_input();
00114 colvars->setup_output();
00115 
00116 // set the same seed as in Measure.C
00117 vmd_srandom(38572111);
00118 
00119 colvarproxy_vmd::setup();
00120 }
00121 
00122 
00123 colvarproxy_vmd::~colvarproxy_vmd()
00124 {}
00125 
00126 
00127 int colvarproxy_vmd::request_deletion()
00128 {
00129 b_delete_requested = true;
00130 return COLVARS_OK;
00131 }
00132 
00133 
00134 int colvarproxy_vmd::setup()
00135 {
00136 int error_code = colvarproxy::setup();
00137 vmdmol = vmd->moleculeList->mol_from_id(vmdmolid);
00138 if (vmdmol) {
00139 set_frame(vmdmol->frame());
00140 } else {
00141 error("Error: requested molecule ("+cvm::to_str(vmdmolid)+") does not exist.\n");
00142 return COLVARS_ERROR;
00143 }
00144 
00145 return colvars->update_engine_parameters();
00146 }
00147 
00148 
00149 cvm::real colvarproxy_vmd::rand_gaussian()
00150 {
00151 return vmd_random_gaussian();
00152 }
00153 
00154 
00155 int colvarproxy_vmd::set_unit_system(std::string const &units_in, bool check_only)
00156 {
00157 // if check_only is specified, just test for compatibility
00158 // cvolvarmodule does that if new units are requested while colvars are already defined
00159 if (check_only) {
00160 if ((units != "" && units_in != units) || (units == "" && units_in != "real")) {
00161 cvm::error("Specified unit system \"" + units_in + "\" is incompatible with previous setting \""
00162 + units + "\".\nReset the Colvars Module or delete all variables to change the unit.\n");
00163 return COLVARS_ERROR;
00164 } else {
00165 return COLVARS_OK;
00166 }
00167 }
00168 
00169 if (units_in == "real") {
00170 angstrom_value_ = 1.;
00171 kcal_mol_value_ = 1.;
00172 } else if (units_in == "metal") {
00173 angstrom_value_ = 1.;
00174 kcal_mol_value_ = 0.0433641017; // eV
00175 // inverse of LAMMPS value is 1/23.060549 = .043364102
00176 } else if (units_in == "electron") {
00177 angstrom_value_ = 1.88972612; // Bohr
00178 kcal_mol_value_ = 0.00159360144; // Hartree
00179 } else if (units_in == "gromacs") {
00180 angstrom_value_ = 0.1; // nm
00181 kcal_mol_value_ = 4.184; // kJ/mol
00182 } else {
00183 cvm::error("Unknown unit system specified: \"" + units_in + "\". Supported are real, metal, electron, and gromacs.\n");
00184 return COLVARS_ERROR;
00185 }
00186 
00187 units = units_in;
00188 return COLVARS_OK;
00189 }
00190 
00191 
00192 int colvarproxy_vmd::update_input()
00193 {
00194 colvarproxy::update_input();
00195 
00196 int error_code = COLVARS_OK;
00197 
00198 // Check that our parent molecule still exists
00199 if (vmd->moleculeList->mol_from_id(vmdmolid) == NULL) {
00200 error("Error: requested molecule ("+cvm::to_str(vmdmolid)+") does not exist.\n");
00201 return COLVARS_ERROR;
00202 }
00203 error_code |= update_atomic_properties();
00204 
00205 size_t i;
00206 // We're not applying any forces but they can be tracked through [cv getatomappliedforces]
00207 // Clear before updating Module
00208 for (i = 0; i < atoms_new_colvar_forces.size(); i++) {
00209 atoms_new_colvar_forces[i].reset();
00210 }
00211 
00212 // Do we still have a valid frame?
00213 if (error_code || vmdmol->get_frame(vmdmol_frame) == NULL) {
00214 error_code |= COLVARS_NO_SUCH_FRAME;
00215 return error_code;
00216 }
00217 
00218 // copy positions in the internal arrays
00219 float *vmdpos = (vmdmol->get_frame(vmdmol_frame))->pos;
00220 for (i = 0; i < atoms_positions.size(); i++) {
00221 atoms_positions[i] = cvm::atom_pos(angstrom_to_internal(vmdpos[atoms_ids[i]*3+0]),
00222 angstrom_to_internal(vmdpos[atoms_ids[i]*3+1]),
00223 angstrom_to_internal(vmdpos[atoms_ids[i]*3+2]));
00224 }
00225 
00226 
00227 Timestep const *ts = vmdmol->get_frame(vmdmol_frame);
00228 {
00229 // Get lattice vectors
00230 float A[3];
00231 float B[3];
00232 float C[3];
00233 ts->get_transform_vectors(A, B, C);
00234 unit_cell_x.set(angstrom_to_internal(A[0]), angstrom_to_internal(A[1]), angstrom_to_internal(A[2]));
00235 unit_cell_y.set(angstrom_to_internal(B[0]), angstrom_to_internal(B[1]), angstrom_to_internal(B[2]));
00236 unit_cell_z.set(angstrom_to_internal(C[0]), angstrom_to_internal(C[1]), angstrom_to_internal(C[2]));
00237 }
00238 
00239 if (ts->a_length + ts->b_length + ts->c_length < 1.0e-6) {
00240 boundaries_type = boundaries_non_periodic;
00241 reset_pbc_lattice();
00242 } else if ((ts->a_length > 1.0e-6) &&
00243 (ts->b_length > 1.0e-6) &&
00244 (ts->c_length > 1.0e-6)) {
00245 if (((ts->alpha-90.0)*(ts->alpha-90.0)) +
00246 ((ts->beta-90.0)*(ts->beta-90.0)) +
00247 ((ts->gamma-90.0)*(ts->gamma-90.0)) < 1.0e-6) {
00248 boundaries_type = boundaries_pbc_ortho;
00249 } else {
00250 boundaries_type = boundaries_pbc_triclinic;
00251 }
00252 colvarproxy_system::update_pbc_lattice();
00253 } else {
00254 boundaries_type = boundaries_unsupported;
00255 }
00256 
00257 return error_code;
00258 }
00259 
00260 
00261 void colvarproxy_vmd::add_energy(cvm::real energy)
00262 {
00263 }
00264 
00265 
00266 int colvarproxy_vmd::update_atomic_properties()
00267 {
00268 float const *masses = vmdmol->mass();
00269 float const *charges = vmdmol->charge();
00270 
00271 int error_code = COLVARS_OK;
00272 
00273 if (masses == NULL) {
00274 error("Error: masses are undefined for the molecule being used.\n");
00275 error_code |= COLVARS_BUG_ERROR;
00276 } else {
00277 for (size_t i = 0; i < atoms_ids.size(); i++) {
00278 atoms_masses[i] = masses[atoms_ids[i]];
00279 }
00280 }
00281 
00282 if (charges == NULL) {
00283 error("Error: charges are undefined for the molecule being used.\n");
00284 error_code |= COLVARS_BUG_ERROR;
00285 } else {
00286 for (size_t i = 0; i < atoms_ids.size(); i++) {
00287 atoms_charges[i] = charges[atoms_ids[i]];
00288 }
00289 }
00290 
00291 return error_code;
00292 }
00293 
00294 
00295 void colvarproxy_vmd::request_total_force(bool yesno)
00296 {
00297 if ((yesno == true) && (total_force_requested == false)) {
00298 cvm::log("Warning: a bias requested total forces, which are undefined in VMD. "
00299 "This is only meaningful when analyzing a simulation where these were used, "
00300 "provided that a state file is loaded.\n");
00301 }
00302 total_force_requested = yesno;
00303 }
00304 
00305 
00306 void colvarproxy_vmd::log(std::string const &message)
00307 {
00308 std::istringstream is(message);
00309 std::string line;
00310 while (std::getline(is, line)) {
00311 msgColvars << line.c_str() << sendmsg;
00312 }
00313 }
00314 
00315 
00316 void colvarproxy_vmd::error(std::string const &message)
00317 {
00318 add_error_msg(message);
00319 log(message);
00320 }
00321 
00322 
00323 int colvarproxy_vmd::get_molid(int &molid)
00324 {
00325 molid = vmdmolid;
00326 return COLVARS_OK;
00327 }
00328 
00329 
00330 int colvarproxy_vmd::get_frame(long int &f)
00331 {
00332 f = vmdmol_frame;
00333 return COLVARS_OK;
00334 }
00335 
00336 
00337 int colvarproxy_vmd::set_frame(long int f)
00338 {
00339 if (vmdmol->get_frame(f) != NULL) {
00340 
00341 vmdmol_frame = f;
00342 colvars->it = f;
00343 
00344 update_input();
00345 
00346 return COLVARS_OK;
00347 } else {
00348 return COLVARS_NO_SUCH_FRAME;
00349 }
00350 }
00351 
00352 
00353 void colvarproxy_vmd::init_tcl_pointers()
00354 {
00355 #ifdef VMDTCL
00356 // Do nothing (when constructing colvarproxy), this will be set by colvarproxy_vmd constructor
00357 #else
00358 colvarproxy::init_tcl_pointers(); // Create dedicated interpreter
00359 #endif
00360 }
00361 
00362 
00363 // Callback functions
00364 
00365 #ifdef VMDTCL
00366 
00367 int colvarproxy_vmd::run_force_callback()
00368 {
00369 return colvarproxy::tcl_run_force_callback();
00370 }
00371 
00372 int colvarproxy_vmd::run_colvar_callback(
00373 std::string const &name,
00374 std::vector<const colvarvalue *> const &cvc_values,
00375 colvarvalue &value)
00376 {
00377 return colvarproxy::tcl_run_colvar_callback(name, cvc_values, value);
00378 }
00379 
00380 int colvarproxy_vmd::run_colvar_gradient_callback(
00381 std::string const &name,
00382 std::vector<const colvarvalue *> const &cvc_values,
00383 std::vector<cvm::matrix2d<cvm::real> > &gradient)
00384 {
00385 return colvarproxy::tcl_run_colvar_gradient_callback(name, cvc_values,
00386 gradient);
00387 }
00388 #endif
00389 
00390 
00391 enum e_pdb_field {
00392 e_pdb_none,
00393 e_pdb_occ,
00394 e_pdb_beta,
00395 e_pdb_x,
00396 e_pdb_y,
00397 e_pdb_z,
00398 e_pdb_ntot
00399 };
00400 
00401 e_pdb_field pdb_field_str2enum(std::string const &pdb_field_str)
00402 {
00403 e_pdb_field pdb_field = e_pdb_none;
00404 
00405 if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00406 colvarparse::to_lower_cppstr("O")) {
00407 pdb_field = e_pdb_occ;
00408 }
00409 
00410 if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00411 colvarparse::to_lower_cppstr("B")) {
00412 pdb_field = e_pdb_beta;
00413 }
00414 
00415 if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00416 colvarparse::to_lower_cppstr("X")) {
00417 pdb_field = e_pdb_x;
00418 }
00419 
00420 if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00421 colvarparse::to_lower_cppstr("Y")) {
00422 pdb_field = e_pdb_y;
00423 }
00424 
00425 if (colvarparse::to_lower_cppstr(pdb_field_str) ==
00426 colvarparse::to_lower_cppstr("Z")) {
00427 pdb_field = e_pdb_z;
00428 }
00429 
00430 if (pdb_field == e_pdb_none) {
00431 cvm::error("Error: unsupported PDB field, \""+
00432 pdb_field_str+"\".\n");
00433 }
00434 
00435 return pdb_field;
00436 }
00437 
00438 
00439 int colvarproxy_vmd::load_coords(char const *pdb_filename,
00440 std::vector<cvm::atom_pos> &pos,
00441 const std::vector<int> &indices,
00442 std::string const &pdb_field_str,
00443 double const pdb_field_value)
00444 {
00445 if (pdb_field_str.size() == 0 && indices.size() == 0) {
00446 cvm::error("Bug alert: either PDB field should be defined or list of "
00447 "atom IDs should be available when loading atom coordinates!\n",
00448 COLVARS_BUG_ERROR);
00449 return COLVARS_ERROR;
00450 }
00451 
00452 e_pdb_field pdb_field_index = e_pdb_none;
00453 bool const use_pdb_field = (pdb_field_str.size() > 0);
00454 if (use_pdb_field) {
00455 pdb_field_index = pdb_field_str2enum(pdb_field_str);
00456 }
00457 
00458 // next index to be looked up in PDB file (if list is supplied)
00459 std::vector<int>::const_iterator current_index = indices.begin();
00460 
00461 FileSpec *tmpspec = new FileSpec();
00462 tmpspec->autobonds = 0;
00463 int tmpmolid = vmd->molecule_load(-1, pdb_filename, "pdb", tmpspec);
00464 delete tmpspec;
00465 if (tmpmolid < 0) {
00466 cvm::error("Error: VMD could not read file \""+std::string(pdb_filename)+"\".\n",
00467 COLVARS_FILE_ERROR);
00468 return COLVARS_ERROR;
00469 }
00470 DrawMolecule *tmpmol = vmd->moleculeList->mol_from_id(tmpmolid);
00471 
00472 vmd->molecule_make_top(vmdmolid);
00473 size_t const pdb_natoms = tmpmol->nAtoms;
00474 
00475 if (pos.size() != pdb_natoms) {
00476 
00477 bool const pos_allocated = (pos.size() > 0);
00478 
00479 size_t ipos = 0, ipdb = 0;
00480 for ( ; ipdb < pdb_natoms; ipdb++) {
00481 
00482 if (use_pdb_field) {
00483 // PDB field mode: skip atoms with wrong value in PDB field
00484 double atom_pdb_field_value = 0.0;
00485 
00486 switch (pdb_field_index) {
00487 case e_pdb_occ:
00488 atom_pdb_field_value = (tmpmol->occupancy())[ipdb];
00489 break;
00490 case e_pdb_beta:
00491 atom_pdb_field_value = (tmpmol->beta())[ipdb];
00492 break;
00493 case e_pdb_x:
00494 atom_pdb_field_value = (tmpmol->get_frame(0)->pos)[ipdb*3];
00495 break;
00496 case e_pdb_y:
00497 atom_pdb_field_value = (tmpmol->get_frame(0)->pos)[ipdb*3+1];
00498 break;
00499 case e_pdb_z:
00500 atom_pdb_field_value = (tmpmol->get_frame(0)->pos)[ipdb*3+2];
00501 break;
00502 default:
00503 break;
00504 }
00505 
00506 if ( (pdb_field_value) &&
00507 (atom_pdb_field_value != pdb_field_value) ) {
00508 continue;
00509 } else if (atom_pdb_field_value == 0.0) {
00510 continue;
00511 }
00512 
00513 } else {
00514 // Atom ID mode: use predefined atom IDs from the atom group
00515 if (((int)ipdb) != *current_index) {
00516 // Skip atoms not in the list
00517 continue;
00518 } else {
00519 current_index++;
00520 }
00521 }
00522 
00523 if (!pos_allocated) {
00524 pos.push_back(cvm::atom_pos(0.0, 0.0, 0.0));
00525 } else if (ipos >= pos.size()) {
00526 cvm::error("Error: the PDB file \""+
00527 std::string(pdb_filename)+
00528 "\" contains coordinates for "
00529 "more atoms than needed.\n", COLVARS_INPUT_ERROR);
00530 return COLVARS_ERROR;
00531 }
00532 
00533 pos[ipos] = cvm::atom_pos(angstrom_to_internal((tmpmol->get_frame(0)->pos)[ipdb*3]),
00534 angstrom_to_internal((tmpmol->get_frame(0)->pos)[ipdb*3+1]),
00535 angstrom_to_internal((tmpmol->get_frame(0)->pos)[ipdb*3+2]));
00536 ipos++;
00537 if (!use_pdb_field && current_index == indices.end())
00538 break;
00539 }
00540 
00541 if (ipos < pos.size() || (!use_pdb_field && current_index != indices.end())) {
00542 size_t n_requested = use_pdb_field ? pos.size() : indices.size();
00543 cvm::error("Error: number of matching records in the PDB file \""+
00544 std::string(pdb_filename)+"\" ("+cvm::to_str(ipos)+
00545 ") does not match the number of requested coordinates ("+
00546 cvm::to_str(n_requested)+").\n", COLVARS_INPUT_ERROR);
00547 return COLVARS_ERROR;
00548 }
00549 } else {
00550 
00551 // when the PDB contains exactly the number of atoms of the array,
00552 // ignore the fields and just read coordinates
00553 for (size_t ia = 0; ia < pos.size(); ia++) {
00554 pos[ia] = cvm::atom_pos(angstrom_to_internal((tmpmol->get_frame(0)->pos)[ia*3]),
00555 angstrom_to_internal((tmpmol->get_frame(0)->pos)[ia*3+1]),
00556 angstrom_to_internal((tmpmol->get_frame(0)->pos)[ia*3+2]));
00557 }
00558 }
00559 
00560 vmd->molecule_delete(tmpmolid);
00561 return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
00562 }
00563 
00564 
00565 int colvarproxy_vmd::load_atoms(char const *pdb_filename,
00566 cvm::atom_group &atoms,
00567 std::string const &pdb_field_str,
00568 double const pdb_field_value)
00569 {
00570 if (pdb_field_str.size() == 0) {
00571 cvm::log("Error: must define which PDB field to use "
00572 "in order to define atoms from a PDB file.\n");
00573 cvm::set_error_bits(COLVARS_INPUT_ERROR);
00574 return COLVARS_ERROR;
00575 }
00576 
00577 FileSpec *tmpspec = new FileSpec();
00578 tmpspec->autobonds = 0;
00579 int tmpmolid = vmd->molecule_load(-1, pdb_filename, "pdb", tmpspec);
00580 delete tmpspec;
00581 
00582 if (tmpmolid < 0) {
00583 cvm::error("Error: VMD could not read file \""+std::string(pdb_filename)+"\".\n",
00584 COLVARS_FILE_ERROR);
00585 return COLVARS_ERROR;
00586 }
00587 DrawMolecule *tmpmol = vmd->moleculeList->mol_from_id(tmpmolid);
00588 
00589 vmd->molecule_make_top(vmdmolid);
00590 size_t const pdb_natoms = tmpmol->nAtoms;
00591 
00592 e_pdb_field pdb_field_index = pdb_field_str2enum(pdb_field_str);
00593 
00594 for (size_t ipdb = 0; ipdb < pdb_natoms; ipdb++) {
00595 
00596 double atom_pdb_field_value = 0.0;
00597 
00598 switch (pdb_field_index) {
00599 case e_pdb_occ:
00600 atom_pdb_field_value = (tmpmol->occupancy())[ipdb];
00601 break;
00602 case e_pdb_beta:
00603 atom_pdb_field_value = (tmpmol->beta())[ipdb];
00604 break;
00605 case e_pdb_x:
00606 atom_pdb_field_value = (tmpmol->get_frame(0)->pos)[ipdb*3];
00607 break;
00608 case e_pdb_y:
00609 atom_pdb_field_value = (tmpmol->get_frame(0)->pos)[ipdb*3+1];
00610 break;
00611 case e_pdb_z:
00612 atom_pdb_field_value = (tmpmol->get_frame(0)->pos)[ipdb*3+2];
00613 break;
00614 default:
00615 break;
00616 }
00617 
00618 if ( (pdb_field_value) &&
00619 (atom_pdb_field_value != pdb_field_value) ) {
00620 continue;
00621 } else if (atom_pdb_field_value == 0.0) {
00622 continue;
00623 }
00624 
00625 atoms.add_atom(cvm::atom(ipdb+1));
00626 }
00627 
00628 vmd->molecule_delete(tmpmolid);
00629 return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
00630 }
00631 
00632 
00633 
00634 int colvarproxy_vmd::check_atom_id(int atom_number)
00635 {
00636 // VMD internal numbering starts from zero
00637 int const aid(atom_number-1);
00638 
00639 if (cvm::debug())
00640 cvm::log("Adding atom "+cvm::to_str(aid+1)+
00641 " for collective variables calculation.\n");
00642 
00643 if ( (aid < 0) || (aid >= vmdmol->nAtoms) ) {
00644 cvm::error("Error: invalid atom number specified, "+
00645 cvm::to_str(atom_number)+"\n");
00646 return COLVARS_INPUT_ERROR;
00647 }
00648 
00649 return aid;
00650 }
00651 
00652 
00653 int colvarproxy_vmd::init_atom(int atom_number)
00654 {
00655 // save time by checking first whether this atom has been requested before
00656 // (this is more common than a non-valid atom number)
00657 int aid = (atom_number-1);
00658 
00659 for (size_t i = 0; i < atoms_ids.size(); i++) {
00660 if (atoms_ids[i] == aid) {
00661 // this atom id was already recorded
00662 atoms_refcount[i] += 1;
00663 return i;
00664 }
00665 }
00666 
00667 aid = check_atom_id(atom_number);
00668 
00669 if (aid < 0) {
00670 return COLVARS_INPUT_ERROR;
00671 }
00672 
00673 int const index = add_atom_slot(aid);
00674 
00675 float const *masses = vmdmol->mass();
00676 float const *charges = vmdmol->charge();
00677 atoms_masses[index] = masses[aid];
00678 atoms_charges[index] = charges[aid];
00679 
00680 return index;
00681 }
00682 
00683 
00684 int colvarproxy_vmd::check_atom_id(cvm::residue_id const &resid,
00685 std::string const &atom_name,
00686 std::string const &segment_id)
00687 {
00688 int aid = -1;
00689 for (int ir = 0; ir < vmdmol->nResidues; ir++) {
00690 Residue *vmdres = vmdmol->residue(ir);
00691 if (vmdres->resid == resid) {
00692 for (int ia = 0; ia < vmdres->atoms.num(); ia++) {
00693 int const resaid = vmdres->atoms[ia];
00694 std::string const sel_segname((vmdmol->segNames).name(vmdmol->atom(resaid)->segnameindex));
00695 std::string const sel_atom_name((vmdmol->atomNames).name(vmdmol->atom(resaid)->nameindex));
00696 if ( ((segment_id.size() == 0) || (segment_id == sel_segname)) &&
00697 (atom_name == sel_atom_name) ) {
00698 aid = resaid;
00699 break;
00700 }
00701 }
00702 }
00703 if (aid >= 0) break;
00704 }
00705 
00706 
00707 if (aid < 0) {
00708 cvm::error("Error: could not find atom \""+
00709 atom_name+"\" in residue "+
00710 cvm::to_str(resid)+
00711 ( (segment_id.size()) ?
00712 (", segment \""+segment_id+"\"") :
00713 ("") )+
00714 "\n", COLVARS_INPUT_ERROR);
00715 return COLVARS_INPUT_ERROR;
00716 }
00717 
00718 return aid;
00719 }
00720 
00721 
00722 int colvarproxy_vmd::init_atom(cvm::residue_id const &resid,
00723 std::string const &atom_name,
00724 std::string const &segment_id)
00725 {
00726 int const aid = check_atom_id(resid, atom_name, segment_id);
00727 
00728 for (size_t i = 0; i < atoms_ids.size(); i++) {
00729 if (atoms_ids[i] == aid) {
00730 // this atom id was already recorded
00731 atoms_refcount[i] += 1;
00732 return i;
00733 }
00734 }
00735 
00736 if (cvm::debug())
00737 cvm::log("Adding atom \""+
00738 atom_name+"\" in residue "+
00739 cvm::to_str(resid)+
00740 " (index "+cvm::to_str(aid)+
00741 ") for collective variables calculation.\n");
00742 
00743 int const index = add_atom_slot(aid);
00744 
00745 float const *masses = vmdmol->mass();
00746 float const *charges = vmdmol->charge();
00747 atoms_masses[index] = masses[aid];
00748 atoms_charges[index] = charges[aid];
00749 
00750 return index;
00751 }
00752 
00753 
00754 int colvarproxy_vmd::init_volmap_by_id(int volmap_id)
00755 {
00756 for (size_t i = 0; i < volmaps_ids.size(); i++) {
00757 if (volmaps_ids[i] == volmap_id) {
00758 // this map has already been requested
00759 volmaps_refcount[i] += 1;
00760 return i;
00761 }
00762 }
00763 
00764 int index = -1;
00765 
00766 int error_code = check_volmap_by_id(volmap_id);
00767 if (error_code == COLVARS_OK) {
00768 index = add_volmap_slot(volmap_id);
00769 }
00770 
00771 return index;
00772 }
00773 
00774 
00775 int colvarproxy_vmd::check_volmap_by_id(int volmap_id)
00776 {
00777 if ((volmap_id < 0) || (volmap_id >= vmdmol->num_volume_data())) {
00778 return cvm::error("Error: invalid numeric ID ("+cvm::to_str(volmap_id)+
00779 ") for map.\n", COLVARS_INPUT_ERROR);
00780 }
00781 return COLVARS_OK;
00782 }
00783 
00784 
00785 void colvarproxy_vmd::clear_volmap(int index)
00786 {
00787 colvarproxy::clear_volmap(index);
00788 }
00789 
00790 
00791 template<int flags>
00792 void colvarproxy_vmd::compute_voldata(VolumetricData const *voldata,
00793 cvm::atom_iter atom_begin,
00794 cvm::atom_iter atom_end,
00795 cvm::real *value,
00796 cvm::real *atom_field)
00797 {
00798 int i = 0;
00799 float coord[3], voxcoord[3], grad[3];
00800 cvm::rvector dV(0.0);
00801 cvm::atom_iter ai = atom_begin;
00802 cvm::atom_pos const origin(0.0, 0.0, 0.0);
00803 
00804 for ( ; ai != atom_end; ai++, i++) {
00805 
00806 // Wrap around the origin
00807 cvm::rvector const wrapped_pos = position_distance(origin, ai->pos);
00808 coord[0] = internal_to_angstrom(wrapped_pos.x);
00809 coord[1] = internal_to_angstrom(wrapped_pos.y);
00810 coord[2] = internal_to_angstrom(wrapped_pos.z);
00811 
00812 // Get the coordinates on the grid and check the boundaries
00813 voldata->voxel_coord_from_cartesian_coord(coord, voxcoord, 0);
00814 int const gx = static_cast<int>(voxcoord[0]);
00815 int const gy = static_cast<int>(voxcoord[1]);
00816 int const gz = static_cast<int>(voxcoord[2]);
00817 int valid_coord = 1;
00818 if ((gx < 0) || (gx >= voldata->xsize) ||
00819 (gy < 0) || (gy >= voldata->ysize) ||
00820 (gz < 0) || (gz >= voldata->zsize)) {
00821 valid_coord = 0;
00822 }
00823 
00824 cvm::real const V = valid_coord ?
00825 static_cast<cvm::real>(voldata->voxel_value_interpolate(voxcoord[0],
00826 voxcoord[1],
00827 voxcoord[2])) :
00828 0.0;
00829 
00830 if (flags & volmap_flag_gradients) {
00831 if (valid_coord) {
00832 voldata->voxel_gradient_interpolate(voxcoord, grad);
00833 // Correct the sign of the gradient
00834 dV = cvm::rvector(-1.0*grad[0], -1.0*grad[1], -1.0*grad[2]);
00835 } else {
00836 dV = cvm::rvector(0.0);
00837 }
00838 }
00839 
00840 if (flags & volmap_flag_use_atom_field) {
00841 *value += V * atom_field[i];
00842 if (flags & volmap_flag_gradients) {
00843 ai->grad += atom_field[i] * dV;
00844 }
00845 } else {
00846 *value += V;
00847 if (flags & volmap_flag_gradients) {
00848 ai->grad += dV;
00849 }
00850 }
00851 }
00852 }
00853 
00854 
00855 int colvarproxy_vmd::compute_volmap(int flags,
00856 int volmap_id,
00857 cvm::atom_iter atom_begin,
00858 cvm::atom_iter atom_end,
00859 cvm::real *value,
00860 cvm::real *atom_field)
00861 {
00862 int error_code = COLVARS_OK;
00863 VolumetricData const *voldata = vmdmol->get_volume_data(volmap_id);
00864 if (voldata != NULL) {
00865 
00866 if (flags & volmap_flag_gradients) {
00867 
00868 if (flags & volmap_flag_use_atom_field) {
00869 int const new_flags = volmap_flag_gradients |
00870 volmap_flag_use_atom_field;
00871 compute_voldata<new_flags>(voldata, atom_begin, atom_end,
00872 value, atom_field);
00873 } else {
00874 int const new_flags = volmap_flag_gradients;
00875 compute_voldata<new_flags>(voldata, atom_begin, atom_end,
00876 value, NULL);
00877 }
00878 
00879 } else {
00880 
00881 if (flags & volmap_flag_use_atom_field) {
00882 int const new_flags = volmap_flag_use_atom_field;
00883 compute_voldata<new_flags>(voldata, atom_begin, atom_end,
00884 value, atom_field);
00885 } else {
00886 int const new_flags = volmap_flag_null;
00887 compute_voldata<new_flags>(voldata, atom_begin, atom_end,
00888 value, NULL);
00889 }
00890 }
00891 } else {
00892 // Error message
00893 error_code |=
00894 const_cast<colvarproxy_vmd *>(this)->check_volmap_by_id(volmap_id);
00895 }
00896 return error_code;
00897 }