From d928c23d21fc61e1ac715b6ef7fb0510929e7e08 Mon Sep 17 00:00:00 2001 From: fabsugar Date: Mon, 7 Nov 2022 15:36:10 -0500 Subject: [PATCH] Modified implementation of interval in calc_energy and calc_forces by using get_colvar_index_bound to treat vaules out of the boundary (which are projected inside the closest boundary). Additionally curr_values are now defined as a general variable to avoid define a new variable every time functions are called --- src/colvarbias_meta.cpp | 82 ++++++++++------------------------------- src/colvarbias_meta.h | 3 ++ src/colvargrid.h | 11 ++++++ 3 files changed, 33 insertions(+), 63 deletions(-) diff --git a/src/colvarbias_meta.cpp b/src/colvarbias_meta.cpp index afa032858..7cbdffe1e 100644 --- a/src/colvarbias_meta.cpp +++ b/src/colvarbias_meta.cpp @@ -575,6 +575,11 @@ int colvarbias_meta::init_interval_params(std::string const &conf) std::vector interval_ulimit_cv; size_t i; size_t j; + // allocate accessory variable of colvar values to be used in energy and forces calculation + curr_values.resize(num_variables()); + for (i = 0; i < num_variables(); i++) { + curr_values[i].type(variables(i)->value()); + } if (get_keyval(conf, "useHillsInterval", use_interval, use_interval)) { if (use_interval) { @@ -1220,11 +1225,7 @@ int colvarbias_meta::calc_energy(std::vector const *values) size_t i; int ii; - cvm::real up_bound_bin_value; - std::vector curr_values(num_variables()); - for (i = 0; i < num_variables(); i++) { - curr_values[i].type(variables(i)->value()); - } + std::vector curr_bin; curr_values = values ? *values : colvar_values; @@ -1235,25 +1236,19 @@ int colvarbias_meta::calc_energy(std::vector const *values) curr_values[i]=interval_llimit[ii]; } ii=which_int_ulimit_cv[i]; - if (ii>-1 && curr_values[i]>=interval_ulimit[ii] ) { - // check if upper border is out of the grid otherwise put it back on the grid - up_bound_bin_value=hills_energy->lower_boundaries[i].real_value+variables(i)->width*(0.5+cvm::floor((interval_ulimit[ii]-hills_energy->lower_boundaries[i].real_value)/variables(i)->width)); - //if (interval_ulimit[ii]==hills_energy->upper_boundaries[i].real_value){ - if (up_bound_bin_value>hills_energy->upper_boundaries[i].real_value) { - curr_values[i]=interval_ulimit[ii]-0.5*(variables(i)->width); // upper border is out of grid; in this way is in - } else { - curr_values[i]=interval_ulimit[ii]; - } + if (ii>-1 && curr_values[i]>interval_ulimit[ii] ) { + curr_values[i]=interval_ulimit[ii]; } } + curr_bin = hills_energy->get_colvars_index_bound(curr_values); + } else { + curr_bin = hills_energy->get_colvars_index(curr_values); } for (ir = 0; ir < replicas.size(); ir++) { replicas[ir]->bias_energy = 0.0; } - std::vector const curr_bin = hills_energy->get_colvars_index(curr_values); - if (hills_energy->index_ok(curr_bin)) { // index is within the grid: get the energy from there for (ir = 0; ir < replicas.size(); ir++) { @@ -1282,20 +1277,6 @@ int colvarbias_meta::calc_energy(std::vector const *values) // now include the hills that have not been binned yet (starting // from new_hills_begin) - if (use_interval) { - curr_values = values ? *values : colvar_values; - for (i = 0; i < num_variables(); i++) { - ii=which_int_llimit_cv[i]; - if (ii>-1 && curr_values[i]-1 && curr_values[i]>interval_ulimit[ii] ) { - curr_values[i]=interval_ulimit[ii]; - } - } - } - for (ir = 0; ir < replicas.size(); ir++) { calc_hills(replicas[ir]->new_hills_begin, replicas[ir]->hills.end(), @@ -1314,11 +1295,7 @@ int colvarbias_meta::calc_forces(std::vector const *values) { size_t ir = 0, ic = 0; int ii; - cvm::real up_bound_bin_value; - std::vector curr_values(num_variables()); - for (ic = 0; ic < num_variables(); ic++) { - curr_values[ic].type(variables(ic)->value()); - } + std::vector curr_bin; curr_values = values ? *values : colvar_values; std::vector add_force(num_variables()); for (ir = 0; ir < replicas.size(); ir++) { @@ -1338,20 +1315,18 @@ int colvarbias_meta::calc_forces(std::vector const *values) } ii=which_int_ulimit_cv[ic]; if (ii>-1) { - if ( curr_values[ic]>=interval_ulimit[ii] ) { + if ( curr_values[ic]>interval_ulimit[ii] ) { add_force[ic]=false; - up_bound_bin_value=hills_energy->lower_boundaries[ic].real_value+variables(ic)->width*(0.5+cvm::floor((interval_ulimit[ii]-hills_energy->lower_boundaries[ic].real_value)/variables(ic)->width)); - //if (interval_ulimit[ii]==hills_energy->upper_boundaries[ic].real_value){ - if (up_bound_bin_value>hills_energy->upper_boundaries[ic].real_value) { - curr_values[ic]=interval_ulimit[ii]-0.5*(variables(ic)->width); // upper border is out of grid; in this way is in - } else { - curr_values[ic]=interval_ulimit[ii]; - } + curr_values[ic]=interval_ulimit[ii]; } } } - std::vector const curr_bin = hills_energy->get_colvars_index(curr_values); + if (use_interval) { + curr_bin = hills_energy->get_colvars_index_bound(curr_values); + } else { + curr_bin = hills_energy->get_colvars_index(curr_values); + } if (hills_energy->index_ok(curr_bin)) { for (ir = 0; ir < replicas.size(); ir++) { @@ -1381,25 +1356,6 @@ int colvarbias_meta::calc_forces(std::vector const *values) // now include the hills that have not been binned yet (starting // from new_hills_begin) - if (use_interval) { - curr_values = values ? *values : colvar_values; - for (ic = 0; ic < num_variables(); ic++) { - ii=which_int_llimit_cv[ic]; - if (ii>-1) { - if ( curr_values[ic]-1) { - if ( curr_values[ic]>interval_ulimit[ii] ) { - curr_values[ic]=interval_ulimit[ii]; - } - } - } - } - - if (cvm::debug()) { cvm::log("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ diff --git a/src/colvarbias_meta.h b/src/colvarbias_meta.h index 3a8035513..2b9c7e3ab 100644 --- a/src/colvarbias_meta.h +++ b/src/colvarbias_meta.h @@ -252,6 +252,9 @@ class colvarbias_meta std::vector interval_llimit; std::vector interval_ulimit; + /// \brief Current value of colvars to be modifed for calculation of energy and forces with interval + std::vector curr_values; + /// Ensemble-biased metadynamics (EBmeta) flag bool ebmeta; diff --git a/src/colvargrid.h b/src/colvargrid.h index ef3d96734..5c476e575 100644 --- a/src/colvargrid.h +++ b/src/colvargrid.h @@ -625,6 +625,17 @@ template class colvar_grid : public colvarparse { /// \brief Get the bin indices corresponding to the provided values of /// the colvars and assign first or last bin if out of boundaries + inline std::vector const get_colvars_index_bound(std::vector const &values) const + { + std::vector index = new_index(); + for (size_t i = 0; i < nd; i++) { + index[i] = value_to_bin_scalar_bound(values[i], i); + } + return index; + } + + /// \brief Get the bin indices corresponding to the current values of + /// the colvars and assign first or last bin if out of boundaries inline std::vector const get_colvars_index_bound() const { std::vector index = new_index();