diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf
index a14d93cd69e7a2aced22919a7755d66714405d47..ad15752107cc006471cc4b3e7a400aa1658f26e2 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/README b/lib/colvars/README
index ce1d3199742f558cbc56c7883c557b5ea05e6132..5df9612dfa3b6f73af79137701fdfae4bb90158b 100644
--- a/lib/colvars/README
+++ b/lib/colvars/README
@@ -32,7 +32,7 @@ where Makefile.g++ uses the GNU C++ compiler and is a good template to start.
**Optional**: if you use the Install.py script provided in this folder, you
can give the machine name as the '-m' argument. This can be the suffix of one
-of the files from either this folder, or from src/MAKE.
+of the files from either this folder, or from src/MAKE/MACHINES.
*This is only supported by the Install.py within the lib/colvars folder*.
When you are done building this library, two files should
@@ -53,10 +53,10 @@ settings in Makefile.common should work.
For the reference manual see:
http://colvars.github.io/colvars-refman-lammps
-A copy of reference manual is also in:
+A copy of the reference manual is also in:
doc/PDF/colvars-refman-lammps.pdf
-Also included is a Doxygen-based developer documentation:
+Also available is a Doxygen-based developer documentation:
http://colvars.github.io/doxygen/html/
The reference article is:
diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h
index 6113e1678bc9f3252f672c03d6ec62c4aae52bc8..dfa9e093a537b797284909e2d7a468571cb8cb98 100644
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@@ -88,7 +88,12 @@ public:
static std::vector<feature *> cv_features;
/// \brief Implementation of the feature list accessor for colvar
- std::vector<feature *> &features() {
+ virtual const std::vector<feature *> &features()
+ {
+ return cv_features;
+ }
+ virtual std::vector<feature *> &modify_features()
+ {
return cv_features;
}
diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h
index dba2890abcc6e7c94d0a78d8285d786595b9c023..6113fb38a938b37a56d54825c51d43252d5b54db 100644
--- a/lib/colvars/colvaratoms.h
+++ b/lib/colvars/colvaratoms.h
@@ -206,7 +206,12 @@ public:
static std::vector<feature *> ag_features;
/// \brief Implementation of the feature list accessor for atom group
- virtual std::vector<feature *> &features() {
+ virtual const std::vector<feature *> &features()
+ {
+ return ag_features;
+ }
+ virtual std::vector<feature *> &modify_features()
+ {
return ag_features;
}
diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp
index e437466be99e89656d803bae042d6e2e49cf923b..636727ca39b4f90b8b9686a64b51845857843c13 100644
--- a/lib/colvars/colvarbias.cpp
+++ b/lib/colvars/colvarbias.cpp
@@ -384,6 +384,7 @@ std::ostream & colvarbias::write_traj(std::ostream &os)
os << " ";
if (b_output_energy)
os << " "
+ << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
<< bias_energy;
return os;
}
diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h
index 205e761cfcb43e22c9d57cf729c5960e5b4ade24..a147cd3210486cc74d80f5e0c0e5f859e0ea1b70 100644
--- a/lib/colvars/colvarbias.h
+++ b/lib/colvars/colvarbias.h
@@ -175,7 +175,11 @@ public:
static std::vector<feature *> cvb_features;
/// \brief Implementation of the feature list accessor for colvarbias
- virtual std::vector<feature *> &features()
+ virtual const std::vector<feature *> &features()
+ {
+ return cvb_features;
+ }
+ virtual std::vector<feature *> &modify_features()
{
return cvb_features;
}
diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp
index bb6d6164e5cbb1abbbaaeaeb099471b364920a1e..6879190968f7921747d8cf5b0b157d5149ab0869 100644
--- a/lib/colvars/colvarbias_restraint.cpp
+++ b/lib/colvars/colvarbias_restraint.cpp
@@ -99,12 +99,9 @@ int colvarbias_restraint_centers::init(std::string const &conf)
if (null_centers) {
// try to initialize the restraint centers for the first time
colvar_centers.resize(num_variables());
- colvar_centers_raw.resize(num_variables());
for (i = 0; i < num_variables(); i++) {
colvar_centers[i].type(variables(i)->value());
colvar_centers[i].reset();
- colvar_centers_raw[i].type(variables(i)->value());
- colvar_centers_raw[i].reset();
}
}
@@ -113,7 +110,6 @@ int colvarbias_restraint_centers::init(std::string const &conf)
if (cvm::debug()) {
cvm::log("colvarbias_restraint: parsing initial centers, i = "+cvm::to_str(i)+".\n");
}
- colvar_centers_raw[i] = colvar_centers[i];
colvar_centers[i].apply_constraints();
}
null_centers = false;
@@ -141,8 +137,6 @@ int colvarbias_restraint_centers::change_configuration(std::string const &conf)
for (size_t i = 0; i < num_variables(); i++) {
colvar_centers[i].type(variables(i)->value());
colvar_centers[i].apply_constraints();
- colvar_centers_raw[i].type(variables(i)->value());
- colvar_centers_raw[i] = colvar_centers[i];
}
}
return COLVARS_OK;
@@ -232,7 +226,6 @@ int colvarbias_restraint_moving::set_state_params(std::string const &conf)
{
if (b_chg_centers || b_chg_force_k) {
if (target_nstages) {
- // cvm::log ("Reading current stage from the restart.\n");
if (!get_keyval(conf, "stage", stage))
cvm::error("Error: current stage is missing from the restart.\n");
}
@@ -265,100 +258,127 @@ int colvarbias_restraint_centers_moving::init(std::string const &conf)
size_t i;
if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) {
- if (colvar_centers.size() != num_variables()) {
+ if (target_centers.size() != num_variables()) {
cvm::error("Error: number of target centers does not match "
- "that of collective variables.\n");
+ "that of collective variables.\n", INPUT_ERROR);
}
b_chg_centers = true;
for (i = 0; i < target_centers.size(); i++) {
target_centers[i].apply_constraints();
+ centers_incr.push_back(colvar_centers[i]);
+ centers_incr[i].reset();
}
}
if (b_chg_centers) {
- // parse moving restraint options
+ // parse moving schedule options
colvarbias_restraint_moving::init(conf);
+ if (initial_centers.size() == 0) {
+ // One-time init
+ initial_centers = colvar_centers;
+ }
+ // Call to check that the definition is correct
+ for (i = 0; i < num_variables(); i++) {
+ colvarvalue const midpoint =
+ colvarvalue::interpolate(initial_centers[i],
+ target_centers[i],
+ 0.5);
+ }
} else {
target_centers.clear();
return COLVARS_OK;
}
get_keyval(conf, "outputCenters", b_output_centers, b_output_centers);
- get_keyval(conf, "outputAccumulatedWork", b_output_acc_work, b_output_acc_work);
+ get_keyval(conf, "outputAccumulatedWork", b_output_acc_work,
+ b_output_acc_work); // TODO this conflicts with stages
return COLVARS_OK;
}
-int colvarbias_restraint_centers_moving::update()
+int colvarbias_restraint_centers_moving::update_centers(cvm::real lambda)
{
- if (b_chg_centers) {
+ if (cvm::debug()) {
+ cvm::log("Updating centers for the restraint bias \""+
+ this->name+"\": "+cvm::to_str(colvar_centers)+".\n");
+ }
+ size_t i;
+ for (i = 0; i < num_variables(); i++) {
+ colvarvalue const c_new = colvarvalue::interpolate(initial_centers[i],
+ target_centers[i],
+ lambda);
+ centers_incr[i] = (c_new).dist2_grad(colvar_centers[i]);
+ colvar_centers[i] = c_new;
+ variables(i)->wrap(colvar_centers[i]);
+ }
+ if (cvm::debug()) {
+ cvm::log("New centers for the restraint bias \""+
+ this->name+"\": "+cvm::to_str(colvar_centers)+".\n");
+ }
+ return cvm::get_error();
+}
- if (cvm::debug()) {
- cvm::log("Updating centers for the restraint bias \""+
- this->name+"\": "+cvm::to_str(colvar_centers)+".\n");
- }
- if (!centers_incr.size()) {
- // if this is the first calculation, calculate the advancement
- // at each simulation step (or stage, if applicable)
- // (take current stage into account: it can be non-zero
- // if we are restarting a staged calculation)
- centers_incr.resize(num_variables());
- for (size_t i = 0; i < num_variables(); i++) {
- centers_incr[i].type(variables(i)->value());
- centers_incr[i] = (target_centers[i] - colvar_centers_raw[i]) /
- cvm::real( target_nstages ? (target_nstages - stage) :
- (target_nsteps - cvm::step_absolute()));
- }
- if (cvm::debug()) {
- cvm::log("Center increment for the restraint bias \""+
- this->name+"\": "+cvm::to_str(centers_incr)+" at stage "+cvm::to_str(stage)+ ".\n");
- }
- }
+int colvarbias_restraint_centers_moving::update()
+{
+ if (b_chg_centers) {
if (target_nstages) {
- if ((cvm::step_relative() > 0)
- && (cvm::step_absolute() % target_nsteps) == 0
- && stage < target_nstages) {
-
+ // Staged update
+ if (stage <= target_nstages) {
+ if ((cvm::step_relative() > 0) &&
+ ((cvm::step_absolute() % target_nsteps) == 1)) {
+ cvm::real const lambda =
+ cvm::real(stage)/cvm::real(target_nstages);
+ update_centers(lambda);
+ stage++;
+ cvm::log("Moving restraint \"" + this->name +
+ "\" stage " + cvm::to_str(stage) +
+ " : setting centers to " + cvm::to_str(colvar_centers) +
+ " at step " + cvm::to_str(cvm::step_absolute()));
+ } else {
+ for (size_t i = 0; i < num_variables(); i++) {
+ centers_incr[i].reset();
+ }
+ }
+ }
+ } else {
+ // Continuous update
+ if (cvm::step_absolute() <= target_nsteps) {
+ cvm::real const lambda =
+ cvm::real(cvm::step_absolute())/cvm::real(target_nsteps);
+ update_centers(lambda);
+ } else {
for (size_t i = 0; i < num_variables(); i++) {
- colvar_centers_raw[i] += centers_incr[i];
- colvar_centers[i] = colvar_centers_raw[i];
- variables(i)->wrap(colvar_centers[i]);
- colvar_centers[i].apply_constraints();
+ centers_incr[i].reset();
}
- stage++;
- cvm::log("Moving restraint \"" + this->name +
- "\" stage " + cvm::to_str(stage) +
- " : setting centers to " + cvm::to_str(colvar_centers) +
- " at step " + cvm::to_str(cvm::step_absolute()));
}
- } else if ((cvm::step_relative() > 0) && (cvm::step_absolute() <= target_nsteps)) {
- // move the restraint centers in the direction of the targets
- // (slow growth)
+ }
+
+ if (cvm::step_relative() == 0) {
for (size_t i = 0; i < num_variables(); i++) {
- colvar_centers_raw[i] += centers_incr[i];
- colvar_centers[i] = colvar_centers_raw[i];
- variables(i)->wrap(colvar_centers[i]);
- colvar_centers[i].apply_constraints();
+ // finite differences are undefined when restarting
+ centers_incr[i].reset();
}
}
if (cvm::debug()) {
- cvm::log("New centers for the restraint bias \""+
- this->name+"\": "+cvm::to_str(colvar_centers)+".\n");
+ cvm::log("Center increment for the restraint bias \""+
+ this->name+"\": "+cvm::to_str(centers_incr)+
+ " at stage "+cvm::to_str(stage)+ ".\n");
}
}
- return COLVARS_OK;
+ return cvm::get_error();
}
int colvarbias_restraint_centers_moving::update_acc_work()
{
if (b_output_acc_work) {
- if ((cvm::step_relative() > 0) || (cvm::step_absolute() == 0)) {
+ if ((cvm::step_relative() > 0) &&
+ (cvm::step_absolute() <= target_nsteps)) {
for (size_t i = 0; i < num_variables(); i++) {
// project forces on the calculated increments at this step
acc_work += colvar_forces[i] * centers_incr[i];
@@ -383,13 +403,6 @@ std::string const colvarbias_restraint_centers_moving::get_state_params() const
<< colvar_centers[i];
}
os << "\n";
- os << "centers_raw ";
- for (i = 0; i < num_variables(); i++) {
- os << " "
- << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
- << colvar_centers_raw[i];
- }
- os << "\n";
if (b_output_acc_work) {
os << "accumulatedWork "
@@ -398,7 +411,7 @@ std::string const colvarbias_restraint_centers_moving::get_state_params() const
}
}
- return colvarbias_restraint_moving::get_state_params() + os.str();
+ return os.str();
}
@@ -410,8 +423,6 @@ int colvarbias_restraint_centers_moving::set_state_params(std::string const &con
// cvm::log ("Reading the updated restraint centers from the restart.\n");
if (!get_keyval(conf, "centers", colvar_centers))
cvm::error("Error: restraint centers are missing from the restart.\n");
- if (!get_keyval(conf, "centers_raw", colvar_centers_raw))
- cvm::error("Error: \"raw\" restraint centers are missing from the restart.\n");
if (b_output_acc_work) {
if (!get_keyval(conf, "accumulatedWork", acc_work))
cvm::error("Error: accumulatedWork is missing from the restart.\n");
@@ -609,7 +620,7 @@ std::string const colvarbias_restraint_k_moving::get_state_params() const
<< std::setprecision(cvm::en_prec)
<< std::setw(cvm::en_width) << force_k << "\n";
}
- return colvarbias_restraint_moving::get_state_params() + os.str();
+ return os.str();
}
@@ -770,6 +781,7 @@ cvm::real colvarbias_restraint_harmonic::d_restraint_potential_dk(size_t i) cons
std::string const colvarbias_restraint_harmonic::get_state_params() const
{
return colvarbias_restraint::get_state_params() +
+ colvarbias_restraint_moving::get_state_params() +
colvarbias_restraint_centers_moving::get_state_params() +
colvarbias_restraint_k_moving::get_state_params();
}
@@ -779,6 +791,7 @@ int colvarbias_restraint_harmonic::set_state_params(std::string const &conf)
{
int error_code = COLVARS_OK;
error_code |= colvarbias_restraint::set_state_params(conf);
+ error_code |= colvarbias_restraint_moving::set_state_params(conf);
error_code |= colvarbias_restraint_centers_moving::set_state_params(conf);
error_code |= colvarbias_restraint_k_moving::set_state_params(conf);
return error_code;
@@ -1037,6 +1050,7 @@ cvm::real colvarbias_restraint_harmonic_walls::d_restraint_potential_dk(size_t i
std::string const colvarbias_restraint_harmonic_walls::get_state_params() const
{
return colvarbias_restraint::get_state_params() +
+ colvarbias_restraint_moving::get_state_params() +
colvarbias_restraint_k_moving::get_state_params();
}
@@ -1045,6 +1059,7 @@ int colvarbias_restraint_harmonic_walls::set_state_params(std::string const &con
{
int error_code = COLVARS_OK;
error_code |= colvarbias_restraint::set_state_params(conf);
+ error_code |= colvarbias_restraint_moving::set_state_params(conf);
error_code |= colvarbias_restraint_k_moving::set_state_params(conf);
return error_code;
}
@@ -1164,6 +1179,7 @@ cvm::real colvarbias_restraint_linear::d_restraint_potential_dk(size_t i) const
std::string const colvarbias_restraint_linear::get_state_params() const
{
return colvarbias_restraint::get_state_params() +
+ colvarbias_restraint_moving::get_state_params() +
colvarbias_restraint_centers_moving::get_state_params() +
colvarbias_restraint_k_moving::get_state_params();
}
@@ -1173,6 +1189,7 @@ int colvarbias_restraint_linear::set_state_params(std::string const &conf)
{
int error_code = COLVARS_OK;
error_code |= colvarbias_restraint::set_state_params(conf);
+ error_code |= colvarbias_restraint_moving::set_state_params(conf);
error_code |= colvarbias_restraint_centers_moving::set_state_params(conf);
error_code |= colvarbias_restraint_k_moving::set_state_params(conf);
return error_code;
diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h
index 98b967abdb1ba6a50754e2f26c3d9600e742557a..8c3a1537fc881cbca3dc8ef66fcf3caeb4d49c90 100644
--- a/lib/colvars/colvarbias_restraint.h
+++ b/lib/colvars/colvarbias_restraint.h
@@ -74,9 +74,6 @@ protected:
/// \brief Restraint centers
std::vector<colvarvalue> colvar_centers;
-
- /// \brief Restraint centers outside the domain of the colvars (no wrapping or constraints applied)
- std::vector<colvarvalue> colvar_centers_raw;
};
@@ -156,10 +153,16 @@ protected:
/// \brief New restraint centers
std::vector<colvarvalue> target_centers;
+ /// \brief Initial value of the restraint centers
+ std::vector<colvarvalue> initial_centers;
+
/// \brief Amplitude of the restraint centers' increment at each step
- /// (or stage) towards the new values (calculated from target_nsteps)
+ /// towards the new values (calculated from target_nsteps)
std::vector<colvarvalue> centers_incr;
+ /// \brief Update the centers by interpolating between initial and target
+ virtual int update_centers(cvm::real lambda);
+
/// Whether to write the current restraint centers to the trajectory file
bool b_output_centers;
diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h
index 2c865a166b594a8f45876c2aa3f426c3196c67a2..3c1ec2495c9c5af29c5549b9c9ddaa355d2f1502 100644
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@@ -132,9 +132,15 @@ public:
static std::vector<feature *> cvc_features;
/// \brief Implementation of the feature list accessor for colvar
- virtual std::vector<feature *> &features() {
+ virtual const std::vector<feature *> &features()
+ {
return cvc_features;
}
+ virtual std::vector<feature *> &modify_features()
+ {
+ return cvc_features;
+ }
+
/// \brief Obtain data needed for the calculation for the backend
virtual void read_data();
diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp
index 5402836f53404b975fee58a0d4c88d00684a1a30..8f241a6255953f2ae0d4815a050daa12c077f866 100644
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@@ -374,8 +374,8 @@ int colvardeps::decr_ref_count(int feature_id) {
}
void colvardeps::init_feature(int feature_id, const char *description, feature_type type) {
- features()[feature_id]->description = description;
- features()[feature_id]->type = type;
+ modify_features()[feature_id]->description = description;
+ modify_features()[feature_id]->type = type;
}
// Shorthand macros for describing dependencies
@@ -401,7 +401,7 @@ void colvardeps::init_cvb_requires() {
int i;
if (features().size() == 0) {
for (i = 0; i < f_cvb_ntot; i++) {
- features().push_back(new feature);
+ modify_features().push_back(new feature);
}
init_feature(f_cvb_active, "active", f_type_dynamic);
@@ -438,7 +438,7 @@ void colvardeps::init_cv_requires() {
size_t i;
if (features().size() == 0) {
for (i = 0; i < f_cv_ntot; i++) {
- features().push_back(new feature);
+ modify_features().push_back(new feature);
}
init_feature(f_cv_active, "active", f_type_dynamic);
@@ -554,7 +554,7 @@ void colvardeps::init_cvc_requires() {
// Initialize static array once and for all
if (features().size() == 0) {
for (i = 0; i < colvardeps::f_cvc_ntot; i++) {
- features().push_back(new feature);
+ modify_features().push_back(new feature);
}
init_feature(f_cvc_active, "active", f_type_dynamic);
@@ -633,7 +633,7 @@ void colvardeps::init_ag_requires() {
// Initialize static array once and for all
if (features().size() == 0) {
for (i = 0; i < f_ag_ntot; i++) {
- features().push_back(new feature);
+ modify_features().push_back(new feature);
}
init_feature(f_ag_active, "active", f_type_dynamic);
diff --git a/lib/colvars/colvardeps.h b/lib/colvars/colvardeps.h
index b810a5fca1f491c415b52d941c4c0848eb2f3ffd..dfb10d00e421f7635fe4bce5ec88b5f45cfc39eb 100644
--- a/lib/colvars/colvardeps.h
+++ b/lib/colvars/colvardeps.h
@@ -135,7 +135,8 @@ public:
// with a non-static array
// Intermediate classes (colvarbias and colvarcomp, which are also base classes)
// implement this as virtual to allow overriding
- virtual std::vector<feature *>&features() = 0;
+ virtual const std::vector<feature *>&features() = 0;
+ virtual std::vector<feature *>&modify_features() = 0;
void add_child(colvardeps *child);
diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h
index e544756428c79abb131930841c473e6ac3553340..312c0fd1a0c86e832ccbac15605867ddd101c824 100644
--- a/lib/colvars/colvars_version.h
+++ b/lib/colvars/colvars_version.h
@@ -1,4 +1,5 @@
-#define COLVARS_VERSION "2017-07-15"
+#ifndef COLVARS_VERSION
+#define COLVARS_VERSION "2017-08-06"
// 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
@@ -6,3 +7,4 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
+#endif
diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp
index 5bb2faae24ad94dc385470ff521958c6ab5b68ce..89302a16a2ead1c29bbc501ec6d8e2aa9d043cb8 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -472,7 +472,7 @@ int colvarscript::proc_features(colvardeps *obj,
}
if ((subcmd == "get") || (subcmd == "set")) {
- std::vector<colvardeps::feature *> &features = obj->features();
+ std::vector<colvardeps::feature *> const &features = obj->features();
std::string const req_feature(obj_to_str(objv[3]));
colvardeps::feature *f = NULL;
int fid = 0;
diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp
index 5200d4d0415838c52c4dd34e101b619312e2c43f..428fe1a4b1d864c725cd7d5b09652216a03091a3 100644
--- a/lib/colvars/colvartypes.cpp
+++ b/lib/colvars/colvartypes.cpp
@@ -19,6 +19,17 @@ bool colvarmodule::rotation::monitor_crossings = false;
cvm::real colvarmodule::rotation::crossing_threshold = 1.0E-02;
+/// Numerical recipes diagonalization
+static int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot);
+
+/// Eigenvector sort
+static int eigsrt(cvm::real *d, cvm::real **v);
+
+/// Transpose the matrix
+static int transpose(cvm::real **v);
+
+
+
std::string cvm::rvector::to_simple_string() const
{
std::ostringstream os;
@@ -286,7 +297,12 @@ void colvarmodule::rotation::diagonalize_matrix(cvm::matrix2d<cvm::real> &S,
// diagonalize
int jac_nrot = 0;
- jacobi(S.c_array(), S_eigval.c_array(), S_eigvec.c_array(), &jac_nrot);
+ if (jacobi(S.c_array(), S_eigval.c_array(), S_eigvec.c_array(), &jac_nrot) !=
+ COLVARS_OK) {
+ cvm::error("Too many iterations in routine jacobi.\n"
+ "This is usually the result of an ill-defined set of atoms for "
+ "rotational alignment (RMSD, rotateReference, etc).\n");
+ }
eigsrt(S_eigval.c_array(), S_eigvec.c_array());
// jacobi saves eigenvectors by columns
transpose(S_eigvec.c_array());
@@ -528,7 +544,7 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector<cvm::atom_pos> co
#define n 4
-void jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot)
+int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot)
{
int j,iq,ip,i;
cvm::real tresh,theta,tau,t,sm,s,h,g,c;
@@ -554,7 +570,7 @@ void jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot)
sm += std::fabs(a[ip][iq]);
}
if (sm == 0.0) {
- return;
+ return COLVARS_OK;
}
if (i < 4)
tresh=0.2*sm/(n*n);
@@ -606,10 +622,11 @@ void jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot)
z[ip]=0.0;
}
}
- cvm::error("Too many iterations in routine jacobi.\n");
+ return COLVARS_ERROR;
}
-void eigsrt(cvm::real *d, cvm::real **v)
+
+int eigsrt(cvm::real *d, cvm::real **v)
{
int k,j,i;
cvm::real p;
@@ -628,9 +645,11 @@ void eigsrt(cvm::real *d, cvm::real **v)
}
}
}
+ return COLVARS_OK;
}
-void transpose(cvm::real **v)
+
+int transpose(cvm::real **v)
{
cvm::real p;
int i,j;
@@ -641,6 +660,7 @@ void transpose(cvm::real **v)
v[j][i]=p;
}
}
+ return COLVARS_OK;
}
#undef n
diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h
index 17c09a509508a428116b8498703086cd5c09bae9..fe3160eb4b3238ef036c2ea56d6b3954d73cb121 100644
--- a/lib/colvars/colvartypes.h
+++ b/lib/colvars/colvartypes.h
@@ -1020,16 +1020,6 @@ inline cvm::rvector operator * (cvm::rmatrix const &m,
}
-/// Numerical recipes diagonalization
-void jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot);
-
-/// Eigenvector sort
-void eigsrt(cvm::real *d, cvm::real **v);
-
-/// Transpose the matrix
-void transpose(cvm::real **v);
-
-
/// \brief 1-dimensional vector of real numbers with four components and
diff --git a/lib/colvars/colvarvalue.cpp b/lib/colvars/colvarvalue.cpp
index 7b498be6d6a39e4a65069e41a8571084fd55e0f8..312d1016035c3ecb3c28101998d62a01d135af8a 100644
--- a/lib/colvars/colvarvalue.cpp
+++ b/lib/colvars/colvarvalue.cpp
@@ -570,6 +570,50 @@ colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const
}
+/// Return the midpoint between x1 and x2, optionally weighted by lambda
+/// (which must be between 0.0 and 1.0)
+colvarvalue const colvarvalue::interpolate(colvarvalue const &x1,
+ colvarvalue const &x2,
+ cvm::real const lambda)
+{
+ colvarvalue::check_types(x1, x2);
+
+ if ((lambda < 0.0) || (lambda > 1.0)) {
+ cvm::error("Error: trying to interpolate between two colvarvalues with a "
+ "lamdba outside [0:1].\n", BUG_ERROR);
+ }
+
+ colvarvalue interp = ((1.0-lambda)*x1 + lambda*x2);
+ cvm::real const d2 = x1.dist2(x2);
+
+ switch (x1.type()) {
+ case colvarvalue::type_scalar:
+ case colvarvalue::type_3vector:
+ case colvarvalue::type_vector:
+ case colvarvalue::type_unit3vectorderiv:
+ case colvarvalue::type_quaternionderiv:
+ return interp;
+ break;
+ case colvarvalue::type_unit3vector:
+ case colvarvalue::type_quaternion:
+ if (interp.norm()/std::sqrt(d2) < 1.0e-6) {
+ cvm::error("Error: interpolation between "+cvm::to_str(x1)+" and "+
+ cvm::to_str(x2)+" with lambda = "+cvm::to_str(lambda)+
+ " is undefined: result = "+cvm::to_str(interp)+"\n",
+ INPUT_ERROR);
+ }
+ interp.apply_constraints();
+ return interp;
+ break;
+ case colvarvalue::type_notset:
+ default:
+ x1.undef_op();
+ break;
+ }
+ return colvarvalue(colvarvalue::type_notset);
+}
+
+
std::string colvarvalue::to_simple_string() const
{
switch (type()) {
diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h
index fce0e1a970a7cf150439f85be501440ac46c9983..41759e92b02a0f497d7c78147b688d824a5783ef 100644
--- a/lib/colvars/colvarvalue.h
+++ b/lib/colvars/colvarvalue.h
@@ -193,6 +193,12 @@ public:
/// Derivative with respect to this \link colvarvalue \endlink of the square distance
colvarvalue dist2_grad(colvarvalue const &x2) const;
+ /// Return the midpoint between x1 and x2, optionally weighted by lambda
+ /// (which must be between 0.0 and 1.0)
+ static colvarvalue const interpolate(colvarvalue const &x1,
+ colvarvalue const &x2,
+ cvm::real const lambda = 0.5);
+
/// Assignment operator (type of x is checked)
colvarvalue & operator = (colvarvalue const &x);
@@ -285,10 +291,10 @@ public:
cvm::real & operator [] (int const i);
/// Ensure that the two types are the same within a binary operator
- int static check_types(colvarvalue const &x1, colvarvalue const &x2);
+ static int check_types(colvarvalue const &x1, colvarvalue const &x2);
/// Ensure that the two types are the same within an assignment, or that the left side is type_notset
- int static check_types_assign(Type const &vt1, Type const &vt2);
+ static int check_types_assign(Type const &vt1, Type const &vt2);
/// Undefined operation
void undef_op() const;
@@ -317,14 +323,14 @@ public:
/// \brief Optimized routine for the inner product of one collective
/// variable with an array
- void static inner_opt(colvarvalue const &x,
+ static void inner_opt(colvarvalue const &x,
std::vector<colvarvalue>::iterator &xv,
std::vector<colvarvalue>::iterator const &xv_end,
std::vector<cvm::real>::iterator &result);
/// \brief Optimized routine for the inner product of one collective
/// variable with an array
- void static inner_opt(colvarvalue const &x,
+ static void inner_opt(colvarvalue const &x,
std::list<colvarvalue>::iterator &xv,
std::list<colvarvalue>::iterator const &xv_end,
std::vector<cvm::real>::iterator &result);
@@ -332,14 +338,14 @@ public:
/// \brief Optimized routine for the second order Legendre
/// polynomial, (3cos^2(w)-1)/2, of one collective variable with an
/// array
- void static p2leg_opt(colvarvalue const &x,
+ static void p2leg_opt(colvarvalue const &x,
std::vector<colvarvalue>::iterator &xv,
std::vector<colvarvalue>::iterator const &xv_end,
std::vector<cvm::real>::iterator &result);
/// \brief Optimized routine for the second order Legendre
/// polynomial of one collective variable with an array
- void static p2leg_opt(colvarvalue const &x,
+ static void p2leg_opt(colvarvalue const &x,
std::list<colvarvalue>::iterator &xv,
std::list<colvarvalue>::iterator const &xv_end,
std::vector<cvm::real>::iterator &result);
diff --git a/src/USER-COLVARS/colvarproxy_lammps_version.h b/src/USER-COLVARS/colvarproxy_lammps_version.h
index 834bd1748a60ee7aa05b9ae0c1abff6979a16cab..0eb6f2d95ac6bbb4c98fe5c75fc2308a7758342c 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-07-15"
+#define COLVARPROXY_VERSION "2017-07-19"
// 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