From c2430939809ee9d8424a7acad79e39bfaa0ca528 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Mon, 24 Oct 2016 17:05:47 -0400 Subject: [PATCH] Fix wall forces and subtractAppliedForce for extended-Lagrangian ABF --- lib/colvars/colvar.cpp | 20 +++++++----- lib/colvars/colvarbias_abf.cpp | 57 +++++++++++++++++----------------- lib/colvars/colvarbias_abf.h | 2 ++ 3 files changed, 43 insertions(+), 36 deletions(-) diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index dbf5fa5219..78aa2608cf 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -1144,9 +1144,9 @@ cvm::real colvar::update_forces_energy() // For a periodic colvar, both walls may be applicable at the same time // in which case we pick the closer one if ( (!is_enabled(f_cv_upper_wall)) || - (this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) { + (this->dist2(x_reported, lower_wall) < this->dist2(x_reported, upper_wall)) ) { - cvm::real const grad = this->dist2_lgrad(x, lower_wall); + cvm::real const grad = this->dist2_lgrad(x_reported, lower_wall); if (grad < 0.0) { fw = -0.5 * lower_wall_k * grad; f += fw; @@ -1157,7 +1157,7 @@ cvm::real colvar::update_forces_energy() } else { - cvm::real const grad = this->dist2_lgrad(x, upper_wall); + cvm::real const grad = this->dist2_lgrad(x_reported, upper_wall); if (grad > 0.0) { fw = -0.5 * upper_wall_k * grad; f += fw; @@ -1177,15 +1177,21 @@ cvm::real colvar::update_forces_energy() // atoms only feel the harmonic force // fr: bias force on extended variable (without harmonic spring), for output in trajectory // f_ext: total force on extended variable (including harmonic spring) - // f: - initially, external biasing force + // f: - initially, external biasing force (including wall forces) // - after this code block, colvar force to be applied to atomic coordinates, ie. spring force fr = f; f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x); f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x); - // The total force acting on the extended variable is f_ext - // This will be used in the next timestep - ft_reported = f_ext; + if (is_enabled(f_cv_subtract_applied_force)) { + // Report a "system" force without the biases on this colvar + // that is, just the spring force + ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x); + } else { + // The total force acting on the extended variable is f_ext + // This will be used in the next timestep + ft_reported = f_ext; + } // leapfrog: starting from x_i, f_i, v_(i-1/2) vr += (0.5 * dt) * f_ext / ext_mass; diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index e3690a4ead..1665ea61a9 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -80,6 +80,7 @@ int colvarbias_abf::init(std::string const &conf) if (update_bias) { // Request calculation of total force (which also checks for availability) + // TODO - change this to a dependency - needs ABF-specific features if(enable(f_cvb_get_total_force)) return cvm::get_error(); } @@ -133,6 +134,10 @@ int colvarbias_abf::init(std::string const &conf) // Data for eABF z-based estimator if (b_extended) { + // CZAR output files for stratified eABF + get_keyval(conf, "writeCZARwindowFile", b_czar_window_file, false, + colvarparse::parse_silent); + z_bin.assign(colvars.size(), 0); z_samples = new colvar_grid_count(colvars); z_samples->request_actual_value(); @@ -241,7 +246,7 @@ int colvarbias_abf::update() for (size_t i = 0; i < colvars.size(); i++) { // get total forces (lagging by 1 timestep) from colvars - // and subtract previous ABF force + // and subtract previous ABF force if necessary update_system_force(i); } gradients->acc_force(force_bin, system_force); @@ -457,28 +462,30 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app 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"); + cvm::error("Error opening eABF z-histogram 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"); + if (b_czar_window_file) { + std::string z_gradients_out_name = prefix + ".zgrad"; + cvm::ofstream z_gradients_os; + + 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 eABF z-gradient file " + z_gradients_out_name + " for writing"); + } + z_gradients->write_multicol(z_gradients_os); + z_gradients_os.close(); } - z_gradients->write_multicol(z_gradients_os); - z_gradients_os.close(); // Calculate CZAR estimator of gradients for (std::vector ix = czar_gradients->new_index(); @@ -490,6 +497,9 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app } } + std::string czar_gradients_out_name = prefix + ".czar.grad"; + cvm::ofstream czar_gradients_os; + 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()) { @@ -499,17 +509,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app 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"; + std::string czar_pmf_out_name = prefix + ".czar.pmf"; if (!append) cvm::backup_file(czar_pmf_out_name.c_str()); cvm::ofstream czar_pmf_os; // Do numerical integration and output a PMF @@ -520,8 +520,6 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app czar_pmf_os.close(); } } - - return; } @@ -559,7 +557,7 @@ void colvarbias_abf::read_gradients_samples() std::ifstream is; - cvm::log("Reading sample count from " + samples_in_name + " and gradients from " + gradients_in_name); + cvm::log("Reading sample count from " + samples_in_name + " and gradient 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); @@ -572,17 +570,18 @@ void colvarbias_abf::read_gradients_samples() is.close(); if (z_gradients) { - cvm::log("Reading z sample count from " + z_samples_in_name + " and z gradients from " + z_gradients_in_name); + // Read eABF z-averaged data for CZAR + cvm::log("Reading z-histogram from " + z_samples_in_name + " and z-gradient 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"); + if (!is.is_open()) cvm::error("Error opening eABF z-histogram 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"); + if (!is.is_open()) cvm::error("Error opening eABF z-gradient file " + z_gradients_in_name + " for reading"); z_gradients->read_multicol(is, true); is.close(); } diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h index a706bbb5ce..ee76240968 100644 --- a/lib/colvars/colvarbias_abf.h +++ b/lib/colvars/colvarbias_abf.h @@ -40,6 +40,8 @@ private: int output_freq; /// Write combined files with a history of all output data? bool b_history_files; + /// Write CZAR output file for stratified eABF (.zgrad) + bool b_czar_window_file; size_t history_freq; /// Cap applied biasing force?