ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Update Colvars to version 2020-07-07

Browse Source 
This update contains several fixes and small new features or usability
improvements.  Descriptions and authorship information can be accessed from
the pull requests listed below.

Skip the zero-step also when multiple run commands are executed in sequence (@giacomofiorin)
https://github.com/Colvars/colvars/pull/357

Do not accumulate data at step 0 (@giacomofiorin)
https://github.com/Colvars/colvars/pull/345

Support for symmetry permutations of atoms in RMSD (@jhenin)
https://github.com/Colvars/colvars/pull/344

Detect new grid parameters (@jhenin)
https://github.com/Colvars/colvars/pull/341

Only access the output streams in non-threaded regions (@giacomofiorin)
https://github.com/Colvars/colvars/pull/338

Fix incomplete setting of default colvarsRestartFrequency (@giacomofiorin)
https://github.com/Colvars/colvars/pull/334

Fix typo (@e-kwsm)
https://github.com/Colvars/colvars/pull/333

Convert the input keyword to lowercase in read_state_data_key (@HanatoK)
https://github.com/Colvars/colvars/pull/332

Implement reflecting b.c. for ext Lagrangian (@jhenin)
https://github.com/Colvars/colvars/pull/329
This commit is contained in:
Giacomo Fiorin
2020-09-04 14:46:27 -04:00
parent d995ed0d87
commit 85c394453c
45 changed files with 3238 additions and 1582 deletions
Download Patch File Download Diff File Expand all files Collapse all files

273
lib/colvars/colvarbias_abf.cpp
View File

@ -8,7 +8,6 @@
// Colvars repository at GitHub.
#include "colvarmodule.h"
#include "colvarproxy.h"
#include "colvar.h"
#include "colvarbias_abf.h"
@ -29,8 +28,14 @@ colvarbias_abf::colvarbias_abf(char const *key)
last_gradients(NULL),
last_samples(NULL)
{
colvarproxy *proxy = cvm::main()->proxy;
if (!proxy->total_forces_same_step()) {
// Samples at step zero can not be collected
feature_states[f_cvb_step_zero_data].available = false;
}
}
int colvarbias_abf::init(std::string const &conf)
{
colvarbias::init(conf);
@ -71,8 +76,19 @@ int colvarbias_abf::init(std::string const &conf)
// full_samples - min_samples >= 1 is guaranteed
get_keyval(conf, "inputPrefix", input_prefix, std::vector<std::string>());
get_keyval(conf, "outputFreq", output_freq, cvm::restart_out_freq);
get_keyval(conf, "historyFreq", history_freq, 0);
if (history_freq != 0) {
if (output_freq == 0) {
cvm::error("Error: historyFreq must be a multiple of outputFreq.\n",
INPUT_ERROR);
} else {
if ((history_freq % output_freq) != 0) {
cvm::error("Error: historyFreq must be a multiple of outputFreq.\n",
INPUT_ERROR);
}
}
}
b_history_files = (history_freq > 0);
// shared ABF
@ -146,6 +162,7 @@ int colvarbias_abf::init(std::string const &conf)
for (i = 0; i < num_variables(); i++) {
if (max_force[i] < 0.0) {
cvm::error("Error: maxForce should be non-negative.");
return COLVARS_ERROR;
}
}
cap_force = true;
@ -184,25 +201,24 @@ int colvarbias_abf::init(std::string const &conf)
czar_gradients = new colvar_grid_gradient(colvars);
}
// For now, we integrate on-the-fly iff the grid is < 3D
if ( num_variables() <= 3 ) {
get_keyval(conf, "integrate", b_integrate, num_variables() <= 3); // Integrate for output if d<=3
if (b_integrate) {
// For now, we integrate on-the-fly iff the grid is < 3D
if ( num_variables() > 3 ) {
cvm::error("Error: cannot integrate free energy in dimension > 3.\n");
return COLVARS_ERROR;
}
pmf = new integrate_potential(colvars, gradients);
if ( b_CZAR_estimator ) {
czar_pmf = new integrate_potential(colvars, czar_gradients);
}
get_keyval(conf, "integrate", b_integrate, true); // Integrate for output
if ( num_variables() > 1 ) {
// Projected ABF
get_keyval(conf, "pABFintegrateFreq", pabf_freq, 0);
// Parameters for integrating initial (and final) gradient data
get_keyval(conf, "integrateInitMaxIterations", integrate_initial_iterations, 1e4);
get_keyval(conf, "integrateInitTol", integrate_initial_tol, 1e-6);
// for updating the integrated PMF on the fly
get_keyval(conf, "integrateMaxIterations", integrate_iterations, 100);
get_keyval(conf, "integrateTol", integrate_tol, 1e-4);
}
} else {
b_integrate = false;
// Parameters for integrating initial (and final) gradient data
get_keyval(conf, "integrateMaxIterations", integrate_iterations, 1e4, colvarparse::parse_silent);
get_keyval(conf, "integrateTol", integrate_tol, 1e-6, colvarparse::parse_silent);
// Projected ABF, updating the integrated PMF on the fly
get_keyval(conf, "pABFintegrateFreq", pabf_freq, 0, colvarparse::parse_silent);
get_keyval(conf, "pABFintegrateMaxIterations", pabf_integrate_iterations, 100, colvarparse::parse_silent);
get_keyval(conf, "pABFintegrateTol", pabf_integrate_tol, 1e-4, colvarparse::parse_silent);
}
// For shared ABF, we store a second set of grids.
@ -330,7 +346,7 @@ int colvarbias_abf::update()
force_bin = bin;
}
if (cvm::step_relative() > 0 || cvm::proxy->total_forces_same_step()) {
if (cvm::step_relative() > 0 || is_enabled(f_cvb_step_zero_data)) {
if (update_bias) {
// if (b_adiabatic_reweighting) {
@ -346,10 +362,10 @@ int colvarbias_abf::update()
// and subtract previous ABF force if necessary
update_system_force(i);
}
gradients->acc_force(force_bin, system_force);
if ( b_integrate ) {
pmf->update_div_neighbors(force_bin);
}
gradients->acc_force(force_bin, system_force);
if ( b_integrate ) {
pmf->update_div_neighbors(force_bin);
}
}
}
@ -370,10 +386,10 @@ int colvarbias_abf::update()
if ( b_integrate ) {
if ( pabf_freq && cvm::step_relative() % pabf_freq == 0 ) {
cvm::real err;
int iter = pmf->integrate(integrate_iterations, integrate_tol, err);
if ( iter == integrate_iterations ) {
cvm::log("Warning: PMF integration did not converge to " + cvm::to_str(integrate_tol)
+ " in " + cvm::to_str(integrate_iterations)
int iter = pmf->integrate(pabf_integrate_iterations, pabf_integrate_tol, err);
if ( iter == pabf_integrate_iterations ) {
cvm::log("Warning: PMF integration did not converge to " + cvm::to_str(pabf_integrate_tol)
+ " in " + cvm::to_str(pabf_integrate_iterations)
+ " steps. Residual error: " + cvm::to_str(err));
}
pmf->set_zero_minimum(); // TODO: do this only when necessary
@ -448,17 +464,6 @@ int colvarbias_abf::update()
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();
@ -487,10 +492,10 @@ int colvarbias_abf::update()
eabf_UI.update(cvm::step_absolute(), x, y);
}
/// Add the bias energy for 1D ABF
bias_energy = calc_energy(NULL);
/// Compute the bias energy
int error_code = calc_energy(NULL);
return COLVARS_OK;
return error_code;
}
@ -570,83 +575,71 @@ int colvarbias_abf::replica_share() {
last_samples->copy_grid(*samples);
shared_last_step = cvm::step_absolute();
if (b_integrate) {
// Update divergence to account for newly shared gradients
pmf->set_div();
}
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);
std::ostream *samples_os =
cvm::proxy->output_stream(samples_out_name, mode);
if (!samples_os) return;
samples->write_multicol(*samples_os);
cvm::proxy->close_output_stream(samples_out_name);
// In dimension higher than 2, dx is easier to handle and visualize
if (num_variables() > 2) {
std::string samples_dx_out_name = prefix + ".count.dx";
std::ostream *samples_dx_os = cvm::proxy->output_stream(samples_dx_out_name, mode);
if (!samples_os) return;
samples->write_opendx(*samples_dx_os);
*samples_dx_os << std::endl;
cvm::proxy->close_output_stream(samples_dx_out_name);
template <class T> int colvarbias_abf::write_grid_to_file(T const *grid,
std::string const &filename,
bool close) {
std::ostream *os = cvm::proxy->output_stream(filename);
if (!os) {
return cvm::error("Error opening file " + filename + " for writing.\n", COLVARS_ERROR | FILE_ERROR);
}
grid->write_multicol(*os);
if (close) {
cvm::proxy->close_output_stream(filename);
} else {
// Insert empty line between frames in history files
*os << std::endl;
cvm::proxy->flush_output_stream(os);
}
std::ostream *gradients_os =
cvm::proxy->output_stream(gradients_out_name, mode);
if (!gradients_os) return;
gradients->write_multicol(*gradients_os);
cvm::proxy->close_output_stream(gradients_out_name);
// In dimension higher than 2, dx is easier to handle and visualize
// but we cannot write multiple frames in a dx file now
// (could be implemented as multiple dx files)
if (num_variables() > 2 && close) {
std::string dx = filename + ".dx";
std::ostream *dx_os = cvm::proxy->output_stream(dx);
if (!dx_os) {
return cvm::error("Error opening file " + dx + " for writing.\n", COLVARS_ERROR | FILE_ERROR);
}
grid->write_opendx(*dx_os);
// if (close) {
cvm::proxy->close_output_stream(dx);
// }
// else {
// // TODO, decide convention for multiple datasets in dx file
// *dx_os << std::endl;
// dx_os->flush();
// }
}
return COLVARS_OK;
}
void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool close)
{
write_grid_to_file<colvar_grid_count>(samples, prefix + ".count", close);
write_grid_to_file<colvar_grid_gradient>(gradients, prefix + ".grad", close);
if (b_integrate) {
// Do numerical integration (to high precision) and output a PMF
cvm::real err;
pmf->integrate(integrate_initial_iterations, integrate_initial_tol, err);
pmf->integrate(integrate_iterations, integrate_tol, err);
pmf->set_zero_minimum();
std::string pmf_out_name = prefix + ".pmf";
std::ostream *pmf_os = cvm::proxy->output_stream(pmf_out_name, mode);
if (!pmf_os) return;
pmf->write_multicol(*pmf_os);
// In dimension higher than 2, dx is easier to handle and visualize
if (num_variables() > 2) {
std::string pmf_dx_out_name = prefix + ".pmf.dx";
std::ostream *pmf_dx_os = cvm::proxy->output_stream(pmf_dx_out_name, mode);
if (!pmf_dx_os) return;
pmf->write_opendx(*pmf_dx_os);
*pmf_dx_os << std::endl;
cvm::proxy->close_output_stream(pmf_dx_out_name);
}
*pmf_os << std::endl;
cvm::proxy->close_output_stream(pmf_out_name);
write_grid_to_file<colvar_grid_scalar>(pmf, prefix + ".pmf", close);
}
if (b_CZAR_estimator) {
// Write eABF CZAR-related quantities
std::string z_samples_out_name = prefix + ".zcount";
std::ostream *z_samples_os =
cvm::proxy->output_stream(z_samples_out_name, mode);
if (!z_samples_os) return;
z_samples->write_multicol(*z_samples_os);
cvm::proxy->close_output_stream(z_samples_out_name);
write_grid_to_file<colvar_grid_count>(z_samples, prefix + ".zcount", close);
if (b_czar_window_file) {
std::string z_gradients_out_name = prefix + ".zgrad";
std::ostream *z_gradients_os =
cvm::proxy->output_stream(z_gradients_out_name, mode);
if (!z_gradients_os) return;
z_gradients->write_multicol(*z_gradients_os);
cvm::proxy->close_output_stream(z_gradients_out_name);
write_grid_to_file<colvar_grid_gradient>(z_gradients, prefix + ".zgrad", close);
}
// Calculate CZAR estimator of gradients
@ -657,39 +650,15 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
- cvm::temperature() * cvm::boltzmann() * z_samples->log_gradient_finite_diff(ix, n), n);
}
}
std::string czar_gradients_out_name = prefix + ".czar.grad";
std::ostream *czar_gradients_os =
cvm::proxy->output_stream(czar_gradients_out_name, mode);
if (!czar_gradients_os) return;
czar_gradients->write_multicol(*czar_gradients_os);
cvm::proxy->close_output_stream(czar_gradients_out_name);
write_grid_to_file<colvar_grid_gradient>(czar_gradients, prefix + ".czar.grad", close);
if (b_integrate) {
// Do numerical integration (to high precision) and output a PMF
cvm::real err;
czar_pmf->set_div();
czar_pmf->integrate(integrate_initial_iterations, integrate_initial_tol, err);
czar_pmf->integrate(integrate_iterations, integrate_tol, err);
czar_pmf->set_zero_minimum();
std::string czar_pmf_out_name = prefix + ".czar.pmf";
std::ostream *czar_pmf_os = cvm::proxy->output_stream(czar_pmf_out_name, mode);
if (!czar_pmf_os) return;
czar_pmf->write_multicol(*czar_pmf_os);
// In dimension higher than 2, dx is easier to handle and visualize
if (num_variables() > 2) {
std::string czar_pmf_dx_out_name = prefix + ".czar.pmf.dx";
std::ostream *czar_pmf_dx_os = cvm::proxy->output_stream(czar_pmf_dx_out_name, mode);
if (!czar_pmf_dx_os) return;
czar_pmf->write_opendx(*czar_pmf_dx_os);
*czar_pmf_dx_os << std::endl;
cvm::proxy->close_output_stream(czar_pmf_dx_out_name);
}
*czar_pmf_os << std::endl;
cvm::proxy->close_output_stream(czar_pmf_out_name);
write_grid_to_file<colvar_grid_scalar>(czar_pmf, prefix + ".czar.pmf", close);
}
}
return;
@ -836,25 +805,56 @@ std::istream & colvarbias_abf::read_state_data(std::istream& is)
return is;
}
int colvarbias_abf::write_output_files()
{
if (cvm::debug()) {
cvm::log("ABF bias trying to write gradients and samples to disk");
}
if (shared_on && cvm::main()->proxy->replica_index() > 0) {
// No need to report the same data as replica 0, let it do the I/O job
return COLVARS_OK;
}
write_gradients_samples(output_prefix);
if (b_history_files) {
if ((cvm::step_absolute() % history_freq) == 0) {
write_gradients_samples(output_prefix + ".hist", false);
}
}
if (b_UI_estimator) {
eabf_UI.calc_pmf();
eabf_UI.write_files();
}
return COLVARS_OK;
}
int colvarbias_abf::calc_energy(std::vector<colvarvalue> const *values)
{
if (values) {
return cvm::error("colvarbias_abf::calc_energy() with an argument "
"is currently not implemented.\n",
COLVARS_NOT_IMPLEMENTED);
}
bias_energy = 0.0; // default value, overridden if a value can be calculated
if (num_variables() != 1) return 0.0;
if (num_variables() > 1 || values != NULL) {
// Use simple estimate: neglect effect of fullSamples,
// return value at center of bin
if (pmf != NULL) {
std::vector<int> const curr_bin = values ?
pmf->get_colvars_index(*values) :
pmf->get_colvars_index();
if (pmf->index_ok(curr_bin)) {
bias_energy = pmf->value(curr_bin);
}
}
return COLVARS_OK;
}
// Get the home bin.
int home0 = gradients->current_bin_scalar(0);
if (home0 < 0) return 0.0;
if (home0 < 0) return COLVARS_OK;
int gradient_len = (int)(gradients->number_of_points(0));
int home = (home0 < gradient_len) ? home0 : (gradient_len-1);
@ -886,5 +886,6 @@ int colvarbias_abf::calc_energy(std::vector<colvarvalue> const *values)
sum += fact*gradients->value(ix)/count*gradients->widths[0]*frac;
// The applied potential is the negative integral of force samples.
return -sum;
bias_energy = -sum;
return COLVARS_OK;
}