diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf
index 53dccba9c6c92d6b980ec05c9d97b89fc90991b7..daa1393269b953d569aff9846a135480dcdd42c7 100644
Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index d1a398ce4dad6392b4587bffb2508f9ae12b0d19..ce76b3b9eb2c8bd5f37b85af332611fa7ca963b3 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -1145,7 +1145,7 @@ int colvar::collect_cvc_values()
if (!cvcs[i]->is_enabled()) continue;
x += (cvcs[i])->sup_coeff *
( ((cvcs[i])->sup_np != 1) ?
- std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np) :
+ cvm::integer_power((cvcs[i])->value().real_value, (cvcs[i])->sup_np) :
(cvcs[i])->value().real_value );
}
} else {
@@ -1226,7 +1226,7 @@ int colvar::collect_cvc_gradients()
if (!cvcs[i]->is_enabled()) continue;
// Coefficient: d(a * x^n) = a * n * x^(n-1) * dx
cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real((cvcs[i])->sup_np) *
- std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1);
+ cvm::integer_power((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1);
for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
@@ -1593,9 +1593,9 @@ void colvar::communicate_forces()
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->is_enabled()) continue;
(cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff *
- cvm::real((cvcs[i])->sup_np) *
- (std::pow((cvcs[i])->value().real_value,
- (cvcs[i])->sup_np-1)) );
+ cvm::real((cvcs[i])->sup_np) *
+ (cvm::integer_power((cvcs[i])->value().real_value,
+ (cvcs[i])->sup_np-1)) );
}
} else {
diff --git a/lib/colvars/colvar_UIestimator.h b/lib/colvars/colvar_UIestimator.h
index 2fb250b71de01bbb9ffee981209dc14f98bf8b52..7fc7f870a10932d176525bb3c6c42d574a572fb1 100644
--- a/lib/colvars/colvar_UIestimator.h
+++ b/lib/colvars/colvar_UIestimator.h
@@ -45,7 +45,7 @@ namespace UIestimator {
this->width = width;
this->dimension = lowerboundary.size();
this->y_size = y_size; // keep in mind the internal (spare) matrix is stored in diagonal form
- this->y_total_size = int(pow(y_size, dimension) + EPSILON);
+ this->y_total_size = int(pow(double(y_size), dimension) + EPSILON);
// the range of the matrix is [lowerboundary, upperboundary]
x_total_size = 1;
@@ -121,7 +121,7 @@ namespace UIestimator {
int index = 0;
for (i = 0; i < dimension; i++) {
if (i + 1 < dimension)
- index += temp[i] * int(pow(y_size, dimension - i - 1) + EPSILON);
+ index += temp[i] * int(pow(double(y_size), dimension - i - 1) + EPSILON);
else
index += temp[i];
}
diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp
index 902ea4605d2f72aec3d79a546463ef2321d5180d..301e83e73015a91f14ff75a82c045954083d4074 100644
--- a/lib/colvars/colvarbias.cpp
+++ b/lib/colvars/colvarbias.cpp
@@ -524,9 +524,18 @@ int colvarbias_ti::update_system_forces(std::vector<colvarvalue> const
cvm::log("Updating system forces for bias "+this->name+"\n");
}
- // Collect total colvar forces from the previous step
+ colvarproxy *proxy = cvm::main()->proxy;
+
size_t i;
- if (cvm::step_relative() > 0) {
+
+ if (proxy->total_forces_same_step()) {
+ for (i = 0; i < num_variables(); i++) {
+ ti_bin[i] = ti_avg_forces->current_bin_scalar(i);
+ }
+ }
+
+ // Collect total colvar forces
+ if ((cvm::step_relative() > 0) || proxy->total_forces_same_step()) {
if (ti_avg_forces->index_ok(ti_bin)) {
for (i = 0; i < num_variables(); i++) {
if (variables(i)->is_enabled(f_cv_subtract_applied_force)) {
@@ -542,9 +551,11 @@ int colvarbias_ti::update_system_forces(std::vector<colvarvalue> const
}
}
- // Set the index to be used in the next iteration, when total forces come in
- for (i = 0; i < num_variables(); i++) {
- ti_bin[i] = ti_avg_forces->current_bin_scalar(i);
+ if (!proxy->total_forces_same_step()) {
+ // Set the index for use in the next iteration, when total forces come in
+ for (i = 0; i < num_variables(); i++) {
+ ti_bin[i] = ti_avg_forces->current_bin_scalar(i);
+ }
}
return COLVARS_OK;
diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp
index 5f5fe149784362e95807c2cfc7e6ef207ca90726..e4aea8eb86e7dcb6712fc6dc2c771adeeb6909f6 100644
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@@ -304,6 +304,10 @@ int colvarbias_abf::update()
// and subtract previous ABF force if necessary
update_system_force(i);
}
+ if (cvm::proxy->total_forces_same_step()) {
+ // e.g. in LAMMPS, total forces are current
+ force_bin = bin;
+ }
gradients->acc_force(force_bin, system_force);
}
if ( z_gradients && update_bias ) {
@@ -321,8 +325,11 @@ int colvarbias_abf::update()
}
}
- // save bin for next timestep
- force_bin = bin;
+ if (!cvm::proxy->total_forces_same_step()) {
+ // e.g. in NAMD, total forces will be available for next timestep
+ // hence we store the current colvar bin
+ force_bin = bin;
+ }
// Reset biasing forces from previous timestep
for (size_t i = 0; i < colvars.size(); i++) {
diff --git a/lib/colvars/colvarcomp_coordnums.cpp b/lib/colvars/colvarcomp_coordnums.cpp
index edfea96795b77d7dbd13fd174594d2b9ff007fc5..c34dc772157c64063627759370a8a936a2d46c13 100644
--- a/lib/colvars/colvarcomp_coordnums.cpp
+++ b/lib/colvars/colvarcomp_coordnums.cpp
@@ -18,6 +18,7 @@
+
template<bool calculate_gradients>
cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
int const &en,
@@ -32,8 +33,8 @@ cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
int const en2 = en/2;
int const ed2 = ed/2;
- cvm::real const xn = std::pow(l2, en2);
- cvm::real const xd = std::pow(l2, ed2);
+ cvm::real const xn = cvm::integer_power(l2, en2);
+ cvm::real const xd = cvm::integer_power(l2, ed2);
cvm::real const func = (1.0-xn)/(1.0-xd);
if (calculate_gradients) {
@@ -62,8 +63,8 @@ cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
int const en2 = en/2;
int const ed2 = ed/2;
- cvm::real const xn = std::pow(l2, en2);
- cvm::real const xd = std::pow(l2, ed2);
+ cvm::real const xn = cvm::integer_power(l2, en2);
+ cvm::real const xd = cvm::integer_power(l2, ed2);
cvm::real const func = (1.0-xn)/(1.0-xd);
if (calculate_gradients) {
@@ -117,11 +118,17 @@ colvar::coordnum::coordnum(std::string const &conf)
if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
}
- get_keyval(conf, "expNumer", en, int(6) );
- get_keyval(conf, "expDenom", ed, int(12));
+ get_keyval(conf, "expNumer", en, 6);
+ get_keyval(conf, "expDenom", ed, 12);
if ( (en%2) || (ed%2) ) {
- cvm::error("Error: odd exponents provided, can only use even ones.\n", INPUT_ERROR);
+ cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
+ INPUT_ERROR);
+ }
+
+ if ( (en <= 0) || (ed <= 0) ) {
+ cvm::error("Error: negative exponent(s) provided.\n",
+ INPUT_ERROR);
}
if (!is_enabled(f_cvc_pbc_minimum_image)) {
@@ -256,8 +263,13 @@ colvar::h_bond::h_bond(std::string const &conf)
get_keyval(conf, "expDenom", ed, 8);
if ( (en%2) || (ed%2) ) {
- cvm::error("Error: odd exponents provided, can only use even ones.\n");
- return;
+ cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
+ INPUT_ERROR);
+ }
+
+ if ( (en <= 0) || (ed <= 0) ) {
+ cvm::error("Error: negative exponent(s) provided.\n",
+ INPUT_ERROR);
}
if (cvm::debug())
@@ -324,12 +336,18 @@ colvar::selfcoordnum::selfcoordnum(std::string const &conf)
group1 = parse_group(conf, "group1");
get_keyval(conf, "cutoff", r0, cvm::real(4.0 * cvm::unit_angstrom()));
- get_keyval(conf, "expNumer", en, int(6) );
- get_keyval(conf, "expDenom", ed, int(12));
+ get_keyval(conf, "expNumer", en, 6);
+ get_keyval(conf, "expDenom", ed, 12);
+
if ( (en%2) || (ed%2) ) {
- cvm::error("Error: odd exponents provided, can only use even ones.\n");
- return;
+ cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
+ INPUT_ERROR);
+ }
+
+ if ( (en <= 0) || (ed <= 0) ) {
+ cvm::error("Error: negative exponent(s) provided.\n",
+ INPUT_ERROR);
}
if (!is_enabled(f_cvc_pbc_minimum_image)) {
@@ -407,12 +425,17 @@ colvar::groupcoordnum::groupcoordnum(std::string const &conf)
if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
}
- get_keyval(conf, "expNumer", en, int(6) );
- get_keyval(conf, "expDenom", ed, int(12));
+ get_keyval(conf, "expNumer", en, 6);
+ get_keyval(conf, "expDenom", ed, 12);
if ( (en%2) || (ed%2) ) {
- cvm::error("Error: odd exponents provided, can only use even ones.\n");
- return;
+ cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
+ INPUT_ERROR);
+ }
+
+ if ( (en <= 0) || (ed <= 0) ) {
+ cvm::error("Error: negative exponent(s) provided.\n",
+ INPUT_ERROR);
}
if (!is_enabled(f_cvc_pbc_minimum_image)) {
@@ -444,8 +467,8 @@ cvm::real colvar::groupcoordnum::switching_function(cvm::real const &r0,
int const en2 = en/2;
int const ed2 = ed/2;
- cvm::real const xn = std::pow(l2, en2);
- cvm::real const xd = std::pow(l2, ed2);
+ cvm::real const xn = cvm::integer_power(l2, en2);
+ cvm::real const xd = cvm::integer_power(l2, ed2);
cvm::real const func = (1.0-xn)/(1.0-xd);
if (calculate_gradients) {
@@ -477,8 +500,8 @@ cvm::real colvar::groupcoordnum::switching_function(cvm::rvector const &r0_vec,
int const en2 = en/2;
int const ed2 = ed/2;
- cvm::real const xn = std::pow(l2, en2);
- cvm::real const xd = std::pow(l2, ed2);
+ cvm::real const xn = cvm::integer_power(l2, en2);
+ cvm::real const xd = cvm::integer_power(l2, ed2);
cvm::real const func = (1.0-xn)/(1.0-xd);
if (calculate_gradients) {
diff --git a/lib/colvars/colvarcomp_protein.cpp b/lib/colvars/colvarcomp_protein.cpp
index b8fc96cfad140177e5adaa2ceccea1332b5702c7..91e47f13d90f22885d318eff8f9dc43a8b7fc21f 100644
--- a/lib/colvars/colvarcomp_protein.cpp
+++ b/lib/colvars/colvarcomp_protein.cpp
@@ -150,8 +150,8 @@ void colvar::alpha_angles::calc_value()
(theta[i])->calc_value();
cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol;
- cvm::real const f = ( (1.0 - std::pow(t, (int) 2)) /
- (1.0 - std::pow(t, (int) 4)) );
+ cvm::real const f = ( (1.0 - (t*t)) /
+ (1.0 - (t*t*t*t)) );
x.real_value += theta_norm * f;
@@ -202,12 +202,12 @@ void colvar::alpha_angles::apply_force(colvarvalue const &force)
for (size_t i = 0; i < theta.size(); i++) {
cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol;
- cvm::real const f = ( (1.0 - std::pow(t, (int) 2)) /
- (1.0 - std::pow(t, (int) 4)) );
+ cvm::real const f = ( (1.0 - (t*t)) /
+ (1.0 - (t*t*t*t)) );
cvm::real const dfdt =
- 1.0/(1.0 - std::pow(t, (int) 4)) *
- ( (-2.0 * t) + (-1.0*f)*(-4.0 * std::pow(t, (int) 3)) );
+ 1.0/(1.0 - (t*t*t*t)) *
+ ( (-2.0 * t) + (-1.0*f)*(-4.0 * (t*t*t)) );
(theta[i])->apply_force(theta_norm *
dfdt * (1.0/theta_tol) *
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 3f13ae88cf07cb438db71f5716f8bc465090c0ae..14e5d56701f75730c7bb1ce07764c89edd0b665b 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -83,6 +83,15 @@ public:
/// Defining an abstract real number allows to switch precision
typedef double real;
+
+ /// Override std::pow with a product for n positive integer
+ static inline real integer_power(real x, int n)
+ {
+ real result = 1.0;
+ for (int i = 0; i < n; i++) result *= x;
+ return result;
+ }
+
/// Residue identifier
typedef int residue_id;
diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h
index 59a15d7ac848b2820ea92b0ba8ae346a32bde036..a92a776f8a3abc88700794228bcf4a1e0ad8da51 100644
--- a/lib/colvars/colvars_version.h
+++ b/lib/colvars/colvars_version.h
@@ -1,5 +1,5 @@
#ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2017-10-19"
+#define COLVARS_VERSION "2017-10-20"
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars
diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp
index 428fe1a4b1d864c725cd7d5b09652216a03091a3..b604606d4628911138f2553215db9d7b3c306630 100644
--- a/lib/colvars/colvartypes.cpp
+++ b/lib/colvars/colvartypes.cpp
@@ -312,7 +312,7 @@ void colvarmodule::rotation::diagonalize_matrix(cvm::matrix2d<cvm::real> &S,
cvm::real norm2 = 0.0;
size_t i;
for (i = 0; i < 4; i++) {
- norm2 += std::pow(S_eigvec[ie][i], int(2));
+ norm2 += S_eigvec[ie][i] * S_eigvec[ie][i];
}
cvm::real const norm = std::sqrt(norm2);
for (i = 0; i < 4; i++) {
diff --git a/src/USER-COLVARS/colvarproxy_lammps.cpp b/src/USER-COLVARS/colvarproxy_lammps.cpp
index ad2366f3d4963162e86f2ba10053aedf1b671455..c5b9e5a60c860355a4a5960089b8f288b9566841 100644
--- a/src/USER-COLVARS/colvarproxy_lammps.cpp
+++ b/src/USER-COLVARS/colvarproxy_lammps.cpp
@@ -137,7 +137,7 @@ void colvarproxy_lammps::init(const char *conf_file)
colvars = new colvarmodule(this);
cvm::log("Using LAMMPS interface, version "+
- cvm::to_str(COLVARPROXY_VERSION)+".\n");
+ cvm::to_str(COLVARPROXY_VERSION)+".\n");
my_angstrom = _lmp->force->angstrom;
my_boltzmann = _lmp->force->boltz;
@@ -149,7 +149,8 @@ void colvarproxy_lammps::init(const char *conf_file)
colvars->setup_output();
if (_lmp->update->ntimestep != 0) {
- cvm::log("Initializing step number as firstTimestep.\n");
+ cvm::log("Setting initial step number from LAMMPS: "+
+ cvm::to_str(_lmp->update->ntimestep)+"\n");
colvars->it = colvars->it_restart = _lmp->update->ntimestep;
}
@@ -166,7 +167,6 @@ colvarproxy_lammps::~colvarproxy_lammps()
{
delete _random;
if (colvars != NULL) {
- colvars->write_output_files();
delete colvars;
colvars = NULL;
}
@@ -182,10 +182,18 @@ int colvarproxy_lammps::setup()
// trigger colvars computation
double colvarproxy_lammps::compute()
{
+ if (cvm::debug()) {
+ log(std::string(cvm::line_marker)+
+ "colvarproxy_lammps step no. "+
+ cvm::to_str(_lmp->update->ntimestep)+" [first - last = "+
+ cvm::to_str(_lmp->update->beginstep)+" - "+
+ cvm::to_str(_lmp->update->endstep)+"]\n");
+ }
+
if (first_timestep) {
first_timestep = false;
} else {
- // Use the time step number inherited from LAMMPS
+ // Use the time step number from LAMMPS Update object
if ( _lmp->update->ntimestep - previous_step == 1 )
colvars->it++;
// Other cases could mean:
@@ -233,8 +241,12 @@ void colvarproxy_lammps::serialize_status(std::string &rst)
std::ostringstream os;
colvars->write_restart(os);
rst = os.str();
+}
- // TODO separate this as its own function?
+void colvarproxy_lammps::write_output_files()
+{
+ // TODO skip output if undefined
+ colvars->write_restart_file(cvm::output_prefix()+".colvars.state");
colvars->write_output_files();
}
diff --git a/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h
index 849af410f216894309d98ba158a4b3e2d1012d9f..af2aa04dfc890a55443588435206279280581b6c 100644
--- a/src/USER-COLVARS/colvarproxy_lammps.h
+++ b/src/USER-COLVARS/colvarproxy_lammps.h
@@ -100,6 +100,10 @@ class colvarproxy_lammps : public colvarproxy {
// set status from string
bool deserialize_status(std::string &);
+ // Write files expected from Colvars (called by post_run())
+ void write_output_files();
+
+
// implementation of pure methods from base class
public:
diff --git a/src/USER-COLVARS/colvarproxy_lammps_version.h b/src/USER-COLVARS/colvarproxy_lammps_version.h
index e89e00ba039e3c03b32557d28f857216780d6519..45ecea867f17c7fdf2b5bbad780ac860287e8836 100644
--- a/src/USER-COLVARS/colvarproxy_lammps_version.h
+++ b/src/USER-COLVARS/colvarproxy_lammps_version.h
@@ -1,5 +1,5 @@
#ifndef COLVARPROXY_VERSION
-#define COLVARPROXY_VERSION "2017-10-11"
+#define COLVARPROXY_VERSION "2017-10-20"
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars
diff --git a/src/USER-COLVARS/fix_colvars.cpp b/src/USER-COLVARS/fix_colvars.cpp
index 5897b9c5f0191ae2bb77d3ca1870fa595d987a0c..956ba6498a518dacc28f4848788fd55427db839e 100644
--- a/src/USER-COLVARS/fix_colvars.cpp
+++ b/src/USER-COLVARS/fix_colvars.cpp
@@ -379,6 +379,7 @@ int FixColvars::setmask()
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= END_OF_STEP;
+ mask |= POST_RUN;
return mask;
}
@@ -935,6 +936,15 @@ void FixColvars::restart(char *buf)
/* ---------------------------------------------------------------------- */
+void FixColvars::post_run()
+{
+ if (me == 0) {
+ proxy->write_output_files();
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
double FixColvars::compute_scalar()
{
return energy;
diff --git a/src/USER-COLVARS/fix_colvars.h b/src/USER-COLVARS/fix_colvars.h
index c00b18aa4668b56ea8a1a82a729a0d02a9b0433a..509eca5de35892b3094cf5749b7d6d088d7b0d18 100644
--- a/src/USER-COLVARS/fix_colvars.h
+++ b/src/USER-COLVARS/fix_colvars.h
@@ -56,6 +56,7 @@ class FixColvars : public Fix {
virtual void post_force(int);
virtual void post_force_respa(int, int, int);
virtual void end_of_step();
+ virtual void post_run();
virtual double compute_scalar();
virtual double memory_usage();