diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index 870079a61..e3690a4ea 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -1,741 +1,734 @@ // -*- c++ -*- #include "colvarmodule.h" #include "colvar.h" #include "colvarbias_abf.h" colvarbias_abf::colvarbias_abf(char const *key) : colvarbias(key), system_force(NULL), gradients(NULL), samples(NULL), z_gradients(NULL), z_samples(NULL), czar_gradients(NULL), last_gradients(NULL), last_samples(NULL) { } int colvarbias_abf::init(std::string const &conf) { colvarbias::init(conf); provide(f_cvb_history_dependent); // TODO relax this in case of VMD plugin if (cvm::temperature() == 0.0) cvm::log("WARNING: ABF should not be run without a thermostat or at 0 Kelvin!\n"); // ************* parsing general ABF options *********************** get_keyval_feature((colvarparse *)this, conf, "applyBias", f_cvb_apply_force, true); if (!is_enabled(f_cvb_apply_force)){ cvm::log("WARNING: ABF biases will *not* be applied!\n"); } get_keyval(conf, "updateBias", update_bias, true); if (update_bias) { enable(f_cvb_history_dependent); } else { cvm::log("WARNING: ABF biases will *not* be updated!\n"); } get_keyval(conf, "hideJacobian", hide_Jacobian, false); if (hide_Jacobian) { cvm::log("Jacobian (geometric) forces will be handled internally.\n"); } else { cvm::log("Jacobian (geometric) forces will be included in reported free energy gradients.\n"); } get_keyval(conf, "fullSamples", full_samples, 200); if ( full_samples <= 1 ) full_samples = 1; min_samples = full_samples / 2; // full_samples - min_samples >= 1 is guaranteed get_keyval(conf, "inputPrefix", input_prefix, std::vector()); get_keyval(conf, "outputFreq", output_freq, cvm::restart_out_freq); get_keyval(conf, "historyFreq", history_freq, 0); b_history_files = (history_freq > 0); // shared ABF get_keyval(conf, "shared", shared_on, false); if (shared_on) { if (!cvm::replica_enabled() || cvm::replica_num() <= 1) cvm::error("Error: shared ABF requires more than one replica."); else cvm::log("shared ABF will be applied among "+ cvm::to_str(cvm::replica_num()) + " replicas.\n"); // If shared_freq is not set, we default to output_freq get_keyval(conf, "sharedFreq", shared_freq, output_freq); } // ************* checking the associated colvars ******************* if (colvars.size() == 0) { cvm::error("Error: no collective variables specified for the ABF bias.\n"); } if (update_bias) { // Request calculation of total force (which also checks for availability) if(enable(f_cvb_get_total_force)) return cvm::get_error(); } bool b_extended = false; for (size_t i = 0; i < colvars.size(); i++) { if (colvars[i]->value().type() != colvarvalue::type_scalar) { cvm::error("Error: ABF bias can only use scalar-type variables.\n"); } colvars[i]->enable(f_cv_grid); if (hide_Jacobian) { colvars[i]->enable(f_cv_hide_Jacobian); } // If any colvar is extended-system, we need to collect the extended // system gradient if (colvars[i]->is_enabled(f_cv_extended_Lagrangian)) b_extended = true; // Here we could check for orthogonality of the Cartesian coordinates // and make it just a warning if some parameter is set? } if (get_keyval(conf, "maxForce", max_force)) { if (max_force.size() != colvars.size()) { cvm::error("Error: Number of parameters to maxForce does not match number of colvars."); } for (size_t i = 0; i < colvars.size(); i++) { if (max_force[i] < 0.0) { cvm::error("Error: maxForce should be non-negative."); } } cap_force = true; } else { cap_force = false; } bin.assign(colvars.size(), 0); force_bin.assign(colvars.size(), 0); system_force = new cvm::real [colvars.size()]; // Construct empty grids based on the colvars if (cvm::debug()) { cvm::log("Allocating count and free energy gradient grids.\n"); } samples = new colvar_grid_count(colvars); gradients = new colvar_grid_gradient(colvars); gradients->samples = samples; samples->has_parent_data = true; // Data for eABF z-based estimator if (b_extended) { z_bin.assign(colvars.size(), 0); z_samples = new colvar_grid_count(colvars); z_samples->request_actual_value(); z_gradients = new colvar_grid_gradient(colvars); z_gradients->request_actual_value(); z_gradients->samples = z_samples; z_samples->has_parent_data = true; czar_gradients = new colvar_grid_gradient(colvars); } // For shared ABF, we store a second set of grids. // This used to be only if "shared" was defined, // but now we allow calling share externally (e.g. from Tcl). last_samples = new colvar_grid_count(colvars); last_gradients = new colvar_grid_gradient(colvars); last_gradients->samples = last_samples; last_samples->has_parent_data = true; shared_last_step = -1; // If custom grids are provided, read them if ( input_prefix.size() > 0 ) { read_gradients_samples(); } cvm::log("Finished ABF setup.\n"); return COLVARS_OK; } /// Destructor colvarbias_abf::~colvarbias_abf() { if (samples) { delete samples; samples = NULL; } if (gradients) { delete gradients; gradients = NULL; } if (z_samples) { delete z_samples; z_samples = NULL; } if (z_gradients) { delete z_gradients; z_gradients = NULL; } if (czar_gradients) { delete czar_gradients; czar_gradients = NULL; } // shared ABF // We used to only do this if "shared" was defined, // but now we can call shared externally if (last_samples) { delete last_samples; last_samples = NULL; } if (last_gradients) { delete last_gradients; last_gradients = NULL; } if (system_force) { delete [] system_force; system_force = NULL; } if (cvm::n_abf_biases > 0) cvm::n_abf_biases -= 1; } /// Update the FE gradient, compute and apply biasing force /// also output data to disk if needed int colvarbias_abf::update() { if (cvm::debug()) cvm::log("Updating ABF bias " + this->name); if (cvm::step_relative() == 0) { // At first timestep, do only: // initialization stuff (file operations relying on n_abf_biases // compute current value of colvars for (size_t i = 0; i < colvars.size(); i++) { bin[i] = samples->current_bin_scalar(i); } } else { for (size_t i = 0; i < colvars.size(); i++) { bin[i] = samples->current_bin_scalar(i); } if ( update_bias && samples->index_ok(force_bin) ) { // Only if requested and within bounds of the grid... for (size_t i = 0; i < colvars.size(); i++) { // get total forces (lagging by 1 timestep) from colvars // and subtract previous ABF force - system_force[i] = colvars[i]->total_force().real_value - - colvar_forces[i].real_value; -// if (cvm::debug()) -// cvm::log("ABF System force calc: cv " + cvm::to_str(i) + -// " fs " + cvm::to_str(system_force[i]) + -// " = ft " + cvm::to_str(colvars[i]->total_force().real_value) + -// " - fa " + cvm::to_str(colvar_forces[i].real_value)); + update_system_force(i); } gradients->acc_force(force_bin, system_force); } if ( z_gradients && update_bias ) { for (size_t i = 0; i < colvars.size(); i++) { z_bin[i] = z_samples->current_bin_scalar(i); } if ( z_samples->index_ok(z_bin) ) { for (size_t i = 0; i < colvars.size(); i++) { // If we are outside the range of xi, the force has not been obtained above // the function is just an accessor, so cheap to call again anyway - system_force[i] = colvars[i]->total_force().real_value - - colvar_forces[i].real_value; + update_system_force(i); } z_gradients->acc_force(z_bin, system_force); } } } // save bin for next timestep force_bin = bin; // Reset biasing forces from previous timestep for (size_t i = 0; i < colvars.size(); i++) { colvar_forces[i].reset(); } // Compute and apply the new bias, if applicable if (is_enabled(f_cvb_apply_force) && samples->index_ok(bin)) { size_t count = samples->value(bin); cvm::real fact = 1.0; // Factor that ensures smooth introduction of the force if ( count < full_samples ) { fact = ( count < min_samples) ? 0.0 : (cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples)); } const cvm::real * grad = &(gradients->value(bin)); if ( fact != 0.0 ) { if ( (colvars.size() == 1) && colvars[0]->periodic_boundaries() ) { // Enforce a zero-mean bias on periodic, 1D coordinates // in other words: boundary condition is that the biasing potential is periodic colvar_forces[0].real_value = fact * (grad[0] / cvm::real(count) - gradients->average()); } else { for (size_t i = 0; i < colvars.size(); i++) { // subtracting the mean force (opposite of the FE gradient) means adding the gradient colvar_forces[i].real_value = fact * grad[i] / cvm::real(count); } } if (cap_force) { for (size_t i = 0; i < colvars.size(); i++) { if ( colvar_forces[i].real_value * colvar_forces[i].real_value > max_force[i] * max_force[i] ) { colvar_forces[i].real_value = (colvar_forces[i].real_value > 0 ? max_force[i] : -1.0 * max_force[i]); } } } } } // update the output prefix; TODO: move later to setup_output() function if ( cvm::n_abf_biases == 1 && cvm::n_meta_biases == 0 ) { // This is the only ABF bias output_prefix = cvm::output_prefix; } else { output_prefix = cvm::output_prefix + "." + this->name; } if (output_freq && (cvm::step_absolute() % output_freq) == 0) { if (cvm::debug()) cvm::log("ABF bias trying to write gradients and samples to disk"); write_gradients_samples(output_prefix); } if (b_history_files && (cvm::step_absolute() % history_freq) == 0) { // file already exists iff cvm::step_relative() > 0 // otherwise, backup and replace write_gradients_samples(output_prefix + ".hist", (cvm::step_relative() > 0)); } if (shared_on && shared_last_step >= 0 && cvm::step_absolute() % shared_freq == 0) { // Share gradients and samples for shared ABF. replica_share(); } // Prepare for the first sharing. if (shared_last_step < 0) { // Copy the current gradient and count values into last. last_gradients->copy_grid(*gradients); last_samples->copy_grid(*samples); shared_last_step = cvm::step_absolute(); cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+"."); } return COLVARS_OK; } int colvarbias_abf::replica_share() { int p; if ( !cvm::replica_enabled() ) { cvm::error("Error: shared ABF: No replicas.\n"); return COLVARS_ERROR; } // We must have stored the last_gradients and last_samples. if (shared_last_step < 0 ) { cvm::error("Error: shared ABF: Tried to apply shared ABF before any sampling had occurred.\n"); return COLVARS_ERROR; } // Share gradients for shared ABF. cvm::log("shared ABF: Sharing gradient and samples among replicas at step "+cvm::to_str(cvm::step_absolute()) ); // Count of data items. size_t data_n = gradients->raw_data_num(); size_t samp_start = data_n*sizeof(cvm::real); size_t msg_total = data_n*sizeof(size_t) + samp_start; char* msg_data = new char[msg_total]; if (cvm::replica_index() == 0) { // Replica 0 collects the delta gradient and count from the others. for (p = 1; p < cvm::replica_num(); p++) { // Receive the deltas. cvm::replica_comm_recv(msg_data, msg_total, p); // Map the deltas from the others into the grids. last_gradients->raw_data_in((cvm::real*)(&msg_data[0])); last_samples->raw_data_in((size_t*)(&msg_data[samp_start])); // Combine the delta gradient and count of the other replicas // with Replica 0's current state (including its delta). gradients->add_grid( *last_gradients ); samples->add_grid( *last_samples ); } // Now we must send the combined gradient to the other replicas. gradients->raw_data_out((cvm::real*)(&msg_data[0])); samples->raw_data_out((size_t*)(&msg_data[samp_start])); for (p = 1; p < cvm::replica_num(); p++) { cvm::replica_comm_send(msg_data, msg_total, p); } } else { // All other replicas send their delta gradient and count. // Calculate the delta gradient and count. last_gradients->delta_grid(*gradients); last_samples->delta_grid(*samples); // Cast the raw char data to the gradient and samples. last_gradients->raw_data_out((cvm::real*)(&msg_data[0])); last_samples->raw_data_out((size_t*)(&msg_data[samp_start])); cvm::replica_comm_send(msg_data, msg_total, 0); // We now receive the combined gradient from Replica 0. cvm::replica_comm_recv(msg_data, msg_total, 0); // We sync to the combined gradient computed by Replica 0. gradients->raw_data_in((cvm::real*)(&msg_data[0])); samples->raw_data_in((size_t*)(&msg_data[samp_start])); } // Without a barrier it's possible that one replica starts // share 2 when other replicas haven't finished share 1. cvm::replica_comm_barrier(); // Done syncing the replicas. delete[] msg_data; // Copy the current gradient and count values into last. last_gradients->copy_grid(*gradients); last_samples->copy_grid(*samples); shared_last_step = cvm::step_absolute(); return COLVARS_OK; } void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool append) { std::string samples_out_name = prefix + ".count"; std::string gradients_out_name = prefix + ".grad"; std::ios::openmode mode = (append ? std::ios::app : std::ios::out); cvm::ofstream samples_os; cvm::ofstream gradients_os; if (!append) cvm::backup_file(samples_out_name.c_str()); samples_os.open(samples_out_name.c_str(), mode); if (!samples_os.is_open()) { cvm::error("Error opening ABF samples file " + samples_out_name + " for writing"); } samples->write_multicol(samples_os); samples_os.close(); if (!append) cvm::backup_file(gradients_out_name.c_str()); gradients_os.open(gradients_out_name.c_str(), mode); if (!gradients_os.is_open()) { cvm::error("Error opening ABF gradient file " + gradients_out_name + " for writing"); } gradients->write_multicol(gradients_os); gradients_os.close(); if (colvars.size() == 1) { std::string pmf_out_name = prefix + ".pmf"; if (!append) cvm::backup_file(pmf_out_name.c_str()); cvm::ofstream pmf_os; // Do numerical integration and output a PMF pmf_os.open(pmf_out_name.c_str(), mode); if (!pmf_os.is_open()) cvm::error("Error opening pmf file " + pmf_out_name + " for writing"); gradients->write_1D_integral(pmf_os); pmf_os << std::endl; pmf_os.close(); } if (z_gradients) { // Write eABF-related quantities std::string z_samples_out_name = prefix + ".zcount"; std::string z_gradients_out_name = prefix + ".zgrad"; std::string czar_gradients_out_name = prefix + ".czar"; cvm::ofstream z_samples_os; cvm::ofstream z_gradients_os; cvm::ofstream czar_gradients_os; if (!append) cvm::backup_file(z_samples_out_name.c_str()); z_samples_os.open(z_samples_out_name.c_str(), mode); if (!z_samples_os.is_open()) { cvm::error("Error opening ABF z sample file " + z_samples_out_name + " for writing"); } z_samples->write_multicol(z_samples_os); z_samples_os.close(); if (!append) cvm::backup_file(z_gradients_out_name.c_str()); z_gradients_os.open(z_gradients_out_name.c_str(), mode); if (!z_gradients_os.is_open()) { cvm::error("Error opening ABF z gradient file " + z_gradients_out_name + " for writing"); } z_gradients->write_multicol(z_gradients_os); z_gradients_os.close(); // Calculate CZAR estimator of gradients for (std::vector ix = czar_gradients->new_index(); czar_gradients->index_ok(ix); czar_gradients->incr(ix)) { for (size_t n = 0; n < czar_gradients->multiplicity(); n++) { czar_gradients->set_value(ix, z_gradients->value_output(ix, n) - cvm::temperature() * cvm::boltzmann() * z_samples->log_gradient_finite_diff(ix, n), n); } } if (!append) cvm::backup_file(czar_gradients_out_name.c_str()); czar_gradients_os.open(czar_gradients_out_name.c_str(), mode); if (!czar_gradients_os.is_open()) { cvm::error("Error opening CZAR gradient file " + czar_gradients_out_name + " for writing"); } czar_gradients->write_multicol(czar_gradients_os); czar_gradients_os.close(); if (colvars.size() == 1) { std::string z_pmf_out_name = prefix + ".zpmf"; if (!append) cvm::backup_file(z_pmf_out_name.c_str()); cvm::ofstream z_pmf_os; // Do numerical integration and output a PMF z_pmf_os.open(z_pmf_out_name.c_str(), mode); if (!z_pmf_os.is_open()) cvm::error("Error opening z pmf file " + z_pmf_out_name + " for writing"); z_gradients->write_1D_integral(z_pmf_os); z_pmf_os << std::endl; z_pmf_os.close(); std::string czar_pmf_out_name = prefix + ".czarpmf"; if (!append) cvm::backup_file(czar_pmf_out_name.c_str()); cvm::ofstream czar_pmf_os; // Do numerical integration and output a PMF czar_pmf_os.open(czar_pmf_out_name.c_str(), mode); if (!czar_pmf_os.is_open()) cvm::error("Error opening CZAR pmf file " + czar_pmf_out_name + " for writing"); czar_gradients->write_1D_integral(czar_pmf_os); czar_pmf_os << std::endl; czar_pmf_os.close(); } } return; } // For Tcl implementation of selection rules. /// Give the total number of bins for a given bias. int colvarbias_abf::bin_num() { return samples->number_of_points(0); } /// Calculate the bin index for a given bias. int colvarbias_abf::current_bin() { return samples->current_bin_scalar(0); } /// Give the count at a given bin index. int colvarbias_abf::bin_count(int bin_index) { if (bin_index < 0 || bin_index >= bin_num()) { cvm::error("Error: Tried to get bin count from invalid bin index "+cvm::to_str(bin_index)); return -1; } std::vector ix(1,(int)bin_index); return samples->value(ix); } void colvarbias_abf::read_gradients_samples() { std::string samples_in_name, gradients_in_name, z_samples_in_name, z_gradients_in_name; for ( size_t i = 0; i < input_prefix.size(); i++ ) { samples_in_name = input_prefix[i] + ".count"; gradients_in_name = input_prefix[i] + ".grad"; z_samples_in_name = input_prefix[i] + ".zcount"; z_gradients_in_name = input_prefix[i] + ".zgrad"; // For user-provided files, the per-bias naming scheme may not apply std::ifstream is; cvm::log("Reading sample count from " + samples_in_name + " and gradients from " + gradients_in_name); is.open(samples_in_name.c_str()); if (!is.is_open()) cvm::error("Error opening ABF samples file " + samples_in_name + " for reading"); samples->read_multicol(is, true); is.close(); is.clear(); is.open(gradients_in_name.c_str()); if (!is.is_open()) cvm::error("Error opening ABF gradient file " + gradients_in_name + " for reading"); gradients->read_multicol(is, true); is.close(); if (z_gradients) { cvm::log("Reading z sample count from " + z_samples_in_name + " and z gradients from " + z_gradients_in_name); is.clear(); is.open(z_samples_in_name.c_str()); if (!is.is_open()) cvm::error("Error opening ABF z samples file " + z_samples_in_name + " for reading"); z_samples->read_multicol(is, true); is.close(); is.clear(); is.open(z_gradients_in_name.c_str()); if (!is.is_open()) cvm::error("Error opening ABF z gradient file " + z_gradients_in_name + " for reading"); z_gradients->read_multicol(is, true); is.close(); } } return; } std::ostream & colvarbias_abf::write_restart(std::ostream& os) { std::ios::fmtflags flags(os.flags()); os << "abf {\n" << " configuration {\n" << " name " << this->name << "\n"; os << " }\n"; os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format os << "\nsamples\n"; samples->write_raw(os, 8); os.flags(flags); os << "\ngradient\n"; gradients->write_raw(os, 8); if (z_gradients) { os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format os << "\nz_samples\n"; z_samples->write_raw(os, 8); os.flags(flags); os << "\nz_gradient\n"; z_gradients->write_raw(os, 8); } os << "}\n\n"; os.flags(flags); return os; } std::istream & colvarbias_abf::read_restart(std::istream& is) { if ( input_prefix.size() > 0 ) { cvm::error("ERROR: cannot provide both inputPrefix and a colvars state file.\n", INPUT_ERROR); } size_t const start_pos = is.tellg(); cvm::log("Restarting ABF bias \""+ this->name+"\".\n"); std::string key, brace, conf; if ( !(is >> key) || !(key == "abf") || !(is >> brace) || !(brace == "{") || !(is >> colvarparse::read_block("configuration", conf)) ) { cvm::log("Error: in reading restart configuration for ABF bias \""+ this->name+"\" at position "+ cvm::to_str(is.tellg())+" in stream.\n"); is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); return is; } std::string name = ""; if ( (colvarparse::get_keyval(conf, "name", name, std::string(""), colvarparse::parse_silent)) && (name != this->name) ) cvm::error("Error: in the restart file, the " "\"abf\" block has wrong name(" + name + ")\n"); if ( name == "" ) { cvm::error("Error: \"abf\" block in the restart file has no name.\n"); } if ( !(is >> key) || !(key == "samples")) { cvm::log("Error: in reading restart configuration for ABF bias \""+ this->name+"\" at position "+ cvm::to_str(is.tellg())+" in stream.\n"); is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); return is; } if (! samples->read_raw(is)) { is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); return is; } if ( !(is >> key) || !(key == "gradient")) { cvm::log("Error: in reading restart configuration for ABF bias \""+ this->name+"\" at position "+ cvm::to_str(is.tellg())+" in stream.\n"); is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); return is; } if (! gradients->read_raw(is)) { is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); return is; } if (z_gradients) { if ( !(is >> key) || !(key == "z_samples")) { cvm::log("Error: in reading restart configuration for ABF bias \""+ this->name+"\" at position "+ cvm::to_str(is.tellg())+" in stream.\n"); is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); return is; } if (! z_samples->read_raw(is)) { is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); return is; } if ( !(is >> key) || !(key == "z_gradient")) { cvm::log("Error: in reading restart configuration for ABF bias \""+ this->name+"\" at position "+ cvm::to_str(is.tellg())+" in stream.\n"); is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); return is; } if (! z_gradients->read_raw(is)) { is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); return is; } } is >> brace; if (brace != "}") { cvm::error("Error: corrupt restart information for ABF bias \""+ this->name+"\": no matching brace at position "+ cvm::to_str(is.tellg())+" in the restart file.\n"); is.setstate(std::ios::failbit); } return is; } diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h index 93079220b..a706bbb5c 100644 --- a/lib/colvars/colvarbias_abf.h +++ b/lib/colvars/colvarbias_abf.h @@ -1,96 +1,113 @@ // -*- c++ -*- #ifndef COLVARBIAS_ABF_H #define COLVARBIAS_ABF_H #include #include #include #include #include "colvarbias.h" #include "colvargrid.h" typedef cvm::real* gradient_t; /// ABF bias class colvarbias_abf : public colvarbias { public: colvarbias_abf(char const *key); virtual int init(std::string const &conf); virtual ~colvarbias_abf(); virtual int update(); private: /// Filename prefix for human-readable gradient/sample count output std::string output_prefix; /// Base filename(s) for reading previous gradient data (replaces data from restart file) std::vector input_prefix; bool update_bias; bool hide_Jacobian; size_t full_samples; size_t min_samples; /// frequency for updating output files int output_freq; /// Write combined files with a history of all output data? bool b_history_files; size_t history_freq; /// Cap applied biasing force? bool cap_force; std::vector max_force; // Internal data and methods std::vector bin, force_bin, z_bin; gradient_t system_force, applied_force; /// n-dim grid of free energy gradients colvar_grid_gradient *gradients; /// n-dim grid of number of samples colvar_grid_count *samples; /// n-dim grid: average force on "real" coordinate for eABF z-based estimator colvar_grid_gradient *z_gradients; /// n-dim grid of number of samples on "real" coordinate for eABF z-based estimator colvar_grid_count *z_samples; /// n-dim grid contining CZAR estimator of "real" free energy gradients colvar_grid_gradient *czar_gradients; + inline int update_system_force(size_t i) + { + if (colvars[i]->is_enabled(f_cv_subtract_applied_force)) { + // this colvar is already subtracting the ABF force + system_force[i] = colvars[i]->total_force().real_value; + } else { + system_force[i] = colvars[i]->total_force().real_value + - colvar_forces[i].real_value; + } + if (cvm::debug()) + cvm::log("ABF System force calc: cv " + cvm::to_str(i) + + " fs " + cvm::to_str(system_force[i]) + + " = ft " + cvm::to_str(colvars[i]->total_force().real_value) + + " - fa " + cvm::to_str(colvar_forces[i].real_value)); + return COLVARS_OK; + } + // shared ABF bool shared_on; size_t shared_freq; int shared_last_step; // Share between replicas -- may be called independently of update virtual int replica_share(); // Store the last set for shared ABF colvar_grid_gradient *last_gradients; colvar_grid_count *last_samples; // For Tcl implementation of selection rules. /// Give the total number of bins for a given bias. virtual int bin_num(); /// Calculate the bin index for a given bias. virtual int current_bin(); //// Give the count at a given bin index. virtual int bin_count(int bin_index); /// Write human-readable FE gradients and sample count void write_gradients_samples(const std::string &prefix, bool append = false); void write_last_gradients_samples(const std::string &prefix, bool append = false); /// Read human-readable FE gradients and sample count (if not using restart) void read_gradients_samples(); std::istream& read_restart(std::istream&); std::ostream& write_restart(std::ostream&); }; #endif diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index aaaa6fac0..050d4a2b9 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -1,1486 +1,1561 @@ // -*- c++ -*- #include #include #include "colvarmodule.h" #include "colvarparse.h" #include "colvarproxy.h" #include "colvar.h" #include "colvarbias.h" #include "colvarbias_abf.h" #include "colvarbias_alb.h" #include "colvarbias_histogram.h" #include "colvarbias_meta.h" #include "colvarbias_restraint.h" #include "colvarscript.h" colvarmodule::colvarmodule(colvarproxy *proxy_in) { // pointer to the proxy object if (proxy == NULL) { proxy = proxy_in; parse = new colvarparse(); } else { // TODO relax this error to handle multiple molecules in VMD // once the module is not static anymore cvm::error("Error: trying to allocate the collective " "variable module twice.\n"); return; } cvm::log(cvm::line_marker); cvm::log("Initializing the collective variables module, version "+ cvm::to_str(COLVARS_VERSION)+".\n"); cvm::log("Please cite Fiorin et al, Mol Phys 2013:\n " "http://dx.doi.org/10.1080/00268976.2013.813594\n" "in any publication based on this calculation.\n"); if (proxy->smp_enabled() == COLVARS_OK) { cvm::log("SMP parallelism is available.\n"); } // set initial default values // "it_restart" will be set by the input state file, if any; // "it" should be updated by the proxy colvarmodule::it = colvarmodule::it_restart = 0; colvarmodule::it_restart_from_state_file = true; colvarmodule::use_scripted_forces = false; colvarmodule::b_analysis = false; colvarmodule::debug_gradients_step_size = 1.0e-07; colvarmodule::rotation::monitor_crossings = false; colvarmodule::rotation::crossing_threshold = 1.0e-02; colvarmodule::cv_traj_freq = 100; colvarmodule::restart_out_freq = proxy->restart_frequency(); // by default overwrite the existing trajectory file colvarmodule::cv_traj_append = false; } int colvarmodule::read_config_file(char const *config_filename) { cvm::log(cvm::line_marker); cvm::log("Reading new configuration from file \""+ std::string(config_filename)+"\":\n"); // open the configfile config_s.open(config_filename); if (!config_s.is_open()) { cvm::error("Error: in opening configuration file \""+ std::string(config_filename)+"\".\n", FILE_ERROR); return COLVARS_ERROR; } // read the config file into a string std::string conf = ""; std::string line; while (colvarparse::getline_nocomments(config_s, line)) { // Delete lines that contain only white space after removing comments if (line.find_first_not_of(colvarparse::white_space) != std::string::npos) conf.append(line+"\n"); } config_s.close(); return parse_config(conf); } int colvarmodule::read_config_string(std::string const &config_str) { cvm::log(cvm::line_marker); cvm::log("Reading new configuration:\n"); std::istringstream config_s(config_str); // strip the comments away std::string conf = ""; std::string line; while (colvarparse::getline_nocomments(config_s, line)) { // Delete lines that contain only white space after removing comments if (line.find_first_not_of(colvarparse::white_space) != std::string::npos) conf.append(line+"\n"); } return parse_config(conf); } int colvarmodule::parse_config(std::string &conf) { // parse global options if (catch_input_errors(parse_global_params(conf))) { return get_error(); } // parse the options for collective variables if (catch_input_errors(parse_colvars(conf))) { return get_error(); } // parse the options for biases if (catch_input_errors(parse_biases(conf))) { return get_error(); } // done parsing known keywords, check that all keywords found were valid ones if (catch_input_errors(parse->check_keywords(conf, "colvarmodule"))) { return get_error(); } cvm::log(cvm::line_marker); cvm::log("Collective variables module (re)initialized.\n"); cvm::log(cvm::line_marker); // update any necessary proxy data proxy->setup(); if (cv_traj_os.is_open()) { // configuration might have changed, better redo the labels write_traj_label(cv_traj_os); } return get_error(); } int colvarmodule::parse_global_params(std::string const &conf) { std::string index_file_name; if (parse->get_keyval(conf, "indexFile", index_file_name)) { read_index_file(index_file_name.c_str()); } parse->get_keyval(conf, "analysis", b_analysis, b_analysis); parse->get_keyval(conf, "debugGradientsStepSize", debug_gradients_step_size, debug_gradients_step_size, colvarparse::parse_silent); parse->get_keyval(conf, "monitorEigenvalueCrossing", colvarmodule::rotation::monitor_crossings, colvarmodule::rotation::monitor_crossings, colvarparse::parse_silent); parse->get_keyval(conf, "eigenvalueCrossingThreshold", colvarmodule::rotation::crossing_threshold, colvarmodule::rotation::crossing_threshold, colvarparse::parse_silent); parse->get_keyval(conf, "colvarsTrajFrequency", cv_traj_freq, cv_traj_freq); parse->get_keyval(conf, "colvarsRestartFrequency", restart_out_freq, restart_out_freq); // if this is true when initializing, it means // we are continuing after a reset(): default to true parse->get_keyval(conf, "colvarsTrajAppend", cv_traj_append, cv_traj_append); parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false); parse->get_keyval(conf, "scriptingAfterBiases", scripting_after_biases, true); if (use_scripted_forces && !proxy->force_script_defined) { cvm::error("User script for scripted colvar forces not found.", INPUT_ERROR); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::parse_colvars(std::string const &conf) { if (cvm::debug()) cvm::log("Initializing the collective variables.\n"); std::string colvar_conf = ""; size_t pos = 0; while (parse->key_lookup(conf, "colvar", colvar_conf, pos)) { if (colvar_conf.size()) { cvm::log(cvm::line_marker); cvm::increase_depth(); colvars.push_back(new colvar(colvar_conf)); if (cvm::get_error() || ((colvars.back())->check_keywords(colvar_conf, "colvar") != COLVARS_OK)) { cvm::log("Error while constructing colvar number " + cvm::to_str(colvars.size()) + " : deleting."); delete colvars.back(); // the colvar destructor updates the colvars array return COLVARS_ERROR; } cvm::decrease_depth(); } else { cvm::error("Error: \"colvar\" keyword found without any configuration.\n", INPUT_ERROR); return COLVARS_ERROR; } cvm::decrease_depth(); colvar_conf = ""; } if (!colvars.size()) { cvm::log("Warning: no collective variables defined.\n"); } if (colvars.size()) cvm::log(cvm::line_marker); cvm::log("Collective variables initialized, "+ cvm::to_str(colvars.size())+ " in total.\n"); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } bool colvarmodule::check_new_bias(std::string &conf, char const *key) { if (cvm::get_error() || (biases.back()->check_keywords(conf, key) != COLVARS_OK)) { cvm::log("Error while constructing bias number " + cvm::to_str(biases.size()) + " : deleting.\n"); delete biases.back(); // the bias destructor updates the biases array return true; } return false; } template int colvarmodule::parse_biases_type(std::string const &conf, char const *keyword, size_t &bias_count) { std::string bias_conf = ""; size_t conf_saved_pos = 0; while (parse->key_lookup(conf, keyword, bias_conf, conf_saved_pos)) { if (bias_conf.size()) { cvm::log(cvm::line_marker); cvm::increase_depth(); biases.push_back(new bias_type(keyword)); biases.back()->init(bias_conf); if (cvm::check_new_bias(bias_conf, keyword) != COLVARS_OK) { return COLVARS_ERROR; } cvm::decrease_depth(); bias_count++; } else { cvm::error("Error: keyword \""+std::string(keyword)+"\" found without configuration.\n", INPUT_ERROR); return COLVARS_ERROR; } bias_conf = ""; } return COLVARS_OK; } int colvarmodule::parse_biases(std::string const &conf) { if (cvm::debug()) cvm::log("Initializing the collective variables biases.\n"); /// initialize ABF instances parse_biases_type(conf, "abf", n_abf_biases); /// initialize adaptive linear biases parse_biases_type(conf, "ALB", n_rest_biases); /// initialize harmonic restraints parse_biases_type(conf, "harmonic", n_rest_biases); /// initialize histograms parse_biases_type(conf, "histogram", n_histo_biases); /// initialize histogram restraints parse_biases_type(conf, "histogramRestraint", n_rest_biases); /// initialize linear restraints parse_biases_type(conf, "linear", n_rest_biases); /// initialize metadynamics instances parse_biases_type(conf, "metadynamics", n_meta_biases); if (use_scripted_forces) { cvm::log(cvm::line_marker); cvm::increase_depth(); cvm::log("User forces script will be run at each bias update."); cvm::decrease_depth(); } size_t i; for (i = 0; i < biases.size(); i++) { biases[i]->enable(colvardeps::f_cvb_active); if (cvm::debug()) biases[i]->print_state(); } size_t n_hist_dep_biases = 0; std::vector hist_dep_biases_names; for (i = 0; i < biases.size(); i++) { if (biases[i]->is_enabled(colvardeps::f_cvb_apply_force) && biases[i]->is_enabled(colvardeps::f_cvb_history_dependent)) { n_hist_dep_biases++; hist_dep_biases_names.push_back(biases[i]->name); } } if (n_hist_dep_biases > 1) { cvm::log("WARNING: there are "+cvm::to_str(n_hist_dep_biases)+ " history-dependent biases with non-zero force parameters:\n"+ cvm::to_str(hist_dep_biases_names)+"\n"+ "Please make sure that their forces do not counteract each other.\n"); } if (biases.size() || use_scripted_forces) { cvm::log(cvm::line_marker); cvm::log("Collective variables biases initialized, "+ cvm::to_str(biases.size())+" in total.\n"); } else { if (!use_scripted_forces) { cvm::log("No collective variables biases were defined.\n"); } } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::catch_input_errors(int result) { if (result != COLVARS_OK || get_error()) { set_error_bits(result); set_error_bits(INPUT_ERROR); parse->init(); return get_error(); } return COLVARS_OK; } colvarbias * colvarmodule::bias_by_name(std::string const &name) { for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { if ((*bi)->name == name) { return (*bi); } } return NULL; } colvar *colvarmodule::colvar_by_name(std::string const &name) { for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if ((*cvi)->name == name) { return (*cvi); } } return NULL; } int colvarmodule::change_configuration(std::string const &bias_name, std::string const &conf) { // This is deprecated; supported strategy is to delete the bias // and parse the new config cvm::increase_depth(); colvarbias *b; b = bias_by_name(bias_name); if (b == NULL) { cvm::error("Error: bias not found: " + bias_name); } b->change_configuration(conf); cvm::decrease_depth(); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } std::string colvarmodule::read_colvar(std::string const &name) { cvm::increase_depth(); colvar *c; std::stringstream ss; c = colvar_by_name(name); if (c == NULL) { cvm::fatal_error("Error: colvar not found: " + name); } ss << c->value(); cvm::decrease_depth(); return ss.str(); } cvm::real colvarmodule::energy_difference(std::string const &bias_name, std::string const &conf) { cvm::increase_depth(); colvarbias *b; cvm::real energy_diff = 0.; b = bias_by_name(bias_name); if (b == NULL) { cvm::fatal_error("Error: bias not found: " + bias_name); } energy_diff = b->energy_difference(conf); cvm::decrease_depth(); return energy_diff; } int colvarmodule::bias_current_bin(std::string const &bias_name) { cvm::increase_depth(); int ret; colvarbias *b = bias_by_name(bias_name); if (b != NULL) { ret = b->current_bin(); } else { cvm::error("Error: bias not found.\n"); ret = COLVARS_ERROR; } cvm::decrease_depth(); return ret; } int colvarmodule::bias_bin_num(std::string const &bias_name) { cvm::increase_depth(); int ret; colvarbias *b = bias_by_name(bias_name); if (b != NULL) { ret = b->bin_num(); } else { cvm::error("Error: bias not found.\n"); ret = COLVARS_ERROR; } cvm::decrease_depth(); return ret; } int colvarmodule::bias_bin_count(std::string const &bias_name, size_t bin_index) { cvm::increase_depth(); int ret; colvarbias *b = bias_by_name(bias_name); if (b != NULL) { ret = b->bin_count(bin_index); } else { cvm::error("Error: bias not found.\n"); ret = COLVARS_ERROR; } cvm::decrease_depth(); return ret; } int colvarmodule::bias_share(std::string const &bias_name) { cvm::increase_depth(); int ret; colvarbias *b = bias_by_name(bias_name); if (b != NULL) { b->replica_share(); ret = COLVARS_OK; } else { cvm::error("Error: bias not found.\n"); ret = COLVARS_ERROR; } cvm::decrease_depth(); return ret; } int colvarmodule::calc() { int error_code = COLVARS_OK; if (cvm::debug()) { cvm::log(cvm::line_marker); cvm::log("Collective variables module, step no. "+ cvm::to_str(cvm::step_absolute())+"\n"); } error_code |= calc_colvars(); // set biasing forces to zero before biases are calculated and summed over for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->reset_bias_force(); } error_code |= calc_biases(); error_code |= update_colvar_forces(); if (cvm::b_analysis) { error_code |= analyze(); } // write trajectory files, if needed if (cv_traj_freq && cv_traj_name.size()) { error_code |= write_traj_files(); } // write restart files, if needed if (restart_out_freq && restart_out_name.size()) { error_code |= write_restart_files(); } return error_code; } int colvarmodule::calc_colvars() { if (cvm::debug()) cvm::log("Calculating collective variables.\n"); // calculate collective variables and their gradients int error_code = COLVARS_OK; std::vector::iterator cvi; // Determine which colvars are active at this time step for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->feature_states[colvardeps::f_cv_active]->enabled = (step_absolute() % (*cvi)->get_time_step_factor() == 0); } // if SMP support is available, split up the work if (proxy->smp_enabled() == COLVARS_OK) { // first, calculate how much work (currently, how many active CVCs) each colvar has colvars_smp.resize(0); colvars_smp_items.resize(0); colvars_smp.reserve(colvars.size()); colvars_smp_items.reserve(colvars.size()); // set up a vector containing all components size_t num_colvar_items = 0; cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if (!(*cvi)->is_enabled()) continue; error_code |= (*cvi)->update_cvc_flags(); size_t num_items = (*cvi)->num_active_cvcs(); colvars_smp.reserve(colvars_smp.size() + num_items); colvars_smp_items.reserve(colvars_smp_items.size() + num_items); for (size_t icvc = 0; icvc < num_items; icvc++) { colvars_smp.push_back(*cvi); colvars_smp_items.push_back(icvc); } num_colvar_items += num_items; } cvm::decrease_depth(); // calculate colvar components in parallel error_code |= proxy->smp_colvars_loop(); cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if (!(*cvi)->is_enabled()) continue; error_code |= (*cvi)->collect_cvc_data(); } cvm::decrease_depth(); } else { // calculate colvars one at a time cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if (!(*cvi)->is_enabled()) continue; error_code |= (*cvi)->calc(); if (cvm::get_error()) { return COLVARS_ERROR; } } cvm::decrease_depth(); } error_code |= cvm::get_error(); return error_code; } int colvarmodule::calc_biases() { // update the biases and communicate their forces to the collective // variables if (cvm::debug() && biases.size()) cvm::log("Updating collective variable biases.\n"); std::vector::iterator bi; int error_code = COLVARS_OK; // if SMP support is available, split up the work if (proxy->smp_enabled() == COLVARS_OK) { if (use_scripted_forces && !scripting_after_biases) { // calculate biases and scripted forces in parallel error_code |= proxy->smp_biases_script_loop(); } else { // calculate biases in parallel error_code |= proxy->smp_biases_loop(); } } else { if (use_scripted_forces && !scripting_after_biases) { error_code |= calc_scripted_forces(); } cvm::increase_depth(); for (bi = biases.begin(); bi != biases.end(); bi++) { error_code |= (*bi)->update(); if (cvm::get_error()) { return COLVARS_ERROR; } } cvm::decrease_depth(); } cvm::real total_bias_energy = 0.0; for (bi = biases.begin(); bi != biases.end(); bi++) { total_bias_energy += (*bi)->get_energy(); } proxy->add_energy(total_bias_energy); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::update_colvar_forces() { int error_code = COLVARS_OK; std::vector::iterator cvi; std::vector::iterator bi; // sum the forces from all biases for each collective variable if (cvm::debug() && biases.size()) cvm::log("Collecting forces from all biases.\n"); cvm::increase_depth(); for (bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->communicate_forces(); if (cvm::get_error()) { return COLVARS_ERROR; } } cvm::decrease_depth(); if (use_scripted_forces && scripting_after_biases) { error_code |= calc_scripted_forces(); } cvm::real total_colvar_energy = 0.0; // sum up the forces for each colvar, including wall forces // and integrate any internal // equation of motion (extended system) if (cvm::debug()) cvm::log("Updating the internal degrees of freedom " "of colvars (if they have any).\n"); cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { // Here we call even inactive colvars, so they accumulate biasing forces // as well as update their extended-system dynamics total_colvar_energy += (*cvi)->update_forces_energy(); if (cvm::get_error()) { return COLVARS_ERROR; } } cvm::decrease_depth(); proxy->add_energy(total_colvar_energy); // make collective variables communicate their forces to their // coupled degrees of freedom (i.e. atoms) if (cvm::debug()) cvm::log("Communicating forces from the colvars to the atoms.\n"); cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if ((*cvi)->is_enabled(colvardeps::f_cv_gradient)) { if (!(*cvi)->is_enabled()) continue; (*cvi)->communicate_forces(); if (cvm::get_error()) { return COLVARS_ERROR; } } } cvm::decrease_depth(); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::calc_scripted_forces() { // Run user force script, if provided, // potentially adding scripted forces to the colvars int res; res = proxy->run_force_callback(); if (res == COLVARS_NOT_IMPLEMENTED) { cvm::error("Colvar forces scripts are not implemented."); return COLVARS_NOT_IMPLEMENTED; } if (res != COLVARS_OK) { cvm::error("Error running user colvar forces script"); return COLVARS_ERROR; } return COLVARS_OK; } int colvarmodule::write_restart_files() { if ( (cvm::step_relative() > 0) && ((cvm::step_absolute() % restart_out_freq) == 0) ) { cvm::log("Writing the state file \""+ restart_out_name+"\".\n"); proxy->backup_file(restart_out_name.c_str()); restart_out_os.open(restart_out_name.c_str()); if (!restart_out_os.is_open() || !write_restart(restart_out_os)) cvm::error("Error: in writing restart file.\n"); restart_out_os.close(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::write_traj_files() { if (!cv_traj_os.is_open()) { open_traj_file(cv_traj_name); } // write labels in the traj file every 1000 lines and at first timestep if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) { write_traj_label(cv_traj_os); } if ((cvm::step_absolute() % cv_traj_freq) == 0) { write_traj(cv_traj_os); } if (restart_out_freq && cv_traj_os.is_open()) { // flush the trajectory file if we are at the restart frequency if ( (cvm::step_relative() > 0) && ((cvm::step_absolute() % restart_out_freq) == 0) ) { cvm::log("Synchronizing (emptying the buffer of) trajectory file \""+ cv_traj_name+"\".\n"); cv_traj_os.flush(); } } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::analyze() { if (cvm::debug()) { cvm::log("colvarmodule::analyze(), step = "+cvm::to_str(it)+".\n"); } if (cvm::step_relative() == 0) cvm::log("Performing analysis.\n"); // perform colvar-specific analysis for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { cvm::increase_depth(); (*cvi)->analyze(); cvm::decrease_depth(); } // perform bias-specific analysis for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { cvm::increase_depth(); (*bi)->analyze(); cvm::decrease_depth(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::setup() { // loop over all components of all colvars to reset masses of all groups for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->setup(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } colvarmodule::~colvarmodule() { if ((proxy->smp_thread_id() == COLVARS_NOT_IMPLEMENTED) || (proxy->smp_thread_id() == 0)) { reset(); delete parse; parse = NULL; proxy = NULL; } } int colvarmodule::reset() { parse->init(); cvm::log("Resetting the Collective Variables Module.\n"); // Iterate backwards because we are deleting the elements as we go for (std::vector::reverse_iterator bi = biases.rbegin(); bi != biases.rend(); bi++) { delete *bi; // the bias destructor updates the biases array } biases.clear(); // Iterate backwards because we are deleting the elements as we go for (std::vector::reverse_iterator cvi = colvars.rbegin(); cvi != colvars.rend(); cvi++) { delete *cvi; // the colvar destructor updates the colvars array } colvars.clear(); index_groups.clear(); index_group_names.clear(); if (cv_traj_os.is_open()) { // Do not close file here, as we might not be done with it yet. cv_traj_os.flush(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::setup_input() { // name of input state file restart_in_name = proxy->input_prefix().size() ? std::string(proxy->input_prefix()+".colvars.state") : std::string("") ; // read the restart configuration, if available if (restart_in_name.size()) { // read the restart file std::ifstream input_is(restart_in_name.c_str()); if (!input_is.good()) { cvm::error("Error: in opening restart file \""+ std::string(restart_in_name)+"\".\n", FILE_ERROR); return COLVARS_ERROR; } else { cvm::log("Restarting from file \""+restart_in_name+"\".\n"); read_restart(input_is); if (cvm::get_error() != COLVARS_OK) { return COLVARS_ERROR; } cvm::log(cvm::line_marker); } } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::setup_output() { int error_code = 0; // output state file (restart) restart_out_name = proxy->restart_output_prefix().size() ? std::string(proxy->restart_output_prefix()+".colvars.state") : std::string(""); if (restart_out_name.size()) { cvm::log("The restart output state file will be \""+restart_out_name+"\".\n"); } output_prefix = proxy->output_prefix(); if (output_prefix.size()) { cvm::log("The final output state file will be \""+ (output_prefix.size() ? std::string(output_prefix+".colvars.state") : std::string("colvars.state"))+"\".\n"); // cvm::log (cvm::line_marker); } cv_traj_name = (output_prefix.size() ? std::string(output_prefix+".colvars.traj") : std::string("")); if (cv_traj_freq && cv_traj_name.size()) { error_code |= open_traj_file(cv_traj_name); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { error_code |= (*bi)->setup_output(); } if (error_code != COLVARS_OK || cvm::get_error()) { set_error_bits(FILE_ERROR); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } std::istream & colvarmodule::read_restart(std::istream &is) { + bool warn_total_forces = false; + { // read global restart information std::string restart_conf; if (is >> colvarparse::read_block("configuration", restart_conf)) { if (it_restart_from_state_file) { parse->get_keyval(restart_conf, "step", it_restart, (size_t) 0, colvarparse::parse_silent); it = it_restart; } std::string restart_version; parse->get_keyval(restart_conf, "version", restart_version, std::string(""), colvarparse::parse_silent); if (restart_version.size() && (restart_version != std::string(COLVARS_VERSION))) { cvm::log("This state file was generated with version "+restart_version+"\n"); } + if ((restart_version.size() == 0) || (restart_version.compare(std::string(COLVARS_VERSION)) < 0)) { + // check for total force change + if (proxy->total_forces_enabled()) { + warn_total_forces = true; + } + } } is.clear(); parse->clear_keyword_registry(); } // colvars restart cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if ( !((*cvi)->read_restart(is)) ) { cvm::error("Error: in reading restart configuration for collective variable \""+ (*cvi)->name+"\".\n", INPUT_ERROR); } } // biases restart for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { if (!((*bi)->read_restart(is))) { cvm::error("Error: in reading restart configuration for bias \""+ (*bi)->name+"\".\n", INPUT_ERROR); } } cvm::decrease_depth(); + if (warn_total_forces) { + cvm::log(cvm::line_marker); + cvm::log("WARNING:\n"); + std::string const warning("### CHANGES IN THE DEFINITION OF SYSTEM FORCES (NOW TOTAL FORCES)\n\ +\n\ +Starting from the version 2016-08-10 of the Colvars module, \n\ +the role of system forces has been replaced by total forces.\n\ +\n\ +These include *all* forces acting on a collective variable, whether they\n\ +come from the force field potential or from external terms\n\ +(e.g. restraints), including forces applied by Colvars itself.\n\ +\n\ +In NAMD, forces applied by Colvars, IMD, SMD, TMD, symmetry\n\ +restraints and tclForces are now all counted in the total force.\n\ +\n\ +In LAMMPS, forces applied by Colvars itself are now counted in the total\n\ +force (all forces from other fixes were being counted already).\n\ +\n\ +\n\ +### WHEN IS THIS CHANGE RELEVANT\n\ +\n\ +This change affects results *only* when (1) outputSystemForce is\n\ +requested or (2) the ABF bias is used. All other usage cases are\n\ +*unaffected* (colvar restraints, metadynamics, etc).\n\ +\n\ +When system forces are reported (flag: outputSystemForce), their values\n\ +in the output may change, but the physical trajectory is never affected.\n\ +The physical results of ABF calculations may be affected in some cases.\n\ +\n\ +\n\ +### CHANGES TO ABF CALCULATIONS\n\ +\n\ +Compared to previous Colvars versions, the ABF method will now attempt\n\ +to cancel external forces (for example, boundary walls) and it may be\n\ +not possible to resume through a state file a simulation that was\n\ +performed with a previous version.\n\ +\n\ +There are three possible scenarios:\n\ +\n\ +1. No external forces are applied to the atoms used by ABF: results are\n\ +unchanged.\n\ +\n\ +2. Some of the atoms used by ABF experience external forces, but these\n\ +forces are not applied directly to the variables used by ABF\n\ +(e.g. another colvar that uses the same atoms, tclForces, etc): in this\n\ +case, we recommend beginning a new simulation.\n\ +\n\ +3. External forces are applied to one or more of the colvars used by\n\ +ABF, but no other forces are applied to their atoms: you may use the\n\ +subtractAppliedForce keyword inside the corresponding colvars to\n\ +continue the previous simulation.\n\n"); + cvm::log(warning); + cvm::log(cvm::line_marker); + + // update this ahead of time in this special case + output_prefix = proxy->output_prefix(); + cvm::log("All output files will now be saved with the prefix \""+output_prefix+".tmp.*\".\n"); + output_prefix = output_prefix+".tmp"; + write_output_files(); + cvm::log(cvm::line_marker); + cvm::log("Please review the important warning above. After that, you may rename:\n\ +\""+output_prefix+".tmp.colvars.state\"\n\ +to:\n\ +\""+output_prefix+".colvars.state\"\n"); + cvm::error("Exiting with error until issue is addressed.\n", FATAL_ERROR); + } + return is; } int colvarmodule::backup_file(char const *filename) { return proxy->backup_file(filename); } int colvarmodule::write_output_files() { // if this is a simulation run (i.e. not a postprocessing), output data // must be written to be able to restart the simulation std::string const out_name = (output_prefix.size() ? std::string(output_prefix+".colvars.state") : std::string("colvars.state")); cvm::log("Saving collective variables state to \""+out_name+"\".\n"); std::ostream * os = proxy->output_stream(out_name); os->setf(std::ios::scientific, std::ios::floatfield); this->write_restart(*os); proxy->close_output_stream(out_name); cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_output_files(); } cvm::decrease_depth(); cvm::increase_depth(); for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->write_output_files(); } cvm::decrease_depth(); if (cv_traj_os.is_open()) { // do not close to avoid problems with multiple NAMD runs cv_traj_os.flush(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::read_traj(char const *traj_filename, long traj_read_begin, long traj_read_end) { cvm::log("Opening trajectory file \""+ std::string(traj_filename)+"\".\n"); std::ifstream traj_is(traj_filename); while (true) { while (true) { std::string line(""); do { if (!colvarparse::getline_nocomments(traj_is, line)) { cvm::log("End of file \""+std::string(traj_filename)+ "\" reached, or corrupted file.\n"); traj_is.close(); return false; } } while (line.find_first_not_of(colvarparse::white_space) == std::string::npos); std::istringstream is(line); if (!(is >> it)) return false; if ( (it < traj_read_begin) ) { if ((it % 1000) == 0) std::cerr << "Skipping trajectory step " << it << " \r"; continue; } else { if ((it % 1000) == 0) std::cerr << "Reading from trajectory, step = " << it << " \r"; if ( (traj_read_end > traj_read_begin) && (it > traj_read_end) ) { std::cerr << "\n"; cvm::error("Reached the end of the trajectory, " "read_end = "+cvm::to_str(traj_read_end)+"\n", FILE_ERROR); return COLVARS_ERROR; } for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if (!(*cvi)->read_traj(is)) { cvm::error("Error: in reading colvar \""+(*cvi)->name+ "\" from trajectory file \""+ std::string(traj_filename)+"\".\n", FILE_ERROR); return COLVARS_ERROR; } } break; } } } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } std::ostream & colvarmodule::write_restart(std::ostream &os) { os.setf(std::ios::scientific, std::ios::floatfield); os << "configuration {\n" << " step " << std::setw(it_width) << it << "\n" << " dt " << dt() << "\n" << " version " << std::string(COLVARS_VERSION) << "\n" << "}\n\n"; cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_restart(os); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->write_restart(os); } cvm::decrease_depth(); return os; } int colvarmodule::open_traj_file(std::string const &file_name) { if (cv_traj_os.is_open()) { return COLVARS_OK; } // (re)open trajectory file if (cv_traj_append) { cvm::log("Appending to colvar trajectory file \""+file_name+ "\".\n"); cv_traj_os.open(file_name.c_str(), std::ios::app); } else { cvm::log("Writing to colvar trajectory file \""+file_name+ "\".\n"); proxy->backup_file(file_name.c_str()); cv_traj_os.open(file_name.c_str()); } if (!cv_traj_os.is_open()) { cvm::error("Error: cannot write to file \""+file_name+"\".\n", FILE_ERROR); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::close_traj_file() { if (cv_traj_os.is_open()) { cv_traj_os.close(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } std::ostream & colvarmodule::write_traj_label(std::ostream &os) { if (!os.good()) { cvm::error("Cannot write to trajectory file."); return os; } os.setf(std::ios::scientific, std::ios::floatfield); os << "# " << cvm::wrap_string("step", cvm::it_width-2) << " "; cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_traj_label(os); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->write_traj_label(os); } os << "\n"; if (cvm::debug()) { os.flush(); } cvm::decrease_depth(); return os; } std::ostream & colvarmodule::write_traj(std::ostream &os) { os.setf(std::ios::scientific, std::ios::floatfield); os << std::setw(cvm::it_width) << it << " "; cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_traj(os); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->write_traj(os); } os << "\n"; if (cvm::debug()) { os.flush(); } cvm::decrease_depth(); return os; } void cvm::log(std::string const &message) { size_t const d = depth(); if (d > 0) proxy->log((std::string(2*d, ' '))+message); else proxy->log(message); } void cvm::increase_depth() { (depth())++; } void cvm::decrease_depth() { if (depth() > 0) { (depth())--; } } size_t & cvm::depth() { // NOTE: do not call log() or error() here, to avoid recursion size_t const nt = proxy->smp_num_threads(); if (proxy->smp_enabled() == COLVARS_OK) { if (depth_v.size() != nt) { // update array of depths proxy->smp_lock(); if (depth_v.size() > 0) { depth_s = depth_v[0]; } depth_v.clear(); depth_v.assign(nt, depth_s); proxy->smp_unlock(); } return depth_v[proxy->smp_thread_id()]; } return depth_s; } void colvarmodule::set_error_bits(int code) { if (code < 0) { cvm::fatal_error("Error: set_error_bits() received negative error code.\n"); return; } proxy->smp_lock(); errorCode |= code | COLVARS_ERROR; proxy->smp_unlock(); } bool colvarmodule::get_error_bit(int code) { return bool(errorCode & code); } void colvarmodule::clear_error() { proxy->smp_lock(); errorCode = COLVARS_OK; proxy->smp_unlock(); } void cvm::error(std::string const &message, int code) { set_error_bits(code); proxy->error(message); } void cvm::fatal_error(std::string const &message) { // TODO once all non-fatal errors have been set to be handled by error(), // set DELETE_COLVARS here for VMD to handle it set_error_bits(FATAL_ERROR); proxy->fatal_error(message); } void cvm::exit(std::string const &message) { proxy->exit(message); } int cvm::read_index_file(char const *filename) { std::ifstream is(filename, std::ios::binary); if (!is.good()) { cvm::error("Error: in opening index file \""+ std::string(filename)+"\".\n", FILE_ERROR); } while (is.good()) { char open, close; std::string group_name; if ( (is >> open) && (open == '[') && (is >> group_name) && (is >> close) && (close == ']') ) { for (std::list::iterator names_i = index_group_names.begin(); names_i != index_group_names.end(); names_i++) { if (*names_i == group_name) { cvm::error("Error: the group name \""+group_name+ "\" appears in multiple index files.\n", FILE_ERROR); } } cvm::index_group_names.push_back(group_name); cvm::index_groups.push_back(std::vector ()); } else { cvm::error("Error: in parsing index file \""+ std::string(filename)+"\".\n", INPUT_ERROR); } int atom_number = 1; size_t pos = is.tellg(); while ( (is >> atom_number) && (atom_number > 0) ) { (cvm::index_groups.back()).push_back(atom_number); pos = is.tellg(); } is.clear(); is.seekg(pos, std::ios::beg); std::string delim; if ( (is >> delim) && (delim == "[") ) { // new group is.clear(); is.seekg(pos, std::ios::beg); } else { break; } } cvm::log("The following index groups were read from the index file \""+ std::string(filename)+"\":\n"); std::list::iterator names_i = index_group_names.begin(); std::list >::iterator lists_i = index_groups.begin(); for ( ; names_i != index_group_names.end() ; names_i++, lists_i++) { cvm::log(" "+(*names_i)+" ("+cvm::to_str(lists_i->size())+" atoms).\n"); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int cvm::load_atoms(char const *file_name, cvm::atom_group &atoms, std::string const &pdb_field, double const pdb_field_value) { return proxy->load_atoms(file_name, atoms, pdb_field, pdb_field_value); } int cvm::load_coords(char const *file_name, std::vector &pos, const std::vector &indices, std::string const &pdb_field, double const pdb_field_value) { // Differentiate between PDB and XYZ files // for XYZ files, use CVM internal parser // otherwise call proxy function for PDB std::string const ext(strlen(file_name) > 4 ? (file_name + (strlen(file_name) - 4)) : file_name); if (colvarparse::to_lower_cppstr(ext) == std::string(".xyz")) { if ( pdb_field.size() > 0 ) { cvm::error("Error: PDB column may not be specified for XYZ coordinate file.\n", INPUT_ERROR); return COLVARS_ERROR; } return cvm::load_coords_xyz(file_name, pos, indices); } else { return proxy->load_coords(file_name, pos, indices, pdb_field, pdb_field_value); } } int cvm::load_coords_xyz(char const *filename, std::vector &pos, const std::vector &indices) { std::ifstream xyz_is(filename); unsigned int natoms; char symbol[256]; std::string line; if ( ! (xyz_is >> natoms) ) { cvm::error("Error: cannot parse XYZ file " + std::string(filename) + ".\n", INPUT_ERROR); } // skip comment line std::getline(xyz_is, line); std::getline(xyz_is, line); xyz_is.width(255); std::vector::iterator pos_i = pos.begin(); if (pos.size() != natoms) { // Use specified indices int next = 0; // indices are zero-based std::vector::const_iterator index = indices.begin(); for ( ; pos_i != pos.end() ; pos_i++, index++) { while ( next < *index ) { std::getline(xyz_is, line); next++; } xyz_is >> symbol; xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2]; } } else { // Use all positions for ( ; pos_i != pos.end() ; pos_i++) { xyz_is >> symbol; xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2]; } } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } // static pointers std::vector colvarmodule::colvars; std::vector colvarmodule::biases; size_t colvarmodule::n_abf_biases = 0; size_t colvarmodule::n_rest_biases = 0; size_t colvarmodule::n_histo_biases = 0; size_t colvarmodule::n_meta_biases = 0; colvarproxy *colvarmodule::proxy = NULL; // static runtime data cvm::real colvarmodule::debug_gradients_step_size = 1.0e-07; int colvarmodule::errorCode = 0; long colvarmodule::it = 0; long colvarmodule::it_restart = 0; size_t colvarmodule::restart_out_freq = 0; size_t colvarmodule::cv_traj_freq = 0; size_t colvarmodule::depth_s = 0; std::vector colvarmodule::depth_v(0); bool colvarmodule::b_analysis = false; std::list colvarmodule::index_group_names; std::list > colvarmodule::index_groups; bool colvarmodule::use_scripted_forces = false; bool colvarmodule::scripting_after_biases = true; // file name prefixes std::string colvarmodule::output_prefix = ""; std::string colvarmodule::restart_in_name = ""; // i/o constants size_t const colvarmodule::it_width = 12; size_t const colvarmodule::cv_prec = 14; size_t const colvarmodule::cv_width = 21; size_t const colvarmodule::en_prec = 14; size_t const colvarmodule::en_width = 21; std::string const colvarmodule::line_marker = "----------------------------------------------------------------------\n"; diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index b33d3ef1a..454fb65b6 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -1,671 +1,671 @@ // -*- c++ -*- #ifndef COLVARMODULE_H #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2016-09-30" +#define COLVARS_VERSION "2016-10-05" #endif #ifndef COLVARS_DEBUG #define COLVARS_DEBUG false #endif /*! \mainpage Main page */ /// \file colvarmodule.h /// \brief Collective variables main module /// /// This file declares the main class for defining and manipulating /// collective variables: there should be only one instance of this /// class, because several variables are made static (i.e. they are /// shared between all object instances) to be accessed from other /// objects. #define COLVARS_OK 0 #define COLVARS_ERROR 1 #define COLVARS_NOT_IMPLEMENTED (1<<1) #define INPUT_ERROR (1<<2) // out of bounds or inconsistent input #define BUG_ERROR (1<<3) // Inconsistent state indicating bug #define FILE_ERROR (1<<4) #define MEMORY_ERROR (1<<5) #define FATAL_ERROR (1<<6) // Should be set, or not, together with other bits #define DELETE_COLVARS (1<<7) // Instruct the caller to delete cvm #define COLVARS_NO_SUCH_FRAME (1<<8) // Cannot load the requested frame #include #include #include #include #include #include #include #include #include #ifdef NAMD_VERSION // use Lustre-friendly wrapper to POSIX write() #include "fstream_namd.h" #endif class colvarparse; class colvar; class colvarbias; class colvarproxy; class colvarscript; /// \brief Collective variables module (main class) /// /// Class to control the collective variables calculation. An object /// (usually one) of this class is spawned from the MD program, /// containing all i/o routines and general interface. /// /// At initialization, the colvarmodule object creates a proxy object /// to provide a transparent interface between the MD program and the /// child objects class colvarmodule { private: /// Impossible to initialize the main object without arguments colvarmodule(); public: friend class colvarproxy; // TODO colvarscript should be unaware of colvarmodule's internals friend class colvarscript; /// Defining an abstract real number allows to switch precision typedef double real; /// Residue identifier typedef int residue_id; class rvector; template class vector1d; template class matrix2d; class quaternion; class rotation; /// \brief Atom position (different type name from rvector, to make /// possible future PBC-transparent implementations) typedef rvector atom_pos; /// \brief 3x3 matrix of real numbers class rmatrix; // allow these classes to access protected data class atom; class atom_group; friend class atom; friend class atom_group; typedef std::vector::iterator atom_iter; typedef std::vector::const_iterator atom_const_iter; /// Module-wide error state /// see constants at the top of this file protected: static int errorCode; public: static void set_error_bits(int code); static bool get_error_bit(int code); static inline int get_error() { return errorCode; } static void clear_error(); /// Current step number static long it; /// Starting step number for this run static long it_restart; /// Return the current step number from the beginning of this run static inline long step_relative() { return it - it_restart; } /// Return the current step number from the beginning of the whole /// calculation static inline long step_absolute() { return it; } /// If true, get it_restart from the state file; if set to false, /// the MD program is providing it bool it_restart_from_state_file; /// \brief Finite difference step size (if there is no dynamics, or /// if gradients need to be tested independently from the size of /// dt) static real debug_gradients_step_size; /// Prefix for all output files for this run static std::string output_prefix; /// input restart file name (determined from input_prefix) static std::string restart_in_name; /// Array of collective variables static std::vector colvars; /* TODO: implement named CVCs /// Array of named (reusable) collective variable components static std::vector cvcs; /// Named cvcs register themselves at initialization time inline void register_cvc(cvc *p) { cvcs.push_back(p); } */ /// Collective variables to be calculated on different threads; /// colvars with multple items (e.g. multiple active CVCs) are duplicated std::vector colvars_smp; /// Indexes of the items to calculate for each colvar std::vector colvars_smp_items; /// Array of collective variable biases static std::vector biases; /// \brief Number of ABF biases initialized (in normal conditions /// should be 1) static size_t n_abf_biases; /// \brief Number of metadynamics biases initialized (in normal /// conditions should be 1) static size_t n_meta_biases; /// \brief Number of restraint biases initialized (no limit on the /// number) static size_t n_rest_biases; /// \brief Number of histograms initialized (no limit on the /// number) static size_t n_histo_biases; /// \brief Whether debug output should be enabled (compile-time option) static inline bool debug() { return COLVARS_DEBUG; } /// \brief Constructor \param config_name Configuration file name /// \param restart_name (optional) Restart file name colvarmodule(colvarproxy *proxy); /// Destructor ~colvarmodule(); /// Actual function called by the destructor int reset(); /// Open a config file, load its contents, and pass it to config_string() int read_config_file(char const *config_file_name); /// \brief Parse a config string assuming it is a complete configuration /// (i.e. calling all parse functions) int read_config_string(std::string const &conf); /// \brief Parse a "clean" config string (no comments) int parse_config(std::string &conf); // Parse functions (setup internal data based on a string) /// Parse the few module's global parameters int parse_global_params(std::string const &conf); /// Parse and initialize collective variables int parse_colvars(std::string const &conf); /// Parse and initialize collective variable biases int parse_biases(std::string const &conf); /// Parse and initialize collective variable biases of a specific type template int parse_biases_type(std::string const &conf, char const *keyword, size_t &bias_count); /// Test error condition and keyword parsing /// on error, delete new bias bool check_new_bias(std::string &conf, char const *key); private: /// Useful wrapper to interrupt parsing if any error occurs int catch_input_errors(int result); public: // "Setup" functions (change internal data based on related data // from the proxy that may change during program execution) // No additional parsing is done within these functions /// (Re)initialize internal data (currently used by LAMMPS) /// Also calls setup() member functions of colvars and biases int setup(); /// (Re)initialize and (re)read the input state file calling read_restart() int setup_input(); /// (Re)initialize the output trajectory and state file (does not write it yet) int setup_output(); #ifdef NAMD_VERSION typedef ofstream_namd ofstream; #else typedef std::ofstream ofstream; #endif /// Read the input restart file std::istream & read_restart(std::istream &is); /// Write the output restart file std::ostream & write_restart(std::ostream &os); /// Open a trajectory file if requested (and leave it open) int open_traj_file(std::string const &file_name); /// Close it int close_traj_file(); /// Write in the trajectory file std::ostream & write_traj(std::ostream &os); /// Write explanatory labels in the trajectory file std::ostream & write_traj_label(std::ostream &os); /// Write all trajectory files int write_traj_files(); /// Write all restart files int write_restart_files(); /// Write all FINAL output files int write_output_files(); /// Backup a file before writing it static int backup_file(char const *filename); /// Look up a bias by name; returns NULL if not found static colvarbias * bias_by_name(std::string const &name); /// Look up a colvar by name; returns NULL if not found static colvar * colvar_by_name(std::string const &name); /// Load new configuration for the given bias - /// currently works for harmonic (force constant and/or centers) int change_configuration(std::string const &bias_name, std::string const &conf); /// Read a colvar value std::string read_colvar(std::string const &name); /// Calculate change in energy from using alt. config. for the given bias - /// currently works for harmonic (force constant and/or centers) real energy_difference(std::string const &bias_name, std::string const &conf); /// Give the total number of bins for a given bias. int bias_bin_num(std::string const &bias_name); /// Calculate the bin index for a given bias. int bias_current_bin(std::string const &bias_name); //// Give the count at a given bin index. int bias_bin_count(std::string const &bias_name, size_t bin_index); //// Share among replicas. int bias_share(std::string const &bias_name); /// Main worker function int calc(); /// Calculate collective variables int calc_colvars(); /// Calculate biases int calc_biases(); /// Integrate bias and restraint forces, send colvar forces to atoms int update_colvar_forces(); /// Perform analysis int analyze(); /// \brief Read a collective variable trajectory (post-processing /// only, not called at runtime) int read_traj(char const *traj_filename, long traj_read_begin, long traj_read_end); /// Quick conversion of an object to a string template static std::string to_str(T const &x, size_t const &width = 0, size_t const &prec = 0); /// Quick conversion of a vector of objects to a string template static std::string to_str(std::vector const &x, size_t const &width = 0, size_t const &prec = 0); /// Reduce the number of characters in a string static inline std::string wrap_string(std::string const &s, size_t const &nchars) { if (!s.size()) return std::string(nchars, ' '); else return ( (s.size() <= size_t(nchars)) ? (s+std::string(nchars-s.size(), ' ')) : (std::string(s, 0, nchars)) ); } /// Number of characters to represent a time step static size_t const it_width; /// Number of digits to represent a collective variables value(s) static size_t const cv_prec; /// Number of characters to represent a collective variables value(s) static size_t const cv_width; /// Number of digits to represent the collective variables energy static size_t const en_prec; /// Number of characters to represent the collective variables energy static size_t const en_width; /// Line separator in the log output static std::string const line_marker; // proxy functions /// \brief Value of the unit for atomic coordinates with respect to /// angstroms (used by some variables for hard-coded default values) static real unit_angstrom(); /// \brief Boltmann constant static real boltzmann(); /// \brief Temperature of the simulation (K) static real temperature(); /// \brief Time step of MD integrator (fs) static real dt(); /// Request calculation of total force from MD engine static void request_total_force(); /// Print a message to the main log static void log(std::string const &message); /// Print a message to the main log and exit with error code static void fatal_error(std::string const &message); /// Print a message to the main log and set global error code static void error(std::string const &message, int code = COLVARS_ERROR); /// Print a message to the main log and exit normally static void exit(std::string const &message); // Replica exchange commands. static bool replica_enabled(); static int replica_index(); static int replica_num(); static void replica_comm_barrier(); static int replica_comm_recv(char* msg_data, int buf_len, int src_rep); static int replica_comm_send(char* msg_data, int msg_len, int dest_rep); /// \brief Get the distance between two atomic positions with pbcs handled /// correctly static rvector position_distance(atom_pos const &pos1, atom_pos const &pos2); /// \brief Get the square distance between two positions (with /// periodic boundary conditions handled transparently) /// /// Note: in the case of periodic boundary conditions, this provides /// an analytical square distance (while taking the square of /// position_distance() would produce leads to a cusp) static real position_dist2(atom_pos const &pos1, atom_pos const &pos2); /// \brief Get the closest periodic image to a reference position /// \param pos The position to look for the closest periodic image /// \param ref_pos (optional) The reference position static void select_closest_image(atom_pos &pos, atom_pos const &ref_pos); /// \brief Perform select_closest_image() on a set of atomic positions /// /// After that, distance vectors can then be calculated directly, /// without using position_distance() static void select_closest_images(std::vector &pos, atom_pos const &ref_pos); /// \brief Names of groups from a Gromacs .ndx file to be read at startup static std::list index_group_names; /// \brief Groups from a Gromacs .ndx file read at startup static std::list > index_groups; /// \brief Read a Gromacs .ndx file static int read_index_file(char const *filename); /// \brief Create atoms from a file \param filename name of the file /// (usually a PDB) \param atoms array of the atoms to be allocated /// \param pdb_field (optiona) if "filename" is a PDB file, use this /// field to determine which are the atoms to be set static int load_atoms(char const *filename, atom_group &atoms, std::string const &pdb_field, double const pdb_field_value = 0.0); /// \brief Load the coordinates for a group of atoms from a file /// (PDB or XYZ) static int load_coords(char const *filename, std::vector &pos, const std::vector &indices, std::string const &pdb_field, double const pdb_field_value = 0.0); /// \brief Load the coordinates for a group of atoms from an /// XYZ file static int load_coords_xyz(char const *filename, std::vector &pos, const std::vector &indices); /// Frequency for collective variables trajectory output static size_t cv_traj_freq; /// \brief True if only analysis is performed and not a run static bool b_analysis; /// Frequency for saving output restarts static size_t restart_out_freq; /// Output restart file name std::string restart_out_name; /// Pseudo-random number with Gaussian distribution static real rand_gaussian(void); protected: /// Configuration file std::ifstream config_s; /// Configuration file parser object colvarparse *parse; /// Name of the trajectory file std::string cv_traj_name; /// Collective variables output trajectory file colvarmodule::ofstream cv_traj_os; /// Appending to the existing trajectory file? bool cv_traj_append; /// Output restart file colvarmodule::ofstream restart_out_os; protected: /// Counter for the current depth in the object hierarchy (useg e.g. in output) static size_t depth_s; /// Thread-specific depth static std::vector depth_v; public: /// Get the current object depth in the hierarchy static size_t & depth(); /// Increase the depth (number of indentations in the output) static void increase_depth(); /// Decrease the depth (number of indentations in the output) static void decrease_depth(); static inline bool scripted_forces() { return use_scripted_forces; } /// Use scripted colvars forces? static bool use_scripted_forces; /// Wait for all biases before calculating scripted forces? static bool scripting_after_biases; /// Calculate the energy and forces of scripted biases int calc_scripted_forces(); /// \brief Pointer to the proxy object, used to retrieve atomic data /// from the hosting program; it is static in order to be accessible /// from static functions in the colvarmodule class static colvarproxy *proxy; }; /// Shorthand for the frequently used type prefix typedef colvarmodule cvm; #include "colvartypes.h" std::ostream & operator << (std::ostream &os, cvm::rvector const &v); std::istream & operator >> (std::istream &is, cvm::rvector &v); template std::string cvm::to_str(T const &x, size_t const &width, size_t const &prec) { std::ostringstream os; if (width) os.width(width); if (prec) { os.setf(std::ios::scientific, std::ios::floatfield); os.precision(prec); } os << x; return os.str(); } template std::string cvm::to_str(std::vector const &x, size_t const &width, size_t const &prec) { if (!x.size()) return std::string(""); std::ostringstream os; if (prec) { os.setf(std::ios::scientific, std::ios::floatfield); } os << "{ "; if (width) os.width(width); if (prec) os.precision(prec); os << x[0]; for (size_t i = 1; i < x.size(); i++) { os << ", "; if (width) os.width(width); if (prec) os.precision(prec); os << x[i]; } os << " }"; return os.str(); } #include "colvarproxy.h" inline cvm::real cvm::unit_angstrom() { return proxy->unit_angstrom(); } inline cvm::real cvm::boltzmann() { return proxy->boltzmann(); } inline cvm::real cvm::temperature() { return proxy->temperature(); } inline cvm::real cvm::dt() { return proxy->dt(); } // Replica exchange commands inline bool cvm::replica_enabled() { return proxy->replica_enabled(); } inline int cvm::replica_index() { return proxy->replica_index(); } inline int cvm::replica_num() { return proxy->replica_num(); } inline void cvm::replica_comm_barrier() { return proxy->replica_comm_barrier(); } inline int cvm::replica_comm_recv(char* msg_data, int buf_len, int src_rep) { return proxy->replica_comm_recv(msg_data,buf_len,src_rep); } inline int cvm::replica_comm_send(char* msg_data, int msg_len, int dest_rep) { return proxy->replica_comm_send(msg_data,msg_len,dest_rep); } inline void cvm::request_total_force() { proxy->request_total_force(true); } inline void cvm::select_closest_image(atom_pos &pos, atom_pos const &ref_pos) { proxy->select_closest_image(pos, ref_pos); } inline void cvm::select_closest_images(std::vector &pos, atom_pos const &ref_pos) { proxy->select_closest_images(pos, ref_pos); } inline cvm::rvector cvm::position_distance(atom_pos const &pos1, atom_pos const &pos2) { return proxy->position_distance(pos1, pos2); } inline cvm::real cvm::position_dist2(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) { return proxy->position_dist2(pos1, pos2); } inline cvm::real cvm::rand_gaussian(void) { return proxy->rand_gaussian(); } #endif diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index bd5ee3da4..3d41cb58a 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -1,611 +1,619 @@ // -*- c++ -*- #ifndef COLVARPROXY_H #define COLVARPROXY_H #include #include #include "colvarmodule.h" #include "colvarvalue.h" // forward declarations class colvarscript; /// \brief Interface between the collective variables module and /// the simulation or analysis program (NAMD, VMD, LAMMPS...). /// This is the base class: each interfaced program is supported by a derived class. /// Only pure virtual functions ("= 0") must be reimplemented to ensure baseline functionality. class colvarproxy { public: /// Pointer to the main object colvarmodule *colvars; /// Default constructor inline colvarproxy() : script(NULL) {} /// Default destructor virtual ~colvarproxy() {} /// (Re)initialize required member data after construction virtual int setup() { return COLVARS_OK; } /// \brief Update data required by the colvars module (e.g. cache atom positions) /// /// TODO Break up colvarproxy_namd and colvarproxy_lammps function into these virtual int update_input() { return COLVARS_OK; } /// \brief Update data based from the results of a module update (e.g. send forces) virtual int update_output() { return COLVARS_OK; } // **************** SIMULATION PARAMETERS **************** /// \brief Value of the unit for atomic coordinates with respect to /// angstroms (used by some variables for hard-coded default values) virtual cvm::real unit_angstrom() = 0; /// \brief Boltzmann constant virtual cvm::real boltzmann() = 0; /// \brief Temperature of the simulation (K) virtual cvm::real temperature() = 0; /// \brief Time step of the simulation (fs) virtual cvm::real dt() = 0; /// \brief Pseudo-random number with Gaussian distribution virtual cvm::real rand_gaussian(void) = 0; /// \brief Get the current frame number // Returns error code virtual int get_frame(long int&) { return COLVARS_NOT_IMPLEMENTED; } /// \brief Set the current frame number (as well as colvarmodule::it) // Returns error code virtual int set_frame(long int) { return COLVARS_NOT_IMPLEMENTED; } /// \brief Prefix to be used for input files (restarts, not /// configuration) std::string input_prefix_str, output_prefix_str, restart_output_prefix_str; inline std::string input_prefix() { return input_prefix_str; } /// \brief Prefix to be used for output restart files inline std::string restart_output_prefix() { return restart_output_prefix_str; } /// \brief Prefix to be used for output files (final system /// configuration) inline std::string output_prefix() { return output_prefix_str; } /// \brief Restarts will be written each time this number of steps has passed virtual size_t restart_frequency() { return 0; } protected: /// \brief Currently opened output files: by default, these are ofstream objects. /// Allows redefinition to implement different output mechanisms std::list output_files; /// \brief Identifiers for output_stream objects: by default, these are the names of the files std::list output_stream_names; public: // ***************** SHARED-MEMORY PARALLELIZATION ***************** /// Whether or not threaded parallelization is available virtual int smp_enabled() { return COLVARS_NOT_IMPLEMENTED; } /// Distribute calculation of colvars (and their components) across threads virtual int smp_colvars_loop() { return COLVARS_NOT_IMPLEMENTED; } /// Distribute calculation of biases across threads virtual int smp_biases_loop() { return COLVARS_NOT_IMPLEMENTED; } /// Distribute calculation of biases across threads 2nd through last, with all scripted biased on 1st thread virtual int smp_biases_script_loop() { return COLVARS_NOT_IMPLEMENTED; } /// Index of this thread virtual int smp_thread_id() { return COLVARS_NOT_IMPLEMENTED; } /// Number of threads sharing this address space virtual int smp_num_threads() { return COLVARS_NOT_IMPLEMENTED; } /// Lock the proxy's shared data for access by a thread, if threads are implemented; if not implemented, does nothing virtual int smp_lock() { return COLVARS_OK; } /// Attempt to lock the proxy's shared data virtual int smp_trylock() { return COLVARS_OK; } /// Release the lock virtual int smp_unlock() { return COLVARS_OK; } // **************** MULTIPLE REPLICAS COMMUNICATION **************** // Replica exchange commands: /// \brief Indicate if multi-replica support is available and active virtual bool replica_enabled() { return false; } /// \brief Index of this replica virtual int replica_index() { return 0; } /// \brief Total number of replica virtual int replica_num() { return 1; } /// \brief Synchronize replica virtual void replica_comm_barrier() {} /// \brief Receive data from other replica virtual int replica_comm_recv(char* msg_data, int buf_len, int src_rep) { return COLVARS_NOT_IMPLEMENTED; } /// \brief Send data to other replica virtual int replica_comm_send(char* msg_data, int msg_len, int dest_rep) { return COLVARS_NOT_IMPLEMENTED; } // **************** SCRIPTING INTERFACE **************** /// Pointer to the scripting interface object /// (does not need to be allocated in a new interface) colvarscript *script; /// is a user force script defined? bool force_script_defined; /// Do we have a scripting interface? bool have_scripts; /// Run a user-defined colvar forces script virtual int run_force_callback() { return COLVARS_NOT_IMPLEMENTED; } virtual int run_colvar_callback(std::string const &name, std::vector const &cvcs, colvarvalue &value) { return COLVARS_NOT_IMPLEMENTED; } virtual int run_colvar_gradient_callback(std::string const &name, std::vector const &cvcs, std::vector > &gradient) { return COLVARS_NOT_IMPLEMENTED; } // **************** INPUT/OUTPUT **************** /// Print a message to the main log virtual void log(std::string const &message) = 0; /// Print a message to the main log and let the rest of the program handle the error virtual void error(std::string const &message) = 0; /// Print a message to the main log and exit with error code virtual void fatal_error(std::string const &message) = 0; /// Print a message to the main log and exit normally virtual void exit(std::string const &message) { cvm::error("Error: exiting without error is not implemented, returning error code.\n", COLVARS_NOT_IMPLEMENTED); } // TODO the following definitions may be moved to a .cpp file /// \brief Returns a reference to the given output channel; /// if this is not open already, then open it virtual std::ostream * output_stream(std::string const &output_name) { std::list::iterator osi = output_files.begin(); std::list::iterator osni = output_stream_names.begin(); for ( ; osi != output_files.end(); osi++, osni++) { if (*osni == output_name) { return *osi; } } output_stream_names.push_back(output_name); std::ofstream * os = new std::ofstream(output_name.c_str()); if (!os->is_open()) { cvm::error("Error: cannot write to file \""+output_name+"\".\n", FILE_ERROR); } output_files.push_back(os); return os; } /// \brief Closes the given output channel virtual int close_output_stream(std::string const &output_name) { std::list::iterator osi = output_files.begin(); std::list::iterator osni = output_stream_names.begin(); for ( ; osi != output_files.end(); osi++, osni++) { if (*osni == output_name) { ((std::ofstream *) (*osi))->close(); output_files.erase(osi); output_stream_names.erase(osni); return COLVARS_OK; } } cvm::error("Error: trying to close an output file or stream that wasn't open.\n", BUG_ERROR); return COLVARS_ERROR; } /// \brief Rename the given file, before overwriting it virtual int backup_file(char const *filename) { return COLVARS_NOT_IMPLEMENTED; } // **************** ACCESS SYSTEM DATA **************** /// Pass restraint energy value for current timestep to MD engine virtual void add_energy(cvm::real energy) = 0; /// Tell the proxy whether total forces are needed (may not always be available) virtual void request_total_force(bool yesno) { if (yesno == true) cvm::error("Error: total forces are currently not implemented.\n", COLVARS_NOT_IMPLEMENTED); } + /// Are total forces being used? + virtual bool total_forces_enabled() const + { + cvm::error("Error: total forces are currently not implemented.\n", + COLVARS_NOT_IMPLEMENTED); + return false; + } + /// \brief Get the PBC-aware distance vector between two positions virtual cvm::rvector position_distance(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) = 0; /// \brief Get the PBC-aware square distance between two positions; /// may need to be reimplemented independently from position_distance() for optimization purposes virtual cvm::real position_dist2(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) { return (position_distance(pos1, pos2)).norm2(); } /// \brief Get the closest periodic image to a reference position /// \param pos The position to look for the closest periodic image /// \param ref_pos The reference position virtual void select_closest_image(cvm::atom_pos &pos, cvm::atom_pos const &ref_pos) { pos = position_distance(ref_pos, pos) + ref_pos; } /// \brief Perform select_closest_image() on a set of atomic positions /// /// After that, distance vectors can then be calculated directly, /// without using position_distance() void select_closest_images(std::vector &pos, cvm::atom_pos const &ref_pos) { for (std::vector::iterator pi = pos.begin(); pi != pos.end(); ++pi) { select_closest_image(*pi, ref_pos); } } // **************** ACCESS ATOMIC DATA **************** protected: /// \brief Array of 0-based integers used to uniquely associate atoms /// within the host program std::vector atoms_ids; /// \brief Keep track of how many times each atom is used by a separate colvar object std::vector atoms_ncopies; /// \brief Masses of the atoms (allow redefinition during a run, as done e.g. in LAMMPS) std::vector atoms_masses; /// \brief Charges of the atoms (allow redefinition during a run, as done e.g. in LAMMPS) std::vector atoms_charges; /// \brief Current three-dimensional positions of the atoms std::vector atoms_positions; /// \brief Most recent total forces on each atom std::vector atoms_total_forces; /// \brief Forces applied from colvars, to be communicated to the MD integrator std::vector atoms_new_colvar_forces; /// Used by all init_atom() functions: create a slot for an atom not requested yet inline int add_atom_slot(int atom_id) { atoms_ids.push_back(atom_id); atoms_ncopies.push_back(1); atoms_masses.push_back(1.0); atoms_charges.push_back(0.0); atoms_positions.push_back(cvm::rvector(0.0, 0.0, 0.0)); atoms_total_forces.push_back(cvm::rvector(0.0, 0.0, 0.0)); atoms_new_colvar_forces.push_back(cvm::rvector(0.0, 0.0, 0.0)); return (atoms_ids.size() - 1); } public: /// Prepare this atom for collective variables calculation, selecting it by numeric index (1-based) virtual int init_atom(int atom_number) = 0; /// Check that this atom number is valid, but do not initialize the corresponding atom yet virtual int check_atom_id(int atom_number) = 0; /// Select this atom for collective variables calculation, using name and residue number. /// Not all programs support this: leave this function as is in those cases. virtual int init_atom(cvm::residue_id const &residue, std::string const &atom_name, std::string const &segment_id) { cvm::error("Error: initializing an atom by name and residue number is currently not supported.\n", COLVARS_NOT_IMPLEMENTED); return COLVARS_NOT_IMPLEMENTED; } /// Check that this atom is valid, but do not initialize it yet virtual int check_atom_id(cvm::residue_id const &residue, std::string const &atom_name, std::string const &segment_id) { cvm::error("Error: initializing an atom by name and residue number is currently not supported.\n", COLVARS_NOT_IMPLEMENTED); return COLVARS_NOT_IMPLEMENTED; } /// \brief Used by the atom class destructor: rather than deleting the array slot /// (costly) set the corresponding atoms_ncopies to zero virtual void clear_atom(int index) { if (((size_t) index) >= atoms_ids.size()) { cvm::error("Error: trying to disable an atom that was not previously requested.\n", INPUT_ERROR); } if (atoms_ncopies[index] > 0) { atoms_ncopies[index] -= 1; } } /// Get the numeric ID of the given atom (for the program) inline int get_atom_id(int index) const { return atoms_ids[index]; } /// Get the mass of the given atom inline cvm::real get_atom_mass(int index) const { return atoms_masses[index]; } /// Get the charge of the given atom inline cvm::real get_atom_charge(int index) const { return atoms_charges[index]; } /// Read the current position of the given atom inline cvm::rvector get_atom_position(int index) const { return atoms_positions[index]; } /// Read the current total force of the given atom inline cvm::rvector get_atom_total_force(int index) const { return atoms_total_forces[index]; } /// Request that this force is applied to the given atom inline void apply_atom_force(int index, cvm::rvector const &new_force) { atoms_new_colvar_forces[index] += new_force; } /// Read the current velocity of the given atom virtual cvm::rvector get_atom_velocity(int index) { cvm::error("Error: reading the current velocity of an atom is not yet implemented.\n", COLVARS_NOT_IMPLEMENTED); return cvm::rvector(0.0); } // useful functions for data management outside this class inline std::vector *modify_atom_ids() { return &atoms_ids; } inline std::vector *modify_atom_masses() { return &atoms_masses; } inline std::vector *modify_atom_charges() { return &atoms_charges; } inline std::vector *modify_atom_positions() { return &atoms_positions; } inline std::vector *modify_atom_total_forces() { return &atoms_total_forces; } inline std::vector *modify_atom_new_colvar_forces() { return &atoms_new_colvar_forces; } /// \brief Read atom identifiers from a file \param filename name of /// the file (usually a PDB) \param atoms array to which atoms read /// from "filename" will be appended \param pdb_field (optiona) if /// "filename" is a PDB file, use this field to determine which are /// the atoms to be set virtual int load_atoms(char const *filename, cvm::atom_group &atoms, std::string const &pdb_field, double const pdb_field_value = 0.0) { cvm::error("Error: loading atom identifiers from a file is currently not implemented.\n", COLVARS_NOT_IMPLEMENTED); return COLVARS_NOT_IMPLEMENTED; } /// \brief Load the coordinates for a group of atoms from a file /// (usually a PDB); if "pos" is already allocated, the number of its /// elements must match the number of atoms in "filename" virtual int load_coords(char const *filename, std::vector &pos, const std::vector &indices, std::string const &pdb_field, double const pdb_field_value = 0.0) { cvm::error("Error: loading atomic coordinates from a file is currently not implemented.\n"); return COLVARS_NOT_IMPLEMENTED; } // **************** ACCESS GROUP DATA **************** protected: /// \brief Array of 0-based integers used to uniquely associate atom groups /// within the host program std::vector atom_groups_ids; /// \brief Keep track of how many times each group is used by a separate cvc std::vector atom_groups_ncopies; /// \brief Total masses of the atom groups std::vector atom_groups_masses; /// \brief Total charges of the atom groups (allow redefinition during a run, as done e.g. in LAMMPS) std::vector atom_groups_charges; /// \brief Current centers of mass of the atom groups std::vector atom_groups_coms; /// \brief Most recently updated total forces on the com of each group std::vector atom_groups_total_forces; /// \brief Forces applied from colvars, to be communicated to the MD integrator std::vector atom_groups_new_colvar_forces; /// TODO Add here containers of handles to cvc objects that are computed in parallel public: /// \brief Whether this proxy implementation has capability for scalable groups virtual int scalable_group_coms() { return COLVARS_NOT_IMPLEMENTED; } /// Used by all init_atom_group() functions: create a slot for an atom group not requested yet // TODO Add a handle to cvc objects inline int add_atom_group_slot(int atom_group_id) { atom_groups_ids.push_back(atom_group_id); atom_groups_ncopies.push_back(1); atom_groups_masses.push_back(1.0); atom_groups_charges.push_back(0.0); atom_groups_coms.push_back(cvm::rvector(0.0, 0.0, 0.0)); atom_groups_total_forces.push_back(cvm::rvector(0.0, 0.0, 0.0)); atom_groups_new_colvar_forces.push_back(cvm::rvector(0.0, 0.0, 0.0)); return (atom_groups_ids.size() - 1); } /// Prepare this group for collective variables calculation, selecting atoms by internal ids (0-based) virtual int init_atom_group(std::vector const &atoms_ids) // TODO Add a handle to cvc objects { cvm::error("Error: initializing a group outside of the colvars module is currently not supported.\n", COLVARS_NOT_IMPLEMENTED); return COLVARS_NOT_IMPLEMENTED; } /// \brief Used by the atom_group class destructor virtual void clear_atom_group(int index) { if (cvm::debug()) { log("Trying to remove/disable atom group number "+cvm::to_str(index)+"\n"); } if (((size_t) index) >= atom_groups_ids.size()) { cvm::error("Error: trying to disable an atom group that was not previously requested.\n", INPUT_ERROR); } if (atom_groups_ncopies[index] > 0) { atom_groups_ncopies[index] -= 1; } } /// Get the numeric ID of the given atom group (for the MD program) inline cvm::real get_atom_group_id(int index) const { return atom_groups_ids[index]; } /// Get the mass of the given atom group inline cvm::real get_atom_group_mass(int index) const { return atom_groups_masses[index]; } /// Get the charge of the given atom group inline cvm::real get_atom_group_charge(int index) const { return atom_groups_charges[index]; } /// Read the current position of the center of mass given atom group inline cvm::rvector get_atom_group_com(int index) const { return atom_groups_coms[index]; } /// Read the current total force of the given atom group inline cvm::rvector get_atom_group_total_force(int index) const { return atom_groups_total_forces[index]; } /// Request that this force is applied to the given atom group inline void apply_atom_group_force(int index, cvm::rvector const &new_force) { atom_groups_new_colvar_forces[index] += new_force; } /// Read the current velocity of the given atom group virtual cvm::rvector get_atom_group_velocity(int index) { cvm::error("Error: reading the current velocity of an atom group is not yet implemented.\n", COLVARS_NOT_IMPLEMENTED); return cvm::rvector(0.0); } }; #endif diff --git a/src/USER-COLVARS/colvarproxy_lammps.cpp b/src/USER-COLVARS/colvarproxy_lammps.cpp index f02fd2395..5e63d07c3 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.cpp +++ b/src/USER-COLVARS/colvarproxy_lammps.cpp @@ -1,476 +1,476 @@ #include #include "lammps.h" #include "atom.h" #include "error.h" #include "output.h" #include "random_park.h" #include "fix_colvars.h" #include "colvarmodule.h" #include "colvar.h" #include "colvarbias.h" #include "colvaratoms.h" #include "colvarproxy.h" #include "colvarproxy_lammps.h" #include #include #include #include #include #include #include #include #include #define HASH_FAIL -1 //////////////////////////////////////////////////////////////////////// // local helper functions // safely move filename to filename.extension static int my_backup_file(const char *filename, const char *extension) { struct stat sbuf; if (stat(filename, &sbuf) == 0) { if (!extension) extension = ".BAK"; char *backup = new char[strlen(filename)+strlen(extension)+1]; strcpy(backup, filename); strcat(backup, extension); #if defined(_WIN32) && !defined(__CYGWIN__) remove(backup); #endif if (rename(filename,backup)) { char *sys_err_msg = strerror(errno); if (!sys_err_msg) sys_err_msg = (char *) "(unknown error)"; fprintf(stderr,"Error renaming file %s to %s: %s\n", filename, backup, sys_err_msg); delete [] backup; return COLVARS_ERROR; } delete [] backup; } return COLVARS_OK; } //////////////////////////////////////////////////////////////////////// colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp, const char *inp_name, const char *out_name, const int seed, const double temp, MPI_Comm root2root) : _lmp(lmp), inter_comm(root2root) { if (cvm::debug()) log("Initializing the colvars proxy object.\n"); _random = new LAMMPS_NS::RanPark(lmp,seed); first_timestep=true; total_force_requested=false; previous_step=-1; t_target=temp; do_exit=false; restart_every=0; // User-scripted forces are not available in LAMMPS force_script_defined = false; have_scripts = false; // set input restart name and strip the extension, if present input_prefix_str = std::string(inp_name ? inp_name : ""); if (input_prefix_str.rfind(".colvars.state") != std::string::npos) input_prefix_str.erase(input_prefix_str.rfind(".colvars.state"), std::string(".colvars.state").size()); // output prefix is always given output_prefix_str = std::string(out_name); // not so for restarts restart_output_prefix_str = std::string("rest"); // check if it is possible to save output configuration if ((!output_prefix_str.size()) && (!restart_output_prefix_str.size())) { fatal_error("Error: neither the final output state file or " "the output restart file could be defined, exiting.\n"); } // try to extract a restart prefix from a potential restart command. LAMMPS_NS::Output *outp = _lmp->output; if ((outp->restart_every_single > 0) && (outp->restart1 != 0)) { restart_output_prefix_str = std::string(outp->restart1); } else if ((outp->restart_every_double > 0) && (outp->restart2a != 0)) { restart_output_prefix_str = std::string(outp->restart2a); } // trim off unwanted stuff from the restart prefix if (restart_output_prefix_str.rfind(".*") != std::string::npos) restart_output_prefix_str.erase(restart_output_prefix_str.rfind(".*"),2); #if defined(_OPENMP) if (smp_thread_id() == 0) { omp_init_lock(&smp_lock_state); } #endif // initialize multi-replica support, if available if (replica_enabled()) { MPI_Comm_rank(inter_comm, &inter_me); MPI_Comm_size(inter_comm, &inter_num); } if (cvm::debug()) log("Done initializing the colvars proxy object.\n"); } void colvarproxy_lammps::init(const char *conf_file) { // create the colvarmodule instance colvars = new colvarmodule(this); cvm::log("Using LAMMPS interface, version "+ cvm::to_str(COLVARPROXY_VERSION)+".\n"); my_angstrom = _lmp->force->angstrom; my_boltzmann = _lmp->force->boltz; my_timestep = _lmp->update->dt * _lmp->force->femtosecond; // TODO move one or more of these to setup() if needed colvars->read_config_file(conf_file); colvars->setup_input(); colvars->setup_output(); if (_lmp->update->ntimestep != 0) { cvm::log("Initializing step number as firstTimestep.\n"); colvars->it = colvars->it_restart = _lmp->update->ntimestep; } if (cvm::debug()) { log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n"); log(cvm::line_marker); log("Info: done initializing the colvars proxy object.\n"); } } colvarproxy_lammps::~colvarproxy_lammps() { delete _random; if (colvars != NULL) { colvars->write_output_files(); delete colvars; colvars = NULL; } } // re-initialize data where needed int colvarproxy_lammps::setup() { my_timestep = _lmp->update->dt * _lmp->force->femtosecond; return colvars->setup(); } // trigger colvars computation double colvarproxy_lammps::compute() { if (first_timestep) { first_timestep = false; } else { // Use the time step number inherited from LAMMPS if ( _lmp->update->ntimestep - previous_step == 1 ) colvars->it++; // Other cases could mean: // - run 0 // - beginning of a new run statement // then the internal counter should not be incremented } previous_step = _lmp->update->ntimestep; if (cvm::debug()) { cvm::log(cvm::line_marker+ "colvarproxy_lammps, step no. "+cvm::to_str(colvars->it)+"\n"+ "Updating internal data.\n"); } // zero the forces on the atoms, so that they can be accumulated by the colvars for (size_t i = 0; i < atoms_new_colvar_forces.size(); i++) { atoms_new_colvar_forces[i].reset(); } bias_energy = 0.0; if (cvm::debug()) { log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n"); log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n"); } // call the collective variable module colvars->calc(); if (cvm::debug()) { log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n"); log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n"); } return bias_energy; } void colvarproxy_lammps::serialize_status(std::string &rst) { std::ostringstream os; colvars->write_restart(os); rst = os.str(); } // set status from string bool colvarproxy_lammps::deserialize_status(std::string &rst) { std::istringstream is; is.str(rst); if (!colvars->read_restart(is)) { return false; } else { return true; } } cvm::rvector colvarproxy_lammps::position_distance(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) { double xtmp = pos2.x - pos1.x; double ytmp = pos2.y - pos1.y; double ztmp = pos2.z - pos1.z; _lmp->domain->minimum_image(xtmp,ytmp,ztmp); return cvm::rvector(xtmp, ytmp, ztmp); } cvm::real colvarproxy_lammps::position_dist2(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) { double xtmp = pos2.x - pos1.x; double ytmp = pos2.y - pos1.y; double ztmp = pos2.z - pos1.z; _lmp->domain->minimum_image(xtmp,ytmp,ztmp); return cvm::real(xtmp*xtmp + ytmp*ytmp + ztmp*ztmp); } void colvarproxy_lammps::select_closest_image(cvm::atom_pos &pos, cvm::atom_pos const &ref) { double xtmp = pos.x - ref.x; double ytmp = pos.y - ref.y; double ztmp = pos.z - ref.z; _lmp->domain->minimum_image(xtmp,ytmp,ztmp); pos.x = ref.x + xtmp; pos.y = ref.y + ytmp; pos.z = ref.z + ztmp; } void colvarproxy_lammps::log(std::string const &message) { std::istringstream is(message); std::string line; while (std::getline(is, line)) { if (_lmp->screen) fprintf(_lmp->screen,"colvars: %s\n",line.c_str()); if (_lmp->logfile) fprintf(_lmp->logfile,"colvars: %s\n",line.c_str()); } } void colvarproxy_lammps::error(std::string const &message) { // In LAMMPS, all errors are fatal fatal_error(message); } void colvarproxy_lammps::fatal_error(std::string const &message) { log(message); - if (!cvm::debug()) - log("If this error message is unclear, try recompiling the " - "colvars library and LAMMPS with -DCOLVARS_DEBUG.\n"); + // if (!cvm::debug()) + // log("If this error message is unclear, try recompiling the " + // "colvars library and LAMMPS with -DCOLVARS_DEBUG.\n"); _lmp->error->one(FLERR, "Fatal error in the collective variables module.\n"); } void colvarproxy_lammps::exit(std::string const &message) { log(message); log("Request to exit the simulation made.\n"); do_exit=true; } int colvarproxy_lammps::backup_file(char const *filename) { if (std::string(filename).rfind(std::string(".colvars.state")) != std::string::npos) { return my_backup_file(filename, ".old"); } else { return my_backup_file(filename, ".BAK"); } } #if defined(_OPENMP) // SMP support int colvarproxy_lammps::smp_enabled() { return COLVARS_OK; } int colvarproxy_lammps::smp_colvars_loop() { colvarmodule *cv = this->colvars; #pragma omp parallel for for (size_t i = 0; i < cv->colvars_smp.size(); i++) { if (cvm::debug()) { cvm::log("Calculating colvar \""+cv->colvars_smp[i]->name+"\" on thread "+cvm::to_str(smp_thread_id())+"\n"); } cv->colvars_smp[i]->calc_cvcs(cv->colvars_smp_items[i], 1); } return cvm::get_error(); } int colvarproxy_lammps::smp_biases_loop() { colvarmodule *cv = this->colvars; #pragma omp parallel for for (size_t i = 0; i < cv->biases.size(); i++) { if (cvm::debug()) { cvm::log("Calculating bias \""+cv->biases[i]->name+"\" on thread "+cvm::to_str(smp_thread_id())+"\n"); } cv->biases[i]->update(); } return cvm::get_error(); } int colvarproxy_lammps::smp_thread_id() { return omp_get_thread_num(); } int colvarproxy_lammps::smp_num_threads() { return omp_get_max_threads(); } int colvarproxy_lammps::smp_lock() { omp_set_lock(&smp_lock_state); return COLVARS_OK; } int colvarproxy_lammps::smp_trylock() { return omp_test_lock(&smp_lock_state) ? COLVARS_OK : COLVARS_ERROR; } int colvarproxy_lammps::smp_unlock() { omp_unset_lock(&smp_lock_state); return COLVARS_OK; } #endif // multi-replica support void colvarproxy_lammps::replica_comm_barrier() { MPI_Barrier(inter_comm); } int colvarproxy_lammps::replica_comm_recv(char* msg_data, int buf_len, int src_rep) { MPI_Status status; int retval; retval = MPI_Recv(msg_data,buf_len,MPI_CHAR,src_rep,0,inter_comm,&status); if (retval == MPI_SUCCESS) { MPI_Get_count(&status, MPI_CHAR, &retval); } else retval = 0; return retval; } int colvarproxy_lammps::replica_comm_send(char* msg_data, int msg_len, int dest_rep) { int retval; retval = MPI_Send(msg_data,msg_len,MPI_CHAR,dest_rep,0,inter_comm); if (retval == MPI_SUCCESS) { retval = msg_len; } else retval = 0; return retval; } int colvarproxy_lammps::check_atom_id(int atom_number) { int const aid = atom_number; if (cvm::debug()) log("Adding atom "+cvm::to_str(atom_number)+ " for collective variables calculation.\n"); // TODO add upper boundary check? if ( (aid < 0) ) { cvm::error("Error: invalid atom number specified, "+ cvm::to_str(atom_number)+"\n", INPUT_ERROR); return INPUT_ERROR; } return aid; } int colvarproxy_lammps::init_atom(int atom_number) { int aid = atom_number; for (size_t i = 0; i < atoms_ids.size(); i++) { if (atoms_ids[i] == aid) { // this atom id was already recorded atoms_ncopies[i] += 1; return i; } } aid = check_atom_id(atom_number); if (aid < 0) { return aid; } int const index = colvarproxy::add_atom_slot(aid); // add entries for the LAMMPS-specific fields atoms_types.push_back(0); return index; } diff --git a/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h index 925fefac0..f1636e112 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.h +++ b/src/USER-COLVARS/colvarproxy_lammps.h @@ -1,174 +1,174 @@ // -*- c++ -*- #ifndef COLVARPROXY_LAMMPS_H #define COLVARPROXY_LAMMPS_H #include "colvarmodule.h" #include "colvarproxy.h" #include "colvarvalue.h" #include "lammps.h" #include "domain.h" #include "force.h" #include "update.h" #include #include #include #if defined(_OPENMP) #include #endif #ifndef COLVARPROXY_VERSION -#define COLVARPROXY_VERSION "2016-04-08" +#define COLVARPROXY_VERSION "2016-10-05" #endif /* struct for packed data communication of coordinates and forces. */ struct commdata { int tag,type; double x,y,z,m,q; }; inline std::ostream & operator<< (std::ostream &out, const commdata &cd) { out << " (" << cd.tag << "/" << cd.type << ": " << cd.x << ", " << cd.y << ", " << cd.z << ") "; return out; } /// \brief Communication between colvars and LAMMPS /// (implementation of \link colvarproxy \endlink) class colvarproxy_lammps : public colvarproxy { // LAMMPS specific data objects and flags protected: // pointers to LAMMPS class instances class LAMMPS_NS::LAMMPS *_lmp; class LAMMPS_NS::RanPark *_random; // state of LAMMPS properties double t_target, my_timestep, my_boltzmann, my_angstrom; double bias_energy; int restart_every; int previous_step; bool first_timestep; bool total_force_requested; bool do_exit; // std::vector colvars_atoms; // std::vector colvars_atoms_ncopies; // std::vector positions; // std::vector total_forces; // std::vector applied_forces; // std::vector previous_applied_forces; std::vector atoms_types; MPI_Comm inter_comm; // MPI comm with 1 root proc from each world int inter_me, inter_num; // rank for the inter replica comm public: friend class cvm::atom; colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp, const char *, const char *, const int, const double, MPI_Comm); virtual ~colvarproxy_lammps(); void init(const char*); int setup(); // disable default and copy constructor private: colvarproxy_lammps() {}; colvarproxy_lammps(const colvarproxy_lammps &) {}; // methods for lammps to move data or trigger actions in the proxy public: void set_temperature(double t) { t_target = t; }; - bool need_total_forces() const { return total_force_requested; }; + bool total_forces_enabled() const { return total_force_requested; }; bool want_exit() const { return do_exit; }; // perform colvars computation. returns biasing energy double compute(); // dump status to string void serialize_status(std::string &); // set status from string bool deserialize_status(std::string &); // implementation of pure methods from base class public: inline cvm::real unit_angstrom() { return my_angstrom; }; inline cvm::real boltzmann() { return my_boltzmann; }; inline cvm::real temperature() { return t_target; }; inline cvm::real dt() { return my_timestep; }; // return _lmp->update->dt * _lmp->force->femtosecond; }; inline size_t restart_frequency() { return restart_every; }; void add_energy(cvm::real energy) { bias_energy += energy; }; void request_total_force(bool yesno) { total_force_requested = yesno; }; void log(std::string const &message); void error(std::string const &message); void fatal_error(std::string const &message); void exit(std::string const &message); cvm::rvector position_distance(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2); cvm::real position_dist2(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2); void select_closest_image(cvm::atom_pos &pos, cvm::atom_pos const &ref_pos); int backup_file(char const *filename); cvm::real rand_gaussian(void) { return _random->gaussian(); }; int init_atom(int atom_number); int check_atom_id(int atom_number); inline std::vector *modify_atom_types() { return &atoms_types; } // implementation of optional methods from base class public: #if defined(_OPENMP) // SMP support int smp_enabled(); int smp_colvars_loop(); int smp_biases_loop(); int smp_thread_id(); int smp_num_threads(); protected: omp_lock_t smp_lock_state; public: int smp_lock(); int smp_trylock(); int smp_unlock(); #endif // Multi-replica support // Indicate if multi-replica support is available and active virtual bool replica_enabled() { return (inter_comm != MPI_COMM_NULL); } // Index of this replica virtual int replica_index() { return inter_me; } // Total number of replica virtual int replica_num() { return inter_num; } // Synchronize replica virtual void replica_comm_barrier(); // Receive data from other replica virtual int replica_comm_recv(char* msg_data, int buf_len, int src_rep); // Send data to other replica virtual int replica_comm_send(char* msg_data, int msg_len, int dest_rep); }; #endif diff --git a/src/USER-COLVARS/fix_colvars.cpp b/src/USER-COLVARS/fix_colvars.cpp index 997762544..6a3613d76 100644 --- a/src/USER-COLVARS/fix_colvars.cpp +++ b/src/USER-COLVARS/fix_colvars.cpp @@ -1,940 +1,940 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Axel Kohlmeyer (Temple U) ------------------------------------------------------------------------- */ #include #include #include #include #include #include "fix_colvars.h" #include "atom.h" #include "comm.h" #include "domain.h" #include "error.h" #include "group.h" #include "memory.h" #include "modify.h" #include "random_park.h" #include "respa.h" #include "universe.h" #include "update.h" #include "citeme.h" #include "colvarproxy_lammps.h" static const char colvars_pub[] = "fix colvars command:\n\n" "@Article{fiorin13,\n" " author = {G.~Fiorin and M.{\\,}L.~Klein and J.~H{\\'e}nin},\n" " title = {Using collective variables to drive molecular" " dynamics simulations},\n" " journal = {Mol.~Phys.},\n" " year = 2013,\n" " note = {doi: 10.1080/00268976.2013.813594}\n" "}\n\n"; /* re-usable integer hash table code with static linkage. */ /** hash table top level data structure */ typedef struct inthash_t { struct inthash_node_t **bucket; /* array of hash nodes */ int size; /* size of the array */ int entries; /* number of entries in table */ int downshift; /* shift cound, used in hash function */ int mask; /* used to select bits for hashing */ } inthash_t; /** hash table node data structure */ typedef struct inthash_node_t { int data; /* data in hash node */ int key; /* key for hash lookup */ struct inthash_node_t *next; /* next node in hash chain */ } inthash_node_t; #define HASH_FAIL -1 #define HASH_LIMIT 0.5 /* initialize new hash table */ static void inthash_init(inthash_t *tptr, int buckets); /* lookup entry in hash table */ static int inthash_lookup(void *tptr, int key); /* insert an entry into hash table. */ static int inthash_insert(inthash_t *tptr, int key, int data); /* delete the hash table */ static void inthash_destroy(inthash_t *tptr); /************************************************************************ * integer hash code: ************************************************************************/ /* inthash() - Hash function returns a hash number for a given key. * tptr: Pointer to a hash table, key: The key to create a hash number for */ static int inthash(const inthash_t *tptr, int key) { int hashvalue; hashvalue = (((key*1103515249)>>tptr->downshift) & tptr->mask); if (hashvalue < 0) { hashvalue = 0; } return hashvalue; } /* * rebuild_table_int() - Create new hash table when old one fills up. * * tptr: Pointer to a hash table */ static void rebuild_table_int(inthash_t *tptr) { inthash_node_t **old_bucket, *old_hash, *tmp; int old_size, h, i; old_bucket=tptr->bucket; old_size=tptr->size; /* create a new table and rehash old buckets */ inthash_init(tptr, old_size<<1); for (i=0; inext; h=inthash(tptr, tmp->key); tmp->next=tptr->bucket[h]; tptr->bucket[h]=tmp; tptr->entries++; } /* while */ } /* for */ /* free memory used by old table */ free(old_bucket); return; } /* * inthash_init() - Initialize a new hash table. * * tptr: Pointer to the hash table to initialize * buckets: The number of initial buckets to create */ void inthash_init(inthash_t *tptr, int buckets) { /* make sure we allocate something */ if (buckets==0) buckets=16; /* initialize the table */ tptr->entries=0; tptr->size=2; tptr->mask=1; tptr->downshift=29; /* ensure buckets is a power of 2 */ while (tptr->sizesize<<=1; tptr->mask=(tptr->mask<<1)+1; tptr->downshift--; } /* while */ /* allocate memory for table */ tptr->bucket=(inthash_node_t **) calloc(tptr->size, sizeof(inthash_node_t *)); return; } /* * inthash_lookup() - Lookup an entry in the hash table and return a pointer to * it or HASH_FAIL if it wasn't found. * * tptr: Pointer to the hash table * key: The key to lookup */ int inthash_lookup(void *ptr, int key) { const inthash_t *tptr = (const inthash_t *) ptr; int h; inthash_node_t *node; /* find the entry in the hash table */ h=inthash(tptr, key); for (node=tptr->bucket[h]; node!=NULL; node=node->next) { if (node->key == key) break; } /* return the entry if it exists, or HASH_FAIL */ return(node ? node->data : HASH_FAIL); } /* * inthash_insert() - Insert an entry into the hash table. If the entry already * exists return a pointer to it, otherwise return HASH_FAIL. * * tptr: A pointer to the hash table * key: The key to insert into the hash table * data: A pointer to the data to insert into the hash table */ int inthash_insert(inthash_t *tptr, int key, int data) { int tmp; inthash_node_t *node; int h; /* check to see if the entry exists */ if ((tmp=inthash_lookup(tptr, key)) != HASH_FAIL) return(tmp); /* expand the table if needed */ while (tptr->entries>=HASH_LIMIT*tptr->size) rebuild_table_int(tptr); /* insert the new entry */ h=inthash(tptr, key); node=(struct inthash_node_t *) malloc(sizeof(inthash_node_t)); node->data=data; node->key=key; node->next=tptr->bucket[h]; tptr->bucket[h]=node; tptr->entries++; return HASH_FAIL; } /* * inthash_destroy() - Delete the entire table, and all remaining entries. * */ void inthash_destroy(inthash_t *tptr) { inthash_node_t *node, *last; int i; for (i=0; isize; i++) { node = tptr->bucket[i]; while (node != NULL) { last = node; node = node->next; free(last); } } /* free the entire array of buckets */ if (tptr->bucket != NULL) { free(tptr->bucket); memset(tptr, 0, sizeof(inthash_t)); } } /***************************************************************/ using namespace LAMMPS_NS; using namespace FixConst; // initialize static class members int FixColvars::instances=0; /*************************************************************** create class and parse arguments in LAMMPS script. Syntax: fix ID group-ID colvars [optional flags...] optional keyword value pairs: input (for restarting/continuing, defaults to NULL, but set to at end) output (defaults to 'out') seed (seed for RNG, defaults to '1966') tstat (label of thermostatting fix) ***************************************************************/ FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { if (narg < 4) error->all(FLERR,"Illegal fix colvars command: too few arguments"); if (instances > 0) error->all(FLERR,"Only one colvars fix can be active at a time"); ++instances; scalar_flag = 1; global_freq = 1; nevery = 1; extscalar = 1; restart_global = 1; me = comm->me; root2root = MPI_COMM_NULL; conf_file = strdup(arg[3]); rng_seed = 1966; unwrap_flag = 1; inp_name = NULL; out_name = NULL; tmp_name = NULL; /* parse optional arguments */ int argsdone = 4; while (argsdone < narg) { // we have keyword/value pairs. check if value is missing if (argsdone+1 == narg) error->all(FLERR,"Missing argument to keyword"); if (0 == strcmp(arg[argsdone], "input")) { inp_name = strdup(arg[argsdone+1]); } else if (0 == strcmp(arg[argsdone], "output")) { out_name = strdup(arg[argsdone+1]); } else if (0 == strcmp(arg[argsdone], "seed")) { rng_seed = force->inumeric(FLERR,arg[argsdone+1]); } else if (0 == strcmp(arg[argsdone], "unwrap")) { if (0 == strcmp(arg[argsdone+1], "yes")) { unwrap_flag = 1; } else if (0 == strcmp(arg[argsdone+1], "no")) { unwrap_flag = 0; } else { error->all(FLERR,"Incorrect fix colvars unwrap flag"); } } else if (0 == strcmp(arg[argsdone], "tstat")) { tmp_name = strdup(arg[argsdone+1]); } else { error->all(FLERR,"Unknown fix colvars parameter"); } ++argsdone; ++argsdone; } if (!out_name) out_name = strdup("out"); /* initialize various state variables. */ tstat_id = -1; energy = 0.0; nlevels_respa = 0; init_flag = 0; num_coords = 0; comm_buf = NULL; force_buf = NULL; proxy = NULL; idmap = NULL; /* storage required to communicate a single coordinate or force. */ size_one = sizeof(struct commdata); if (lmp->citeme) lmp->citeme->add(colvars_pub); } /********************************* * Clean up on deleting the fix. * *********************************/ FixColvars::~FixColvars() { memory->sfree(conf_file); memory->sfree(inp_name); memory->sfree(out_name); memory->sfree(tmp_name); memory->sfree(comm_buf); if (proxy) { delete proxy; inthash_t *hashtable = (inthash_t *)idmap; inthash_destroy(hashtable); delete hashtable; } if (root2root != MPI_COMM_NULL) MPI_Comm_free(&root2root); --instances; } /* ---------------------------------------------------------------------- */ int FixColvars::setmask() { int mask = 0; mask |= THERMO_ENERGY; mask |= MIN_POST_FORCE; mask |= POST_FORCE; mask |= POST_FORCE_RESPA; mask |= END_OF_STEP; return mask; } /* ---------------------------------------------------------------------- */ // initial checks for colvars run. void FixColvars::init() { if (atom->tag_enable == 0) error->all(FLERR,"Cannot use fix colvars without atom IDs"); if (atom->map_style == 0) error->all(FLERR,"Fix colvars requires an atom map, see atom_modify"); if ((me == 0) && (update->whichflag == 2)) error->warning(FLERR,"Using fix colvars with minimization"); if (strstr(update->integrate_style,"respa")) nlevels_respa = ((Respa *) update->integrate)->nlevels; } /* ---------------------------------------------------------------------- */ void FixColvars::one_time_init() { int i,tmp; if (init_flag) return; init_flag = 1; if (universe->nworlds > 1) { // create inter root communicator int color = 1; if (me == 0) color = 0; MPI_Comm_split(universe->uworld,color,universe->iworld,&root2root); } // create and initialize the colvars proxy if (me == 0) { if (screen) fputs("colvars: Creating proxy instance\n",screen); if (logfile) fputs("colvars: Creating proxy instance\n",logfile); #ifdef LAMMPS_BIGBIG if (screen) fputs("colvars: cannot handle atom ids > 2147483647\n",screen); if (logfile) fputs("colvars: cannot handle atom ids > 2147483647\n",logfile); #endif if (inp_name) { if (strcmp(inp_name,"NULL") == 0) { memory->sfree(inp_name); inp_name = NULL; } } // try to determine thermostat target temperature double t_target = 0.0; if (tmp_name) { if (strcmp(tmp_name,"NULL") == 0) { tstat_id = -1; } else { tstat_id = modify->find_fix(tmp_name); if (tstat_id < 0) error->one(FLERR,"Could not find tstat fix ID"); double *tt = (double*)modify->fix[tstat_id]->extract("t_target",tmp); if (tt) t_target = *tt; } } proxy = new colvarproxy_lammps(lmp,inp_name,out_name, rng_seed,t_target,root2root); proxy->init(conf_file); num_coords = (proxy->modify_atom_positions()->size()); } // send the list of all colvar atom IDs to all nodes. // also initialize and build hashtable on master. MPI_Bcast(&num_coords, 1, MPI_INT, 0, world); memory->create(taglist,num_coords,"colvars:taglist"); memory->create(force_buf,3*num_coords,"colvars:force_buf"); if (me == 0) { std::vector &tl = *(proxy->modify_atom_ids()); inthash_t *hashtable=new inthash_t; inthash_init(hashtable, num_coords); idmap = (void *)hashtable; for (i=0; i < num_coords; ++i) { taglist[i] = tl[i]; inthash_insert(hashtable, tl[i], i); } } MPI_Bcast(taglist, num_coords, MPI_LMP_TAGINT, 0, world); } /* ---------------------------------------------------------------------- */ void FixColvars::setup(int vflag) { const tagint * const tag = atom->tag; const int * const type = atom->type; int i,nme,tmp,ndata; int nlocal = atom->nlocal; MPI_Status status; MPI_Request request; one_time_init(); // determine size of comm buffer nme=0; for (i=0; i < num_coords; ++i) { const tagint k = atom->map(taglist[i]); if ((k >= 0) && (k < nlocal)) ++nme; } MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); memory->create(comm_buf,nmax,"colvars:comm_buf"); const double * const * const x = atom->x; const imageint * const image = atom->image; const double xprd = domain->xprd; const double yprd = domain->yprd; const double zprd = domain->zprd; const double xy = domain->xy; const double xz = domain->xz; const double yz = domain->yz; if (me == 0) { std::vector &id = *(proxy->modify_atom_ids()); std::vector &tp = *(proxy->modify_atom_types()); std::vector &cd = *(proxy->modify_atom_positions()); std::vector &of = *(proxy->modify_atom_total_forces()); std::vector &m = *(proxy->modify_atom_masses()); std::vector &q = *(proxy->modify_atom_charges()); // store coordinate data in holding array, clear old forces for (i=0; imap(taglist[i]); if ((k >= 0) && (k < nlocal)) { of[i].x = of[i].y = of[i].z = 0.0; if (unwrap_flag) { const int ix = (image[k] & IMGMASK) - IMGMAX; const int iy = (image[k] >> IMGBITS & IMGMASK) - IMGMAX; const int iz = (image[k] >> IMG2BITS) - IMGMAX; cd[i].x = x[k][0] + ix * xprd + iy * xy + iz * xz; cd[i].y = x[k][1] + iy * yprd + iz * yz; cd[i].z = x[k][2] + iz * zprd; } else { cd[i].x = x[k][0]; cd[i].y = x[k][1]; cd[i].z = x[k][2]; } if (atom->rmass_flag) { m[i] = atom->rmass[k]; } else { m[i] = atom->mass[type[k]]; } if (atom->q_flag) { q[i] = atom->q[k]; } } } // loop over procs to receive and apply remote data for (i=1; i < comm->nprocs; ++i) { int maxbuf = nmax*size_one; MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request); MPI_Send(&tmp, 0, MPI_INT, i, 0, world); MPI_Wait(&request, &status); MPI_Get_count(&status, MPI_BYTE, &ndata); ndata /= size_one; for (int k=0; kmap(taglist[i]); if ((k >= 0) && (k < nlocal)) { comm_buf[nme].tag = tag[k]; comm_buf[nme].type = type[k]; if (unwrap_flag) { const int ix = (image[k] & IMGMASK) - IMGMAX; const int iy = (image[k] >> IMGBITS & IMGMASK) - IMGMAX; const int iz = (image[k] >> IMG2BITS) - IMGMAX; comm_buf[nme].x = x[k][0] + ix * xprd + iy * xy + iz * xz; comm_buf[nme].y = x[k][1] + iy * yprd + iz * yz; comm_buf[nme].z = x[k][2] + iz * zprd; } else { comm_buf[nme].x = x[k][0]; comm_buf[nme].y = x[k][1]; comm_buf[nme].z = x[k][2]; } if (atom->rmass_flag) { comm_buf[nme].m = atom->rmass[k]; } else { comm_buf[nme].m = atom->mass[type[k]]; } if (atom->q_flag) { comm_buf[nme].q = atom->q[k]; } ++nme; } } /* blocking receive to wait until it is our turn to send data. */ MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, MPI_STATUS_IGNORE); MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world); } // run pre-run setup in colvarproxy if (me == 0) proxy->setup(); // initialize forces if (strstr(update->integrate_style,"verlet") || (update->whichflag == 2)) post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } /* ---------------------------------------------------------------------- */ /* Main colvars handler: * Send coodinates and add colvar forces to atoms. */ void FixColvars::post_force(int vflag) { // some housekeeping: update status of the proxy as needed. if (me == 0) { if (proxy->want_exit()) error->one(FLERR,"Run aborted on request from colvars module.\n"); if (tstat_id < 0) { proxy->set_temperature(0.0); } else { int tmp; // get thermostat target temperature from corresponding fix, // if the fix supports extraction. double *tt = (double *) modify->fix[tstat_id]->extract("t_target",tmp); if (tt) proxy->set_temperature(*tt); else proxy->set_temperature(0.0); } } const tagint * const tag = atom->tag; const double * const * const x = atom->x; double * const * const f = atom->f; const imageint * const image = atom->image; const double xprd = domain->xprd; const double yprd = domain->yprd; const double zprd = domain->zprd; const double xy = domain->xy; const double xz = domain->xz; const double yz = domain->yz; const int nlocal = atom->nlocal; /* check and potentially grow local communication buffers. */ int i,nmax_new,nme=0; for (i=0; i < num_coords; ++i) { const tagint k = atom->map(taglist[i]); if ((k >= 0) && (k < nlocal)) ++nme; } MPI_Allreduce(&nme,&nmax_new,1,MPI_INT,MPI_MAX,world); if (nmax_new > nmax) { nmax = nmax_new; memory->grow(comm_buf,nmax,"colvars:comm_buf"); } MPI_Status status; MPI_Request request; int tmp, ndata; if (me == 0) { std::vector &cd = *(proxy->modify_atom_positions()); // store coordinate data for (i=0; imap(taglist[i]); if ((k >= 0) && (k < nlocal)) { if (unwrap_flag) { const int ix = (image[k] & IMGMASK) - IMGMAX; const int iy = (image[k] >> IMGBITS & IMGMASK) - IMGMAX; const int iz = (image[k] >> IMG2BITS) - IMGMAX; cd[i].x = x[k][0] + ix * xprd + iy * xy + iz * xz; cd[i].y = x[k][1] + iy * yprd + iz * yz; cd[i].z = x[k][2] + iz * zprd; } else { cd[i].x = x[k][0]; cd[i].y = x[k][1]; cd[i].z = x[k][2]; } } } /* loop over procs to receive remote data */ for (i=1; i < comm->nprocs; ++i) { int maxbuf = nmax*size_one; MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request); MPI_Send(&tmp, 0, MPI_INT, i, 0, world); MPI_Wait(&request, &status); MPI_Get_count(&status, MPI_BYTE, &ndata); ndata /= size_one; for (int k=0; kmap(taglist[i]); if ((k >= 0) && (k < nlocal)) { comm_buf[nme].tag = tag[k]; if (unwrap_flag) { const int ix = (image[k] & IMGMASK) - IMGMAX; const int iy = (image[k] >> IMGBITS & IMGMASK) - IMGMAX; const int iz = (image[k] >> IMG2BITS) - IMGMAX; comm_buf[nme].x = x[k][0] + ix * xprd + iy * xy + iz * xz; comm_buf[nme].y = x[k][1] + iy * yprd + iz * yz; comm_buf[nme].z = x[k][2] + iz * zprd; } else { comm_buf[nme].x = x[k][0]; comm_buf[nme].y = x[k][1]; comm_buf[nme].z = x[k][2]; } ++nme; } } /* blocking receive to wait until it is our turn to send data. */ MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, MPI_STATUS_IGNORE); MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world); } //////////////////////////////////////////////////////////////////////// // call our workhorse and retrieve additional information. if (me == 0) { energy = proxy->compute(); - store_forces = proxy->need_total_forces(); + store_forces = proxy->total_forces_enabled(); } //////////////////////////////////////////////////////////////////////// // broadcast store_forces flag and energy data to all processors MPI_Bcast(&energy, 1, MPI_DOUBLE, 0, world); MPI_Bcast(&store_forces, 1, MPI_INT, 0, world); // broadcast and apply biasing forces if (me == 0) { std::vector &fo = *(proxy->modify_atom_new_colvar_forces()); double *fbuf = force_buf; for (int j=0; j < num_coords; ++j) { *fbuf++ = fo[j].x; *fbuf++ = fo[j].y; *fbuf++ = fo[j].z; } } MPI_Bcast(force_buf, 3*num_coords, MPI_DOUBLE, 0, world); for (int i=0; i < num_coords; ++i) { const tagint k = atom->map(taglist[i]); if ((k >= 0) && (k < nlocal)) { f[k][0] += force_buf[3*i+0]; f[k][1] += force_buf[3*i+1]; f[k][2] += force_buf[3*i+2]; } } } /* ---------------------------------------------------------------------- */ void FixColvars::min_post_force(int vflag) { post_force(vflag); } /* ---------------------------------------------------------------------- */ void FixColvars::post_force_respa(int vflag, int ilevel, int iloop) { /* only process colvar forces on the outmost RESPA level. */ if (ilevel == nlevels_respa-1) post_force(vflag); return; } /* ---------------------------------------------------------------------- */ void FixColvars::end_of_step() { if (store_forces) { const tagint * const tag = atom->tag; double * const * const f = atom->f; const int nlocal = atom->nlocal; /* check and potentially grow local communication buffers. */ int i,nmax_new,nme=0; for (i=0; i < num_coords; ++i) { const tagint k = atom->map(taglist[i]); if ((k >= 0) && (k < nlocal)) ++nme; } MPI_Allreduce(&nme,&nmax_new,1,MPI_INT,MPI_MAX,world); if (nmax_new > nmax) { nmax = nmax_new; memory->grow(comm_buf,nmax,"colvars:comm_buf"); } MPI_Status status; MPI_Request request; int tmp, ndata; if (me == 0) { // store old force data std::vector &of = *(proxy->modify_atom_total_forces()); for (i=0; imap(taglist[i]); if ((k >= 0) && (k < nlocal)) { const int j = inthash_lookup(idmap, tag[k]); if (j != HASH_FAIL) { of[j].x = f[k][0]; of[j].y = f[k][1]; of[j].z = f[k][2]; } } } /* loop over procs to receive remote data */ for (i=1; i < comm->nprocs; ++i) { int maxbuf = nmax*size_one; MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request); MPI_Send(&tmp, 0, MPI_INT, i, 0, world); MPI_Wait(&request, &status); MPI_Get_count(&status, MPI_BYTE, &ndata); ndata /= size_one; for (int k=0; kmap(taglist[i]); if ((k >= 0) && (k < nlocal)) { comm_buf[nme].tag = tag[k]; comm_buf[nme].x = f[k][0]; comm_buf[nme].y = f[k][1]; comm_buf[nme].z = f[k][2]; ++nme; } } /* blocking receive to wait until it is our turn to send data. */ MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, MPI_STATUS_IGNORE); MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world); } } } /* ---------------------------------------------------------------------- */ void FixColvars::write_restart(FILE *fp) { if (me == 0) { std::string rest_text(""); proxy->serialize_status(rest_text); const char *cvm_state = rest_text.c_str(); int len = strlen(cvm_state) + 1; // need to include terminating NULL byte. fwrite(&len,sizeof(int),1,fp); fwrite(cvm_state,1,len,fp); } } /* ---------------------------------------------------------------------- */ void FixColvars::restart(char *buf) { one_time_init(); if (me == 0) { std::string rest_text(buf); proxy->deserialize_status(rest_text); } } /* ---------------------------------------------------------------------- */ double FixColvars::compute_scalar() { return energy; } /* ---------------------------------------------------------------------- */ /* local memory usage. approximately. */ double FixColvars::memory_usage(void) { double bytes = (double) (num_coords * (2*sizeof(int)+3*sizeof(double))); bytes += (double) (nmax*size_one) + sizeof(this); return bytes; }