jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Update Colvars to version 2019-08-01 

One new feature (geometric path variables) and bugfixes.

Minimize occurrence of file-read errors in multiple-walker metadynamics
https://github.com/Colvars/colvars/pull/276

Better error messages for RMSD
https://github.com/Colvars/colvars/pull/272

Ensemble-biased metadynamics documentation
https://github.com/Colvars/colvars/pull/261

Fix bug on multiple walkers metadynamics, pmf included twice
https://github.com/Colvars/colvars/pull/259

Implementation of the geometric path collective variables (C++11-only feature)
https://github.com/Colvars/colvars/pull/249
This commit is contained in:
Giacomo Fiorin 2019-08-01 15:12:42 -04:00
parent 7f342b1cd0
commit 52e2db44a1
21 changed files with 1687 additions and 136 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/src/PDF/colvars-refman-lammps.pdf
View File

Binary file not shown.

1
lib/colvars/Makefile.common
View File

@ -34,6 +34,7 @@ COLVARS_SRCS = \
colvarcomp_coordnums.cpp \
colvarcomp.cpp \
colvarcomp_distances.cpp \
colvarcomp_gpath.cpp \
colvarcomp_protein.cpp \
colvarcomp_rotations.cpp \
colvar.cpp \

20
lib/colvars/colvar.cpp
View File

@ -791,6 +791,11 @@ int colvar::init_components(std::string const &conf)
"inertia", "inertia");
error_code |= init_components_type<inertia_z>(conf, "moment of inertia around an axis", "inertiaZ");
error_code |= init_components_type<eigenvector>(conf, "eigenvector", "eigenvector");
error_code |= init_components_type<gspath>(conf, "geometrical path collective variables (s)", "gspath");
error_code |= init_components_type<gzpath>(conf, "geometrical path collective variables (z)", "gzpath");
error_code |= init_components_type<linearCombination>(conf, "linear combination of other collective variables", "subColvar");
error_code |= init_components_type<gspathCV>(conf, "geometrical path collective variables (s) for other CVs", "gspathCV");
error_code |= init_components_type<gzpathCV>(conf, "geometrical path collective variables (z) for other CVs", "gzpathCV");
if (!cvcs.size() || (error_code != COLVARS_OK)) {
cvm::error("Error: no valid components were provided "
@ -1495,6 +1500,8 @@ int colvar::calc_colvar_properties()
// calculate the velocity by finite differences
if (cvm::step_relative() == 0) {
x_old = x;
v_fdiff.reset(); // Do not pretend we know anything about the actual velocity
// eg. upon restarting. That would require saving v_fdiff or x_old to the state file
} else {
v_fdiff = fdiff_velocity(x_old, x);
v_reported = v_fdiff;
@ -1516,8 +1523,9 @@ int colvar::calc_colvar_properties()
x_ext = prev_x_ext;
v_ext = prev_v_ext;
}
// report the restraint center as "value"
// These position and velocities come from integration at the _previous timestep_ in update_forces_energy()
// But we report values at the beginning of the timestep (value at t=0 on the first timestep)
x_reported = x_ext;
v_reported = v_ext;
// the "total force" with the extended Lagrangian is
@ -1651,8 +1659,6 @@ cvm::real colvar::update_forces_energy()
// is equal to the actual coordinate
x_ext = x;
}
// Report extended value
x_reported = x_ext;
}
// Now adding the force on the actual colvar (for those biases that
@ -1975,9 +1981,8 @@ std::istream & colvar::read_restart(std::istream &is)
}
if (is_enabled(f_cv_extended_Lagrangian)) {
if ( !(get_keyval(conf, "extended_x", x_ext,
colvarvalue(x.type()), colvarparse::parse_silent)) &&
colvarvalue(x.type()), colvarparse::parse_silent)) ||
!(get_keyval(conf, "extended_v", v_ext,
colvarvalue(x.type()), colvarparse::parse_silent)) ) {
cvm::log("Error: restart file does not contain "
@ -2079,11 +2084,11 @@ std::ostream & colvar::write_restart(std::ostream &os) {
os << " extended_x "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< x_ext << "\n"
<< x_reported << "\n"
<< " extended_v "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< v_ext << "\n";
<< v_reported << "\n";
}
os << "}\n\n";
@ -2150,7 +2155,6 @@ std::ostream & colvar::write_traj_label(std::ostream & os)
std::ostream & colvar::write_traj(std::ostream &os)
{
os << " ";
if (is_enabled(f_cv_output_value)) {
if (is_enabled(f_cv_extended_Lagrangian)) {

8
lib/colvars/colvar.h
View File

@ -569,6 +569,14 @@ public:
class alpha_dihedrals;
class alpha_angles;
class dihedPC;
class componentDisabled;
class CartesianBasedPath;
class gspath;
class gzpath;
class linearCombination;
class CVBasedPath;
class gspathCV;
class gzpathCV;
// non-scalar components
class distance_vec;

270
lib/colvars/colvar_geometricpath.h Normal file
View File

@ -0,0 +1,270 @@
#ifndef GEOMETRICPATHCV_H
#define GEOMETRICPATHCV_H
// 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
// Please update all Colvars source files before making any changes.
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include <vector>
#include <cmath>
#include <iostream>
#include <algorithm>
#include <string>
#include <map>
namespace GeometricPathCV {
enum path_sz {S, Z};
template <typename element_type, typename scalar_type, path_sz path_type>
class GeometricPathBase {
private:
struct doCompareFrameDistance {
doCompareFrameDistance(const GeometricPathBase& obj): m_obj(obj) {}
const GeometricPathBase& m_obj;
bool operator()(const size_t& i1, const size_t& i2) {
return m_obj.frame_distances[i1] < m_obj.frame_distances[i2];
}
};
protected:
scalar_type v1v1;
scalar_type v2v2;
scalar_type v3v3;
scalar_type v4v4;
scalar_type v1v3;
scalar_type v1v4;
scalar_type f;
scalar_type dx;
scalar_type s;
scalar_type z;
scalar_type zz;
std::vector<element_type> v1;
std::vector<element_type> v2;
std::vector<element_type> v3;
std::vector<element_type> v4;
std::vector<element_type> dfdv1;
std::vector<element_type> dfdv2;
std::vector<element_type> dzdv1;
std::vector<element_type> dzdv2;
std::vector<scalar_type> frame_distances;
std::vector<size_t> frame_index;
bool use_second_closest_frame;
bool use_third_closest_frame;
bool use_z_square;
long min_frame_index_1;
long min_frame_index_2;
long min_frame_index_3;
long sign;
double M;
double m;
public:
GeometricPathBase(size_t vector_size, const element_type& element = element_type(), size_t total_frames = 1, bool p_use_second_closest_frame = true, bool p_use_third_closest_frame = false, bool p_use_z_square = false);
GeometricPathBase(size_t vector_size, const std::vector<element_type>& elements, size_t total_frames = 1, bool p_use_second_closest_frame = true, bool p_use_third_closest_frame = false, bool p_use_z_square = false);
GeometricPathBase() {}
virtual ~GeometricPathBase() {}
virtual void initialize(size_t vector_size, const element_type& element = element_type(), size_t total_frames = 1, bool p_use_second_closest_frame = true, bool p_use_third_closest_frame = false, bool p_use_z_square = false);
virtual void initialize(size_t vector_size, const std::vector<element_type>& elements, size_t total_frames = 1, bool p_use_second_closest_frame = true, bool p_use_third_closest_frame = false, bool p_use_z_square = false);
virtual void prepareVectors();
virtual void updateReferenceDistances();
virtual void compute();
virtual void determineClosestFrames();
virtual void computeValue();
virtual void computeDerivatives();
};
template <typename element_type, typename scalar_type, path_sz path_type>
GeometricPathBase<element_type, scalar_type, path_type>::GeometricPathBase(size_t vector_size, const element_type& element, size_t total_frames, bool p_use_second_closest_frame, bool p_use_third_closest_frame, bool p_use_z_square) {
initialize(vector_size, element, total_frames, p_use_second_closest_frame, p_use_third_closest_frame, p_use_z_square);
}
template <typename element_type, typename scalar_type, path_sz path_type>
GeometricPathBase<element_type, scalar_type, path_type>::GeometricPathBase(size_t vector_size, const std::vector<element_type>& elements, size_t total_frames, bool p_use_second_closest_frame, bool p_use_third_closest_frame, bool p_use_z_square) {
initialize(vector_size, elements, total_frames, p_use_second_closest_frame, p_use_third_closest_frame, p_use_z_square);
}
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::initialize(size_t vector_size, const element_type& element, size_t total_frames, bool p_use_second_closest_frame, bool p_use_third_closest_frame, bool p_use_z_square) {
v1v1 = scalar_type();
v2v2 = scalar_type();
v3v3 = scalar_type();
v4v4 = scalar_type();
v1v3 = scalar_type();
v1v4 = scalar_type();
f = scalar_type();
dx = scalar_type();
z = scalar_type();
zz = scalar_type();
sign = 0;
v1.resize(vector_size, element);
v2.resize(vector_size, element);
v3.resize(vector_size, element);
v4.resize(vector_size, element);
dfdv1.resize(vector_size, element);
dfdv2.resize(vector_size, element);
dzdv1.resize(vector_size, element);
dzdv2.resize(vector_size, element);
frame_distances.resize(total_frames);
frame_index.resize(total_frames);
for (size_t i_frame = 0; i_frame < frame_index.size(); ++i_frame) {
frame_index[i_frame] = i_frame;
}
use_second_closest_frame = p_use_second_closest_frame;
use_third_closest_frame = p_use_third_closest_frame;
use_z_square = p_use_z_square;
M = static_cast<scalar_type>(total_frames - 1);
m = static_cast<scalar_type>(1.0);
}
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::initialize(size_t vector_size, const std::vector<element_type>& elements, size_t total_frames, bool p_use_second_closest_frame, bool p_use_third_closest_frame, bool p_use_z_square) {
v1v1 = scalar_type();
v2v2 = scalar_type();
v3v3 = scalar_type();
v4v4 = scalar_type();
v1v3 = scalar_type();
v1v4 = scalar_type();
f = scalar_type();
dx = scalar_type();
z = scalar_type();
zz = scalar_type();
sign = 0;
v1 = elements;
v2 = elements;
v3 = elements;
v4 = elements;
dfdv1 = elements;
dfdv2 = elements;
dzdv1 = elements;
dzdv2 = elements;
frame_distances.resize(total_frames);
frame_index.resize(total_frames);
for (size_t i_frame = 0; i_frame < frame_index.size(); ++i_frame) {
frame_index[i_frame] = i_frame;
}
use_second_closest_frame = p_use_second_closest_frame;
use_third_closest_frame = p_use_third_closest_frame;
use_z_square = p_use_z_square;
M = static_cast<scalar_type>(total_frames - 1);
m = static_cast<scalar_type>(1.0);
}
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::prepareVectors() {
std::cout << "Warning: you should not call the prepareVectors() in base class!\n";
std::cout << std::flush;
}
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::updateReferenceDistances() {
std::cout << "Warning: you should not call the updateReferenceDistances() in base class!\n";
std::cout << std::flush;
}
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::compute() {
computeValue();
computeDerivatives();
}
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::determineClosestFrames() {
// Find the closest and the second closest frames
std::sort(frame_index.begin(), frame_index.end(), doCompareFrameDistance(*this));
// Determine the sign
sign = static_cast<long>(frame_index[0]) - static_cast<long>(frame_index[1]);
if (sign > 1) {
// sigma(z) is on the left side of the closest frame
sign = 1;
} else if (sign < -1) {
// sigma(z) is on the right side of the closest frame
sign = -1;
}
if (std::abs(static_cast<long>(frame_index[0]) - static_cast<long>(frame_index[1])) > 1) {
std::cout << "Warning: Geometrical pathCV relies on the assumption that the second closest frame is the neighbouring frame\n";
std::cout << " Please check your configuration or increase restraint on z(σ)\n";
for (size_t i_frame = 0; i_frame < frame_index.size(); ++i_frame) {
std::cout << "Frame index: " << frame_index[i_frame] << " ; optimal RMSD = " << frame_distances[frame_index[i_frame]] << "\n";
}
}
min_frame_index_1 = frame_index[0]; // s_m
min_frame_index_2 = use_second_closest_frame ? frame_index[1] : min_frame_index_1 - sign; // s_(m-1)
min_frame_index_3 = use_third_closest_frame ? frame_index[2] : min_frame_index_1 + sign; // s_(m+1)
m = static_cast<double>(frame_index[0]);
}
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::computeValue() {
updateReferenceDistances();
determineClosestFrames();
prepareVectors();
v1v1 = scalar_type();
v2v2 = scalar_type();
v3v3 = scalar_type();
v1v3 = scalar_type();
if (path_type == Z) {
v1v4 = scalar_type();
v4v4 = scalar_type();
}
for (size_t i_elem = 0; i_elem < v1.size(); ++i_elem) {
v1v1 += v1[i_elem] * v1[i_elem];
v2v2 += v2[i_elem] * v2[i_elem];
v3v3 += v3[i_elem] * v3[i_elem];
v1v3 += v1[i_elem] * v3[i_elem];
if (path_type == Z) {
v1v4 += v1[i_elem] * v4[i_elem];
v4v4 += v4[i_elem] * v4[i_elem];
}
}
f = (std::sqrt(v1v3 * v1v3 - v3v3 * (v1v1 - v2v2)) - v1v3) / v3v3;
if (path_type == Z) {
dx = 0.5 * (f - 1);
zz = v1v1 + 2 * dx * v1v4 + dx * dx * v4v4;
if (use_z_square) {
z = zz;
} else {
z = std::sqrt(std::fabs(zz));
}
}
if (path_type == S) {
s = m/M + static_cast<double>(sign) * ((f - 1) / (2 * M));
}
}
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::computeDerivatives() {
const scalar_type factor1 = 1.0 / (2.0 * v3v3 * std::sqrt(v1v3 * v1v3 - v3v3 * (v1v1 - v2v2)));
const scalar_type factor2 = 1.0 / v3v3;
for (size_t i_elem = 0; i_elem < v1.size(); ++i_elem) {
// Compute the derivative of f with vector v1
dfdv1[i_elem] = factor1 * (2.0 * v1v3 * v3[i_elem] - 2.0 * v3v3 * v1[i_elem]) - factor2 * v3[i_elem];
// Compute the derivative of f with respect to vector v2
dfdv2[i_elem] = factor1 * (2.0 * v3v3 * v2[i_elem]);
// dZ(v1(r), v2(r), v3) / dr = ∂Z/∂v1 * dv1/dr + ∂Z/∂v2 * dv2/dr
// dv1/dr = [fitting matrix 1][-1, ..., -1]
// dv2/dr = [fitting matrix 2][1, ..., 1]
// ∂Z/∂v1 = 1/(2*z) * (2v1 + (f-1)v4 + (v1⋅v4)∂f/∂v1 + v4^2 * 1/4 * 2(f-1) * ∂f/∂v1)
// ∂Z/∂v2 = 1/(2*z) * ((v1⋅v4)∂f/∂v2 + v4^2 * 1/4 * 2(f-1) * ∂f/∂v2)
if (path_type == Z) {
if (use_z_square) {
dzdv1[i_elem] = 2.0 * v1[i_elem] + (f-1) * v4[i_elem] + v1v4 * dfdv1[i_elem] + v4v4 * 0.25 * 2.0 * (f-1) * dfdv1[i_elem];
dzdv2[i_elem] = v1v4 * dfdv2[i_elem] + v4v4 * 0.25 * 2.0 * (f-1) * dfdv2[i_elem];
} else {
if (z > static_cast<scalar_type>(0)) {
dzdv1[i_elem] = (1.0 / (2.0 * z)) * (2.0 * v1[i_elem] + (f-1) * v4[i_elem] + v1v4 * dfdv1[i_elem] + v4v4 * 0.25 * 2.0 * (f-1) * dfdv1[i_elem]);
dzdv2[i_elem] = (1.0 / (2.0 * z)) * (v1v4 * dfdv2[i_elem] + v4v4 * 0.25 * 2.0 * (f-1) * dfdv2[i_elem]);
} else {
// workaround at z = 0
dzdv1[i_elem] = 0;
dzdv2[i_elem] = 0;
}
}
}
}
}
}
#endif // GEOMETRICPATHCV_H

202
lib/colvars/colvarbias_meta.cpp
View File

@ -37,7 +37,6 @@ colvarbias_meta::colvarbias_meta(char const *key)
{
new_hills_begin = hills.end();
hills_traj_os = NULL;
replica_hills_os = NULL;
ebmeta_equil_steps = 0L;
}
@ -203,23 +202,37 @@ int colvarbias_meta::init_ebmeta_params(std::string const &conf)
target_dist = new colvar_grid_scalar();
target_dist->init_from_colvars(colvars);
std::string target_dist_file;
get_keyval(conf, "targetdistfile", target_dist_file);
get_keyval(conf, "targetDistFile", target_dist_file);
std::ifstream targetdiststream(target_dist_file.c_str());
target_dist->read_multicol(targetdiststream);
cvm::real min_val = target_dist->minimum_value();
cvm::real max_val = target_dist->maximum_value();
if(min_val<0){
cvm::error("Error: Target distribution of ebMeta "
cvm::error("Error: Target distribution of EBMetaD "
"has negative values!.\n", INPUT_ERROR);
}
cvm::real min_pos_val = target_dist->minimum_pos_value();
if(min_pos_val<=0){
cvm::error("Error: Target distribution of ebMeta has negative "
"or zero minimum positive value!.\n", INPUT_ERROR);
}
if(min_val==0){
cvm::log("WARNING: Target distribution has zero values.\n");
cvm::log("Zeros will be converted to the minimum positive value.\n");
target_dist->remove_zeros(min_pos_val);
cvm::real target_dist_min_val;
get_keyval(conf, "targetDistMinVal", target_dist_min_val, 1/1000000.0);
if(target_dist_min_val>0 && target_dist_min_val<1){
target_dist_min_val=max_val*target_dist_min_val;
target_dist->remove_small_values(target_dist_min_val);
} else {
if (target_dist_min_val==0) {
cvm::log("NOTE: targetDistMinVal is set to zero, the minimum value of the target \n");
cvm::log(" distribution will be set as the minimum positive value.\n");
cvm::real min_pos_val = target_dist->minimum_pos_value();
if(min_pos_val<=0){
cvm::error("Error: Target distribution of EBMetaD has negative "
"or zero minimum positive value!.\n", INPUT_ERROR);
}
if(min_val==0){
cvm::log("WARNING: Target distribution has zero values.\n");
cvm::log("Zeros will be converted to the minimum positive value.\n");
target_dist->remove_small_values(min_pos_val);
}
} else {
cvm::error("Error: targetDistMinVal must be a value between 0 and 1!.\n", INPUT_ERROR);
}
}
// normalize target distribution and multiply by effective volume = exp(differential entropy)
target_dist->multiply_constant(1.0/target_dist->integral());
@ -235,14 +248,14 @@ int colvarbias_meta::init_ebmeta_params(std::string const &conf)
colvarbias_meta::~colvarbias_meta()
{
colvarbias_meta::clear_state_data();
colvarproxy *proxy = cvm::proxy;
if (replica_hills_os) {
cvm::proxy->close_output_stream(replica_hills_file);
replica_hills_os = NULL;
if (proxy->get_output_stream(replica_hills_file)) {
proxy->close_output_stream(replica_hills_file);
}
if (hills_traj_os) {
cvm::proxy->close_output_stream(hills_traj_file_name());
proxy->close_output_stream(hills_traj_file_name());
hills_traj_os = NULL;
}
@ -523,6 +536,8 @@ int colvarbias_meta::update_bias()
case multiple_replicas:
create_hill(hill(hill_weight*hills_scale, colvars, hill_width, replica_id));
std::ostream *replica_hills_os =
cvm::proxy->get_output_stream(replica_hills_file);
if (replica_hills_os) {
*replica_hills_os << hills.back();
} else {
@ -921,13 +936,16 @@ void colvarbias_meta::recount_hills_off_grid(colvarbias_meta::hill_iter h_first
int colvarbias_meta::replica_share()
{
colvarproxy *proxy = cvm::proxy;
// sync with the other replicas (if needed)
if (comm == multiple_replicas) {
// reread the replicas registry
update_replicas_registry();
// empty the output buffer
std::ostream *replica_hills_os =
proxy->get_output_stream(replica_hills_file);
if (replica_hills_os) {
cvm::proxy->flush_output_stream(replica_hills_os);
proxy->flush_output_stream(replica_hills_os);
}
read_replica_files();
}
@ -1089,11 +1107,7 @@ void colvarbias_meta::read_replica_files()
(replicas[ir])->replica_state_file+"\".\n");
std::ifstream is((replicas[ir])->replica_state_file.c_str());
if (! (replicas[ir])->read_state(is)) {
cvm::log("Reading from file \""+(replicas[ir])->replica_state_file+
"\" failed or incomplete: will try again in "+
cvm::to_str(replica_update_freq)+" steps.\n");
} else {
if ((replicas[ir])->read_state(is)) {
// state file has been read successfully
(replicas[ir])->replica_state_file_in_sync = true;
(replicas[ir])->update_status = 0;
@ -1550,15 +1564,11 @@ int colvarbias_meta::setup_output()
// for the others to read
// open the "hills" buffer file
if (!replica_hills_os) {
cvm::proxy->backup_file(replica_hills_file);
replica_hills_os = cvm::proxy->output_stream(replica_hills_file);
if (!replica_hills_os) return cvm::get_error();
replica_hills_os->setf(std::ios::scientific, std::ios::floatfield);
}
reopen_replica_buffer_file();
// write the state file (so that there is always one available)
write_replica_state_file();
// schedule to read the state files of the other replicas
for (size_t ir = 0; ir < replicas.size(); ir++) {
(replicas[ir])->replica_state_file_in_sync = false;
@ -1661,15 +1671,16 @@ std::ostream & colvarbias_meta::write_state_data(std::ostream& os)
int colvarbias_meta::write_state_to_replicas()
{
int error_code = COLVARS_OK;
if (comm != single_replica) {
write_replica_state_file();
// schedule to reread the state files of the other replicas (they
// have also rewritten them)
error_code |= write_replica_state_file();
error_code |= reopen_replica_buffer_file();
// schedule to reread the state files of the other replicas
for (size_t ir = 0; ir < replicas.size(); ir++) {
(replicas[ir])->replica_state_file_in_sync = false;
}
}
return COLVARS_OK;
return error_code;
}
@ -1693,6 +1704,20 @@ void colvarbias_meta::write_pmf()
// output the PMF from this instance or replica
pmf->reset();
pmf->add_grid(*hills_energy);
if (ebmeta) {
int nt_points=pmf->number_of_points();
for (int i = 0; i < nt_points; i++) {
cvm:: real pmf_val=0.0;
cvm:: real target_val=target_dist->value(i);
if (target_val>0) {
pmf_val=pmf->value(i);
pmf_val=pmf_val+cvm::temperature() * cvm::boltzmann() * std::log(target_val);
}
pmf->set_value(i,pmf_val);
}
}
cvm::real const max = pmf->maximum_value();
pmf->add_constant(-1.0 * max);
pmf->multiply_constant(-1.0);
@ -1716,10 +1741,24 @@ void colvarbias_meta::write_pmf()
if (comm != single_replica) {
// output the combined PMF from all replicas
pmf->reset();
pmf->add_grid(*hills_energy);
// current replica already included in the pools of replicas
for (size_t ir = 0; ir < replicas.size(); ir++) {
pmf->add_grid(*(replicas[ir]->hills_energy));
}
if (ebmeta) {
int nt_points=pmf->number_of_points();
for (int i = 0; i < nt_points; i++) {
cvm:: real pmf_val=0.0;
cvm:: real target_val=target_dist->value(i);
if (target_val>0) {
pmf_val=pmf->value(i);
pmf_val=pmf_val+cvm::temperature() * cvm::boltzmann() * std::log(target_val);
}
pmf->set_value(i,pmf_val);
}
}
cvm::real const max = pmf->maximum_value();
pmf->add_constant(-1.0 * max);
pmf->multiply_constant(-1.0);
@ -1744,74 +1783,47 @@ void colvarbias_meta::write_pmf()
int colvarbias_meta::write_replica_state_file()
{
colvarproxy *proxy = cvm::proxy;
if (cvm::debug()) {
cvm::log("Writing replica state file for bias \""+name+"\"\n");
}
// write down also the restart for the other replicas
cvm::backup_file(replica_state_file.c_str());
std::ostream *rep_state_os = cvm::proxy->output_stream(replica_state_file);
if (rep_state_os == NULL) {
cvm::error("Error: in opening file \""+
replica_state_file+"\" for writing.\n", FILE_ERROR);
return FILE_ERROR;
int error_code = COLVARS_OK;
// Write to temporary state file
std::string const tmp_state_file(replica_state_file+".tmp");
error_code |= proxy->remove_file(tmp_state_file);
std::ostream *rep_state_os = cvm::proxy->output_stream(tmp_state_file);
if (rep_state_os) {
if (!write_state(*rep_state_os)) {
error_code |= cvm::error("Error: in writing to temporary file \""+
tmp_state_file+"\".\n", FILE_ERROR);
}
}
error_code |= proxy->close_output_stream(tmp_state_file);
rep_state_os->setf(std::ios::scientific, std::ios::floatfield);
error_code |= proxy->rename_file(tmp_state_file, replica_state_file);
if (!write_state(*rep_state_os)) {
cvm::error("Error: in writing to file \""+
replica_state_file+"\".\n", FILE_ERROR);
cvm::proxy->close_output_stream(replica_state_file);
return FILE_ERROR;
return error_code;
}
int colvarbias_meta::reopen_replica_buffer_file()
{
int error_code = COLVARS_OK;
colvarproxy *proxy = cvm::proxy;
if (proxy->get_output_stream(replica_hills_file) != NULL) {
error_code |= proxy->close_output_stream(replica_hills_file);
}
cvm::proxy->close_output_stream(replica_state_file);
// rep_state_os.setf(std::ios::scientific, std::ios::floatfield);
// rep_state_os << "\n"
// << "metadynamics {\n"
// << " configuration {\n"
// << " name " << this->name << "\n"
// << " step " << cvm::step_absolute() << "\n";
// if (this->comm != single_replica) {
// rep_state_os << " replicaID " << this->replica_id << "\n";
// }
// rep_state_os << " }\n\n";
// rep_state_os << " hills_energy\n";
// rep_state_os << std::setprecision(cvm::cv_prec)
// << std::setw(cvm::cv_width);
// hills_energy->write_restart(rep_state_os);
// rep_state_os << " hills_energy_gradients\n";
// rep_state_os << std::setprecision(cvm::cv_prec)
// << std::setw(cvm::cv_width);
// hills_energy_gradients->write_restart(rep_state_os);
// if ( (!use_grids) || keep_hills ) {
// // write all hills currently in memory
// for (std::list<hill>::const_iterator h = this->hills.begin();
// h != this->hills.end();
// h++) {
// rep_state_os << *h;
// }
// } else {
// // write just those that are near the grid boundaries
// for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
// h != this->hills_off_grid.end();
// h++) {
// rep_state_os << *h;
// }
// }
// rep_state_os << "}\n\n";
// rep_state_os.close();
// reopen the hills file
cvm::proxy->close_output_stream(replica_hills_file);
cvm::proxy->backup_file(replica_hills_file);
replica_hills_os = cvm::proxy->output_stream(replica_hills_file);
if (!replica_hills_os) return cvm::get_error();
replica_hills_os->setf(std::ios::scientific, std::ios::floatfield);
return COLVARS_OK;
error_code |= proxy->remove_file(replica_hills_file);
std::ostream *replica_hills_os = proxy->output_stream(replica_hills_file);
if (replica_hills_os) {
replica_hills_os->setf(std::ios::scientific, std::ios::floatfield);
} else {
error_code |= FILE_ERROR;
}
return error_code;
}
@ -1883,5 +1895,3 @@ std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h)
return os;
}

8
lib/colvars/colvarbias_meta.h
View File

@ -215,9 +215,12 @@ protected:
/// \brief Read new data from replicas' files
virtual void read_replica_files();
/// \brief Write data to other replicas
/// Write full state information to be read by other replicas
virtual int write_replica_state_file();
/// Call this after write_replica_state_file()
virtual int reopen_replica_buffer_file();
/// \brief Additional, "mirror" metadynamics biases, to collect info
/// from the other replicas
///
@ -251,9 +254,6 @@ protected:
/// This file becomes empty after replica_state_file is rewritten
std::string replica_hills_file;
/// \brief Output stream corresponding to replica_hills_file
std::ostream *replica_hills_os;
/// Position within replica_hills_file (when reading it)
int replica_hills_file_pos;

3
lib/colvars/colvarbias_restraint.cpp
View File

@ -1295,7 +1295,8 @@ colvarvalue const colvarbias_restraint_linear::restraint_force(size_t i) const
cvm::real colvarbias_restraint_linear::d_restraint_potential_dk(size_t i) const
{
return 1.0 / variables(i)->width * (variables(i)->value() - colvar_centers[i]);
return 1.0 / variables(i)->width * (variables(i)->value() -
colvar_centers[i]).sum();
}

2
lib/colvars/colvarcomp.cpp
View File

@ -7,6 +7,8 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include <algorithm>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvar.h"

226
lib/colvars/colvarcomp.h
View File

@ -24,6 +24,13 @@
#include "colvar.h"
#include "colvaratoms.h"
#if (__cplusplus >= 201103L)
#include "colvar_geometricpath.h"
#include <functional>
#endif // C++11 checking
#include <map>
/// \brief Colvar component (base class for collective variables)
///
@ -688,9 +695,6 @@ protected:
/// Reference coordinates
std::vector<cvm::atom_pos> ref_pos;
/// Geometric center of the reference coordinates
cvm::atom_pos ref_pos_center;
/// Eigenvector (of a normal or essential mode): will always have zero center
std::vector<cvm::rvector> eigenvec;
@ -1383,6 +1387,222 @@ public:
};
class colvar::componentDisabled
: public colvar::cvc
{
public:
componentDisabled(std::string const &conf) {
cvm::error("Error: this component is not enabled in the current build; please see https://colvars.github.io/README-c++11.html");
}
virtual ~componentDisabled() {}
virtual void calc_value() {}
virtual void calc_gradients() {}
virtual void apply_force(colvarvalue const &force) {}
};
#if (__cplusplus >= 201103L)
class colvar::CartesianBasedPath
: public colvar::cvc
{
protected:
virtual void computeReferenceDistance(std::vector<cvm::real>& result);
/// Selected atoms
cvm::atom_group *atoms;
/// Fitting options
bool has_user_defined_fitting;
/// Reference frames
std::vector<std::vector<cvm::atom_pos>> reference_frames;
std::vector<std::vector<cvm::atom_pos>> reference_fitting_frames;
/// Atom groups for RMSD calculation together with reference frames
std::vector<cvm::atom_group*> comp_atoms;
/// Total number of reference frames
size_t total_reference_frames;
public:
CartesianBasedPath(std::string const &conf);
virtual ~CartesianBasedPath();
virtual void calc_value() = 0;
virtual void apply_force(colvarvalue const &force) = 0;
};
/// \brief Colvar component: alternative path collective variable using geometry, variable s
/// For more information see https://plumed.github.io/doc-v2.5/user-doc/html/_p_a_t_h.html
/// Díaz Leines, G.; Ensing, B. Path Finding on High-Dimensional Free Energy Landscapes. Phys. Rev. Lett. 2012, 109 (2), 020601. https://doi.org/10.1103/PhysRevLett.109.020601.
class colvar::gspath
: public colvar::CartesianBasedPath, public GeometricPathCV::GeometricPathBase<cvm::atom_pos, cvm::real, GeometricPathCV::path_sz::S>
{
private:
// Optimal rotation for compute v3
cvm::rotation rot_v3;
protected:
virtual void prepareVectors();
virtual void updateReferenceDistances();
public:
gspath(std::string const &conf);
virtual ~gspath() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
};
/// \brief Colvar component: alternative path collective variable using geometry, variable z
/// This should be merged with gspath in the same class by class inheritance or something else
class colvar::gzpath
: public colvar::CartesianBasedPath, public GeometricPathCV::GeometricPathBase<cvm::atom_pos, cvm::real, GeometricPathCV::path_sz::Z>
{
private:
// Optimal rotation for compute v3, v4
cvm::rotation rot_v3;
cvm::rotation rot_v4;
protected:
virtual void prepareVectors();
virtual void updateReferenceDistances();
public:
gzpath(std::string const &conf);
virtual ~gzpath() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
};
/// Current only linear combination of sub-CVCs is available
class colvar::linearCombination
: public colvar::cvc
{
protected:
/// Map from string to the types of colvar components
std::map<std::string, std::function<colvar::cvc* (const std::string& subcv_conf)>> string_cv_map;
/// Sub-colvar components
std::vector<colvar::cvc*> cv;
/// If all sub-cvs use explicit gradients then we also use it
bool use_explicit_gradients;
protected:
cvm::real getPolynomialFactorOfCVGradient(size_t i_cv) const;
public:
linearCombination(std::string const &conf);
virtual ~linearCombination();
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
};
class colvar::CVBasedPath
: public colvar::cvc
{
protected:
/// Map from string to the types of colvar components
std::map<std::string, std::function<colvar::cvc* (const std::string& subcv_conf)>> string_cv_map;
/// Sub-colvar components
std::vector<colvar::cvc*> cv;
/// Refernce colvar values from path
std::vector<std::vector<colvarvalue>> ref_cv;
/// If all sub-cvs use explicit gradients then we also use it
bool use_explicit_gradients;
/// Total number of reference frames
size_t total_reference_frames;
protected:
virtual void computeReferenceDistance(std::vector<cvm::real>& result);
cvm::real getPolynomialFactorOfCVGradient(size_t i_cv) const;
public:
CVBasedPath(std::string const &conf);
virtual ~CVBasedPath();
virtual void calc_value() = 0;
virtual void apply_force(colvarvalue const &force) = 0;
};
/// \brief Colvar component: alternative path collective variable using geometry, variable s
/// Allow any combination of existing (scalar) CVs
/// For more information see https://plumed.github.io/doc-v2.5/user-doc/html/_p_a_t_h.html
/// Díaz Leines, G.; Ensing, B. Path Finding on High-Dimensional Free Energy Landscapes. Phys. Rev. Lett. 2012, 109 (2), 020601. https://doi.org/10.1103/PhysRevLett.109.020601.
class colvar::gspathCV
: public colvar::CVBasedPath, public GeometricPathCV::GeometricPathBase<colvarvalue, cvm::real, GeometricPathCV::path_sz::S>
{
protected:
virtual void updateReferenceDistances();
virtual void prepareVectors();
public:
gspathCV(std::string const &conf);
virtual ~gspathCV();
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
};
class colvar::gzpathCV
: public colvar::CVBasedPath, public GeometricPathCV::GeometricPathBase<colvarvalue, cvm::real, GeometricPathCV::path_sz::Z>
{
protected:
virtual void updateReferenceDistances();
virtual void prepareVectors();
public:
gzpathCV(std::string const &conf);
virtual ~gzpathCV();
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
};
#else // if the compiler doesn't support C++11
class colvar::linearCombination
: public colvar::componentDisabled
{
public:
linearCombination(std::string const &conf) : componentDisabled(conf) {}
};
class colvar::CartesianBasedPath
: public colvar::componentDisabled
{
public:
CartesianBasedPath(std::string const &conf) : componentDisabled(conf) {}
};
class colvar::CVBasedPath
: public colvar::componentDisabled
{
public:
CVBasedPath(std::string const &conf) : componentDisabled(conf) {}
};
class colvar::gspath
: public colvar::componentDisabled
{
public:
gspath(std::string const &conf) : componentDisabled(conf) {}
};
class colvar::gzpath
: public colvar::componentDisabled
{
public:
gzpath(std::string const &conf) : componentDisabled(conf) {}
};
class colvar::gspathCV
: public colvar::componentDisabled
{
public:
gspathCV(std::string const &conf) : componentDisabled(conf) {}
};
class colvar::gzpathCV
: public colvar::componentDisabled
{
public:
gzpathCV(std::string const &conf) : componentDisabled(conf) {}
};
#endif // C++11 checking
// metrics functions for cvc implementations
// simple definitions of the distance functions; these are useful only

8
lib/colvars/colvarcomp_distances.cpp
View File

@ -1013,8 +1013,7 @@ colvar::rmsd::rmsd(std::string const &conf)
cvm::to_str(atoms->size())+").\n");
return;
}
}
{
} else { // Only look for ref pos file if ref positions not already provided
std::string ref_pos_file;
if (get_keyval(conf, "refPositionsFile", ref_pos_file, std::string(""))) {
@ -1041,12 +1040,15 @@ colvar::rmsd::rmsd(std::string const &conf)
cvm::load_coords(ref_pos_file.c_str(), &ref_pos, atoms,
ref_pos_col, ref_pos_col_value);
} else {
cvm::error("Error: no reference positions for RMSD; use either refPositions of refPositionsFile.");
return;
}
}
if (ref_pos.size() != atoms->size()) {
cvm::error("Error: found " + cvm::to_str(ref_pos.size()) +
" reference positions; expected " + cvm::to_str(atoms->size()));
" reference positions for RMSD; expected " + cvm::to_str(atoms->size()));
return;
}

910
lib/colvars/colvarcomp_gpath.cpp Normal file
View File

@ -0,0 +1,910 @@
#if (__cplusplus >= 201103L)
// 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
// Please update all Colvars source files before making any changes.
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include <numeric>
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <limits>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
namespace GeometricPathCV {
void init_string_cv_map(std::map<std::string, std::function<colvar::cvc* (const std::string& conf)>>& string_cv_map);
}
bool compareColvarComponent(colvar::cvc *i, colvar::cvc *j)
{
return i->name < j->name;
}
colvar::CartesianBasedPath::CartesianBasedPath(std::string const &conf): cvc(conf), atoms(nullptr), reference_frames(0) {
// Parse selected atoms
atoms = parse_group(conf, "atoms");
has_user_defined_fitting = false;
std::string fitting_conf;
if (key_lookup(conf, "fittingAtoms", &fitting_conf)) {
has_user_defined_fitting = true;
}
// Lookup reference column of PDB
// Copied from the RMSD class
std::string reference_column;
double reference_column_value;
if (get_keyval(conf, "refPositionsCol", reference_column, std::string(""))) {
bool found = get_keyval(conf, "refPositionsColValue", reference_column_value, 0.0);
if (found && reference_column_value == 0.0) {
cvm::error("Error: refPositionsColValue, "
"if provided, must be non-zero.\n");
return;
}
}
// Lookup all reference frames
bool has_frames = true;
total_reference_frames = 0;
while (has_frames) {
std::string reference_position_file_lookup = "refPositionsFile" + cvm::to_str(total_reference_frames + 1);
if (key_lookup(conf, reference_position_file_lookup.c_str())) {
std::string reference_position_filename;
get_keyval(conf, reference_position_file_lookup.c_str(), reference_position_filename, std::string(""));
std::vector<cvm::atom_pos> reference_position(atoms->size());
cvm::load_coords(reference_position_filename.c_str(), &reference_position, atoms, reference_column, reference_column_value);
reference_frames.push_back(reference_position);
++total_reference_frames;
} else {
has_frames = false;
}
}
// Setup alignment to compute RMSD with respect to reference frames
for (size_t i_frame = 0; i_frame < reference_frames.size(); ++i_frame) {
cvm::atom_group* tmp_atoms = parse_group(conf, "atoms");
if (!has_user_defined_fitting) {
// Swipe from the rmsd class
tmp_atoms->b_center = true;
tmp_atoms->b_rotate = true;
tmp_atoms->ref_pos = reference_frames[i_frame];
tmp_atoms->center_ref_pos();
tmp_atoms->enable(f_ag_fit_gradients);
tmp_atoms->rot.request_group1_gradients(tmp_atoms->size());
tmp_atoms->rot.request_group2_gradients(tmp_atoms->size());
comp_atoms.push_back(tmp_atoms);
} else {
// parse a group of atoms for fitting
std::string fitting_group_name = std::string("fittingAtoms") + cvm::to_str(i_frame);
cvm::atom_group* tmp_fitting_atoms = new cvm::atom_group(fitting_group_name.c_str());
tmp_fitting_atoms->parse(fitting_conf);
tmp_fitting_atoms->disable(f_ag_scalable);
tmp_fitting_atoms->disable(f_ag_scalable_com);
tmp_fitting_atoms->fit_gradients.assign(tmp_fitting_atoms->size(), cvm::atom_pos(0.0, 0.0, 0.0));
std::string reference_position_file_lookup = "refPositionsFile" + cvm::to_str(i_frame + 1);
std::string reference_position_filename;
get_keyval(conf, reference_position_file_lookup.c_str(), reference_position_filename, std::string(""));
std::vector<cvm::atom_pos> reference_fitting_position(tmp_fitting_atoms->size());
cvm::load_coords(reference_position_filename.c_str(), &reference_fitting_position, tmp_fitting_atoms, reference_column, reference_column_value);
// setup the atom group for calculating
tmp_atoms->b_center = true;
tmp_atoms->b_rotate = true;
tmp_atoms->b_user_defined_fit = true;
tmp_atoms->disable(f_ag_scalable);
tmp_atoms->disable(f_ag_scalable_com);
tmp_atoms->ref_pos = reference_fitting_position;
tmp_atoms->center_ref_pos();
tmp_atoms->enable(f_ag_fit_gradients);
tmp_atoms->enable(f_ag_fitting_group);
tmp_atoms->fitting_group = tmp_fitting_atoms;
tmp_atoms->rot.request_group1_gradients(tmp_fitting_atoms->size());
tmp_atoms->rot.request_group2_gradients(tmp_fitting_atoms->size());
reference_fitting_frames.push_back(reference_fitting_position);
comp_atoms.push_back(tmp_atoms);
}
}
x.type(colvarvalue::type_scalar);
// Don't use implicit gradient
enable(f_cvc_explicit_gradient);
}
colvar::CartesianBasedPath::~CartesianBasedPath() {
for (auto it_comp_atoms = comp_atoms.begin(); it_comp_atoms != comp_atoms.end(); ++it_comp_atoms) {
if (*it_comp_atoms != nullptr) {
delete (*it_comp_atoms);
(*it_comp_atoms) = nullptr;
}
}
}
void colvar::CartesianBasedPath::computeReferenceDistance(std::vector<cvm::real>& result) {
for (size_t i_frame = 0; i_frame < reference_frames.size(); ++i_frame) {
cvm::real frame_rmsd = 0.0;
for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
frame_rmsd += ((*(comp_atoms[i_frame]))[i_atom].pos - reference_frames[i_frame][i_atom]).norm2();
}
frame_rmsd /= cvm::real(atoms->size());
frame_rmsd = cvm::sqrt(frame_rmsd);
result[i_frame] = frame_rmsd;
}
}
colvar::gspath::gspath(std::string const &conf): CartesianBasedPath(conf) {
function_type = "gspath";
get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true);
if (use_second_closest_frame == true) {
cvm::log(std::string("Geometric path s(σ) will use the second closest frame to compute s_(m-1)\n"));
} else {
cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m-1)\n"));
}
get_keyval(conf, "useThirdClosestFrame", use_third_closest_frame, false);
if (use_third_closest_frame == true) {
cvm::log(std::string("Geometric path s(σ) will use the third closest frame to compute s_(m+1)\n"));
} else {
cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m+1)\n"));
}
GeometricPathCV::GeometricPathBase<cvm::atom_pos, cvm::real, GeometricPathCV::path_sz::S>::initialize(atoms->size(), cvm::atom_pos(), total_reference_frames, use_second_closest_frame, use_third_closest_frame);
cvm::log(std::string("Geometric pathCV(s) is initialized.\n"));
cvm::log(std::string("Geometric pathCV(s) loaded ") + cvm::to_str(reference_frames.size()) + std::string(" frames.\n"));
}
void colvar::gspath::updateReferenceDistances() {
computeReferenceDistance(frame_distances);
}
void colvar::gspath::prepareVectors() {
size_t i_atom;
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
// v1 = s_m - z
v1[i_atom] = reference_frames[min_frame_index_1][i_atom] - (*(comp_atoms[min_frame_index_1]))[i_atom].pos;
// v2 = z - s_(m-1)
v2[i_atom] = (*(comp_atoms[min_frame_index_2]))[i_atom].pos - reference_frames[min_frame_index_2][i_atom];
}
if (min_frame_index_3 < 0 || min_frame_index_3 > M) {
cvm::atom_pos reference_cog_1, reference_cog_2;
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
reference_cog_1 += reference_frames[min_frame_index_1][i_atom];
reference_cog_2 += reference_frames[min_frame_index_2][i_atom];
}
reference_cog_1 /= reference_frames[min_frame_index_1].size();
reference_cog_2 /= reference_frames[min_frame_index_2].size();
std::vector<cvm::atom_pos> tmp_reference_frame_1(reference_frames[min_frame_index_1].size());
std::vector<cvm::atom_pos> tmp_reference_frame_2(reference_frames[min_frame_index_2].size());
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
tmp_reference_frame_1[i_atom] = reference_frames[min_frame_index_1][i_atom] - reference_cog_1;
tmp_reference_frame_2[i_atom] = reference_frames[min_frame_index_2][i_atom] - reference_cog_2;
}
if (has_user_defined_fitting) {
cvm::atom_pos reference_fitting_cog_1, reference_fitting_cog_2;
for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
reference_fitting_cog_1 += reference_fitting_frames[min_frame_index_1][i_atom];
reference_fitting_cog_2 += reference_fitting_frames[min_frame_index_2][i_atom];
}
reference_fitting_cog_1 /= reference_fitting_frames[min_frame_index_1].size();
reference_fitting_cog_2 /= reference_fitting_frames[min_frame_index_2].size();
std::vector<cvm::atom_pos> tmp_reference_fitting_frame_1(reference_fitting_frames[min_frame_index_1].size());
std::vector<cvm::atom_pos> tmp_reference_fitting_frame_2(reference_fitting_frames[min_frame_index_2].size());
for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
tmp_reference_fitting_frame_1[i_atom] = reference_fitting_frames[min_frame_index_1][i_atom] - reference_fitting_cog_1;
tmp_reference_fitting_frame_2[i_atom] = reference_fitting_frames[min_frame_index_2][i_atom] - reference_fitting_cog_2;
}
rot_v3.calc_optimal_rotation(tmp_reference_fitting_frame_1, tmp_reference_fitting_frame_2);
} else {
rot_v3.calc_optimal_rotation(tmp_reference_frame_1, tmp_reference_frame_2);
}
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
v3[i_atom] = rot_v3.q.rotate(tmp_reference_frame_1[i_atom]) - tmp_reference_frame_2[i_atom];
}
} else {
cvm::atom_pos reference_cog_1, reference_cog_3;
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
reference_cog_1 += reference_frames[min_frame_index_1][i_atom];
reference_cog_3 += reference_frames[min_frame_index_3][i_atom];
}
reference_cog_1 /= reference_frames[min_frame_index_1].size();
reference_cog_3 /= reference_frames[min_frame_index_3].size();
std::vector<cvm::atom_pos> tmp_reference_frame_1(reference_frames[min_frame_index_1].size());
std::vector<cvm::atom_pos> tmp_reference_frame_3(reference_frames[min_frame_index_3].size());
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
tmp_reference_frame_1[i_atom] = reference_frames[min_frame_index_1][i_atom] - reference_cog_1;
tmp_reference_frame_3[i_atom] = reference_frames[min_frame_index_3][i_atom] - reference_cog_3;
}
if (has_user_defined_fitting) {
cvm::atom_pos reference_fitting_cog_1, reference_fitting_cog_3;
for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
reference_fitting_cog_1 += reference_fitting_frames[min_frame_index_1][i_atom];
reference_fitting_cog_3 += reference_fitting_frames[min_frame_index_3][i_atom];
}
reference_fitting_cog_1 /= reference_fitting_frames[min_frame_index_1].size();
reference_fitting_cog_3 /= reference_fitting_frames[min_frame_index_3].size();
std::vector<cvm::atom_pos> tmp_reference_fitting_frame_1(reference_fitting_frames[min_frame_index_1].size());
std::vector<cvm::atom_pos> tmp_reference_fitting_frame_3(reference_fitting_frames[min_frame_index_3].size());
for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
tmp_reference_fitting_frame_1[i_atom] = reference_fitting_frames[min_frame_index_1][i_atom] - reference_fitting_cog_1;
tmp_reference_fitting_frame_3[i_atom] = reference_fitting_frames[min_frame_index_3][i_atom] - reference_fitting_cog_3;
}
rot_v3.calc_optimal_rotation(tmp_reference_fitting_frame_1, tmp_reference_fitting_frame_3);
} else {
rot_v3.calc_optimal_rotation(tmp_reference_frame_1, tmp_reference_frame_3);
}
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
// v3 = s_(m+1) - s_m
v3[i_atom] = tmp_reference_frame_3[i_atom] - rot_v3.q.rotate(tmp_reference_frame_1[i_atom]);
}
}
}
void colvar::gspath::calc_value() {
computeValue();
x = s;
}
void colvar::gspath::calc_gradients() {
computeDerivatives();
cvm::rvector tmp_atom_grad_v1, tmp_atom_grad_v2;
// dS(v1, v2(r), v3) / dr = ∂S/∂v1 * dv1/dr + ∂S/∂v2 * dv2/dr
// dv1/dr = [fitting matrix 1][-1, ..., -1]
// dv2/dr = [fitting matrix 2][1, ..., 1]
// ∂S/∂v1 = ± (∂f/∂v1) / (2M)
// ∂S/∂v2 = ± (∂f/∂v2) / (2M)
// dS(v1, v2(r), v3) / dr = -1.0 * ± (∂f/∂v1) / (2M) + ± (∂f/∂v2) / (2M)
for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
tmp_atom_grad_v1[0] = -1.0 * sign * 0.5 * dfdv1[i_atom][0] / M;
tmp_atom_grad_v1[1] = -1.0 * sign * 0.5 * dfdv1[i_atom][1] / M;
tmp_atom_grad_v1[2] = -1.0 * sign * 0.5 * dfdv1[i_atom][2] / M;
tmp_atom_grad_v2[0] = sign * 0.5 * dfdv2[i_atom][0] / M;
tmp_atom_grad_v2[1] = sign * 0.5 * dfdv2[i_atom][1] / M;
tmp_atom_grad_v2[2] = sign * 0.5 * dfdv2[i_atom][2] / M;
(*(comp_atoms[min_frame_index_1]))[i_atom].grad += tmp_atom_grad_v1;
(*(comp_atoms[min_frame_index_2]))[i_atom].grad += tmp_atom_grad_v2;
}
}
void colvar::gspath::apply_force(colvarvalue const &force) {
// The force applied to this CV is scalar type
cvm::real const &F = force.real_value;
(*(comp_atoms[min_frame_index_1])).apply_colvar_force(F);
(*(comp_atoms[min_frame_index_2])).apply_colvar_force(F);
}
colvar::gzpath::gzpath(std::string const &conf): CartesianBasedPath(conf) {
function_type = "gzpath";
get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true);
if (use_second_closest_frame == true) {
cvm::log(std::string("Geometric path z(σ) will use the second closest frame to compute s_(m-1)\n"));
} else {
cvm::log(std::string("Geometric path z(σ) will use the neighbouring frame to compute s_(m-1)\n"));
}
get_keyval(conf, "useThirdClosestFrame", use_third_closest_frame, false);
if (use_third_closest_frame == true) {
cvm::log(std::string("Geometric path z(σ) will use the third closest frame to compute s_(m+1)\n"));
} else {
cvm::log(std::string("Geometric path z(σ) will use the neighbouring frame to compute s_(m+1)\n"));
}
bool b_use_z_square = false;
get_keyval(conf, "useZsquare", b_use_z_square, false);
if (b_use_z_square == true) {
cvm::log(std::string("Geometric path z(σ) will use the square of distance from current frame to path compute z\n"));
}
GeometricPathCV::GeometricPathBase<cvm::atom_pos, cvm::real, GeometricPathCV::path_sz::Z>::initialize(atoms->size(), cvm::atom_pos(), total_reference_frames, use_second_closest_frame, use_third_closest_frame, b_use_z_square);
// Logging
cvm::log(std::string("Geometric pathCV(z) is initialized.\n"));
cvm::log(std::string("Geometric pathCV(z) loaded ") + cvm::to_str(reference_frames.size()) + std::string(" frames.\n"));
}
void colvar::gzpath::updateReferenceDistances() {
computeReferenceDistance(frame_distances);
}
void colvar::gzpath::prepareVectors() {
cvm::atom_pos reference_cog_1, reference_cog_2;
size_t i_atom;
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
reference_cog_1 += reference_frames[min_frame_index_1][i_atom];
reference_cog_2 += reference_frames[min_frame_index_2][i_atom];
}
reference_cog_1 /= reference_frames[min_frame_index_1].size();
reference_cog_2 /= reference_frames[min_frame_index_2].size();
std::vector<cvm::atom_pos> tmp_reference_frame_1(reference_frames[min_frame_index_1].size());
std::vector<cvm::atom_pos> tmp_reference_frame_2(reference_frames[min_frame_index_2].size());
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
tmp_reference_frame_1[i_atom] = reference_frames[min_frame_index_1][i_atom] - reference_cog_1;
tmp_reference_frame_2[i_atom] = reference_frames[min_frame_index_2][i_atom] - reference_cog_2;
}
std::vector<cvm::atom_pos> tmp_reference_fitting_frame_1;
std::vector<cvm::atom_pos> tmp_reference_fitting_frame_2;
if (has_user_defined_fitting) {
cvm::atom_pos reference_fitting_cog_1, reference_fitting_cog_2;
for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
reference_fitting_cog_1 += reference_fitting_frames[min_frame_index_1][i_atom];
reference_fitting_cog_2 += reference_fitting_frames[min_frame_index_2][i_atom];
}
reference_fitting_cog_1 /= reference_fitting_frames[min_frame_index_1].size();
reference_fitting_cog_2 /= reference_fitting_frames[min_frame_index_2].size();
tmp_reference_fitting_frame_1.resize(reference_fitting_frames[min_frame_index_1].size());
tmp_reference_fitting_frame_2.resize(reference_fitting_frames[min_frame_index_2].size());
for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
tmp_reference_fitting_frame_1[i_atom] = reference_fitting_frames[min_frame_index_1][i_atom] - reference_fitting_cog_1;
tmp_reference_fitting_frame_2[i_atom] = reference_fitting_frames[min_frame_index_2][i_atom] - reference_fitting_cog_2;
}
rot_v4.calc_optimal_rotation(tmp_reference_fitting_frame_1, tmp_reference_fitting_frame_2);
} else {
rot_v4.calc_optimal_rotation(tmp_reference_frame_1, tmp_reference_frame_2);
}
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
v1[i_atom] = reference_frames[min_frame_index_1][i_atom] - (*(comp_atoms[min_frame_index_1]))[i_atom].pos;
v2[i_atom] = (*(comp_atoms[min_frame_index_2]))[i_atom].pos - reference_frames[min_frame_index_2][i_atom];
// v4 only computes in gzpath
// v4 = s_m - s_(m-1)
v4[i_atom] = rot_v4.q.rotate(tmp_reference_frame_1[i_atom]) - tmp_reference_frame_2[i_atom];
}
if (min_frame_index_3 < 0 || min_frame_index_3 > M) {
v3 = v4;
} else {
cvm::atom_pos reference_cog_3;
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
reference_cog_3 += reference_frames[min_frame_index_3][i_atom];
}
reference_cog_3 /= reference_frames[min_frame_index_3].size();
std::vector<cvm::atom_pos> tmp_reference_frame_3(reference_frames[min_frame_index_3].size());
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
tmp_reference_frame_3[i_atom] = reference_frames[min_frame_index_3][i_atom] - reference_cog_3;
}
if (has_user_defined_fitting) {
cvm::atom_pos reference_fitting_cog_3;
for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_3].size(); ++i_atom) {
reference_fitting_cog_3 += reference_fitting_frames[min_frame_index_3][i_atom];
}
reference_fitting_cog_3 /= reference_fitting_frames[min_frame_index_3].size();
std::vector<cvm::atom_pos> tmp_reference_fitting_frame_3(reference_fitting_frames[min_frame_index_3].size());
for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_3].size(); ++i_atom) {
tmp_reference_fitting_frame_3[i_atom] = reference_fitting_frames[min_frame_index_3][i_atom] - reference_fitting_cog_3;
}
rot_v3.calc_optimal_rotation(tmp_reference_fitting_frame_1, tmp_reference_fitting_frame_3);
} else {
rot_v3.calc_optimal_rotation(tmp_reference_frame_1, tmp_reference_frame_3);
}
for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
// v3 = s_(m+1) - s_m
v3[i_atom] = tmp_reference_frame_3[i_atom] - rot_v3.q.rotate(tmp_reference_frame_1[i_atom]);
}
}
}
void colvar::gzpath::calc_value() {
computeValue();
x = z;
}
void colvar::gzpath::calc_gradients() {
computeDerivatives();
cvm::rvector tmp_atom_grad_v1, tmp_atom_grad_v2;
for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
tmp_atom_grad_v1 = -1.0 * dzdv1[i_atom];
tmp_atom_grad_v2 = dzdv2[i_atom];
(*(comp_atoms[min_frame_index_1]))[i_atom].grad += tmp_atom_grad_v1;
(*(comp_atoms[min_frame_index_2]))[i_atom].grad += tmp_atom_grad_v2;
}
}
void colvar::gzpath::apply_force(colvarvalue const &force) {
// The force applied to this CV is scalar type
cvm::real const &F = force.real_value;
(*(comp_atoms[min_frame_index_1])).apply_colvar_force(F);
(*(comp_atoms[min_frame_index_2])).apply_colvar_force(F);
}
colvar::linearCombination::linearCombination(std::string const &conf): cvc(conf) {
GeometricPathCV::init_string_cv_map(string_cv_map);
// Lookup all available sub-cvcs
for (auto it_cv_map = string_cv_map.begin(); it_cv_map != string_cv_map.end(); ++it_cv_map) {
if (key_lookup(conf, it_cv_map->first.c_str())) {
std::vector<std::string> sub_cvc_confs;
get_key_string_multi_value(conf, it_cv_map->first.c_str(), sub_cvc_confs);
for (auto it_sub_cvc_conf = sub_cvc_confs.begin(); it_sub_cvc_conf != sub_cvc_confs.end(); ++it_sub_cvc_conf) {
cv.push_back((it_cv_map->second)(*(it_sub_cvc_conf)));
}
}
}
// Sort all sub CVs by their names
std::sort(cv.begin(), cv.end(), compareColvarComponent);
for (auto it_sub_cv = cv.begin(); it_sub_cv != cv.end(); ++it_sub_cv) {
for (auto it_atom_group = (*it_sub_cv)->atom_groups.begin(); it_atom_group != (*it_sub_cv)->atom_groups.end(); ++it_atom_group) {
register_atom_group(*it_atom_group);
}
}
x.type(cv[0]->value());
x.reset();
use_explicit_gradients = true;
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
if (!cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
use_explicit_gradients = false;
}
}
if (!use_explicit_gradients) {
disable(f_cvc_explicit_gradient);
}
}
cvm::real colvar::linearCombination::getPolynomialFactorOfCVGradient(size_t i_cv) const {
cvm::real factor_polynomial = 1.0;
if (cv[i_cv]->value().type() == colvarvalue::type_scalar) {
factor_polynomial = cv[i_cv]->sup_coeff * cv[i_cv]->sup_np * cvm::pow(cv[i_cv]->value().real_value, cv[i_cv]->sup_np - 1);
} else {
factor_polynomial = cv[i_cv]->sup_coeff;
}
return factor_polynomial;
}
colvar::linearCombination::~linearCombination() {
for (auto it = cv.begin(); it != cv.end(); ++it) {
delete (*it);
}
}
void colvar::linearCombination::calc_value() {
x.reset();
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
cv[i_cv]->calc_value();
colvarvalue current_cv_value(cv[i_cv]->value());
// polynomial combination allowed
if (current_cv_value.type() == colvarvalue::type_scalar) {
x += cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np));
} else {
x += cv[i_cv]->sup_coeff * current_cv_value;
}
}
}
void colvar::linearCombination::calc_gradients() {
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
cv[i_cv]->calc_gradients();
if ( cv[i_cv]->is_enabled(f_cvc_explicit_gradient) &&
!cv[i_cv]->is_enabled(f_cvc_scalable) &&
!cv[i_cv]->is_enabled(f_cvc_scalable_com)) {
cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) {
for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) {
(*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = factor_polynomial * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad;
}
}
}
}
}
}
void colvar::linearCombination::apply_force(colvarvalue const &force) {
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
// If this CV us explicit gradients, then atomic gradients is already calculated
// We can apply the force to atom groups directly
if ( cv[i_cv]->is_enabled(f_cvc_explicit_gradient) &&
!cv[i_cv]->is_enabled(f_cvc_scalable) &&
!cv[i_cv]->is_enabled(f_cvc_scalable_com)
) {
for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
(cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value);
}
} else {
// Compute factors for polynomial combinations
cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
colvarvalue cv_force = force.real_value * factor_polynomial;
cv[i_cv]->apply_force(cv_force);
}
}
}
colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) {
GeometricPathCV::init_string_cv_map(string_cv_map);
// Lookup all available sub-cvcs
for (auto it_cv_map = string_cv_map.begin(); it_cv_map != string_cv_map.end(); ++it_cv_map) {
if (key_lookup(conf, it_cv_map->first.c_str())) {
std::vector<std::string> sub_cvc_confs;
get_key_string_multi_value(conf, it_cv_map->first.c_str(), sub_cvc_confs);
for (auto it_sub_cvc_conf = sub_cvc_confs.begin(); it_sub_cvc_conf != sub_cvc_confs.end(); ++it_sub_cvc_conf) {
cv.push_back((it_cv_map->second)(*(it_sub_cvc_conf)));
}
}
}
// Sort all sub CVs by their names
std::sort(cv.begin(), cv.end(), compareColvarComponent);
// Register atom groups and determine the colvar type for reference
std::vector<colvarvalue> tmp_cv;
for (auto it_sub_cv = cv.begin(); it_sub_cv != cv.end(); ++it_sub_cv) {
for (auto it_atom_group = (*it_sub_cv)->atom_groups.begin(); it_atom_group != (*it_sub_cv)->atom_groups.end(); ++it_atom_group) {
register_atom_group(*it_atom_group);
}
colvarvalue tmp_i_cv((*it_sub_cv)->value());
tmp_i_cv.reset();
tmp_cv.push_back(tmp_i_cv);
}
// Read path file
// Lookup all reference CV values
std::string path_filename;
get_keyval(conf, "pathFile", path_filename);
cvm::log(std::string("Reading path file: ") + path_filename + std::string("\n"));
std::ifstream ifs_path(path_filename);
if (!ifs_path.is_open()) {
cvm::error("Error: failed to open path file.\n");
}
std::string line;
const std::string token(" ");
total_reference_frames = 0;
while (std::getline(ifs_path, line)) {
std::vector<std::string> fields;
split_string(line, token, fields);
size_t num_value_required = 0;
for (size_t i_cv = 0; i_cv < tmp_cv.size(); ++i_cv) {
const size_t value_size = tmp_cv[i_cv].size();
num_value_required += value_size;
cvm::log(std::string("Reading CV ") + cv[i_cv]->name + std::string(" with ") + cvm::to_str(value_size) + std::string(" value(s)\n"));
if (num_value_required <= fields.size()) {
size_t start_index = num_value_required - value_size;
for (size_t i = start_index; i < num_value_required; ++i) {
tmp_cv[i_cv][i] = std::atof(fields[i].c_str());
cvm::log(fields[i] + std::string(" "));
}
cvm::log(std::string("\n"));
} else {
cvm::error("Error: incorrect format of path file.\n");
}
}
if (!fields.empty()) {
ref_cv.push_back(tmp_cv);
++total_reference_frames;
}
}
x.type(colvarvalue::type_scalar);
use_explicit_gradients = true;
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
if (!cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
use_explicit_gradients = false;
}
}
if (!use_explicit_gradients) {
disable(f_cvc_explicit_gradient);
}
}
void colvar::CVBasedPath::computeReferenceDistance(std::vector<cvm::real>& result) {
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
cv[i_cv]->calc_value();
}
for (size_t i_frame = 0; i_frame < ref_cv.size(); ++i_frame) {
cvm::real rmsd_i = 0.0;
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
colvarvalue ref_cv_value(ref_cv[i_frame][i_cv]);
colvarvalue current_cv_value(cv[i_cv]->value());
// polynomial combination allowed
if (current_cv_value.type() == colvarvalue::type_scalar) {
// wrapping is already in dist2
rmsd_i += cv[i_cv]->dist2(cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)), ref_cv_value.real_value);
} else {
rmsd_i += cv[i_cv]->dist2(cv[i_cv]->sup_coeff * current_cv_value, ref_cv_value);
}
}
rmsd_i /= cvm::real(cv.size());
rmsd_i = cvm::sqrt(rmsd_i);
result[i_frame] = rmsd_i;
}
}
cvm::real colvar::CVBasedPath::getPolynomialFactorOfCVGradient(size_t i_cv) const {
cvm::real factor_polynomial = 1.0;
if (cv[i_cv]->value().type() == colvarvalue::type_scalar) {
factor_polynomial = cv[i_cv]->sup_coeff * cv[i_cv]->sup_np * cvm::pow(cv[i_cv]->value().real_value, cv[i_cv]->sup_np - 1);
} else {
factor_polynomial = cv[i_cv]->sup_coeff;
}
return factor_polynomial;
}
colvar::CVBasedPath::~CVBasedPath() {
for (auto it = cv.begin(); it != cv.end(); ++it) {
delete (*it);
}
}
colvar::gspathCV::gspathCV(std::string const &conf): CVBasedPath(conf) {
function_type = "gspathCV";
cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n"));
// Initialize variables for future calculation
get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true);
if (use_second_closest_frame == true) {
cvm::log(std::string("Geometric path s(σ) will use the second closest frame to compute s_(m-1)\n"));
} else {
cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m-1)\n"));
}
get_keyval(conf, "useThirdClosestFrame", use_third_closest_frame, false);
if (use_third_closest_frame == true) {
cvm::log(std::string("Geometric path s(σ) will use the third closest frame to compute s_(m+1)\n"));
} else {
cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m+1)\n"));
}
GeometricPathCV::GeometricPathBase<colvarvalue, cvm::real, GeometricPathCV::path_sz::S>::initialize(cv.size(), ref_cv[0], total_reference_frames, use_second_closest_frame, use_third_closest_frame);
x.type(colvarvalue::type_scalar);
use_explicit_gradients = true;
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
if (!cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
use_explicit_gradients = false;
}
}
if (!use_explicit_gradients) {
cvm::log("Geometric path s(σ) will use implicit gradients.\n");
disable(f_cvc_explicit_gradient);
}
}
colvar::gspathCV::~gspathCV() {}
void colvar::gspathCV::updateReferenceDistances() {
computeReferenceDistance(frame_distances);
}
void colvar::gspathCV::prepareVectors() {
// Compute v1, v2 and v3
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
// values of sub-cvc are computed in update_distances
// cv[i_cv]->calc_value();
colvarvalue f1_ref_cv_i_value(ref_cv[min_frame_index_1][i_cv]);
colvarvalue f2_ref_cv_i_value(ref_cv[min_frame_index_2][i_cv]);
colvarvalue current_cv_value(cv[i_cv]->value());
// polynomial combination allowed
if (current_cv_value.type() == colvarvalue::type_scalar) {
v1[i_cv] = f1_ref_cv_i_value.real_value - cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np));
v2[i_cv] = cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)) - f2_ref_cv_i_value.real_value;
} else {
v1[i_cv] = f1_ref_cv_i_value - cv[i_cv]->sup_coeff * current_cv_value;
v2[i_cv] = cv[i_cv]->sup_coeff * current_cv_value - f2_ref_cv_i_value;
}
cv[i_cv]->wrap(v1[i_cv]);
cv[i_cv]->wrap(v2[i_cv]);
}
if (min_frame_index_3 < 0 || min_frame_index_3 > M) {
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
v3[i_cv] = ref_cv[min_frame_index_1][i_cv] - ref_cv[min_frame_index_2][i_cv];
cv[i_cv]->wrap(v3[i_cv]);
}
} else {
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
v3[i_cv] = ref_cv[min_frame_index_3][i_cv] - ref_cv[min_frame_index_1][i_cv];
cv[i_cv]->wrap(v3[i_cv]);
}
}
}
void colvar::gspathCV::calc_value() {
computeValue();
x = s;
}
void colvar::gspathCV::calc_gradients() {
computeDerivatives();
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
// No matter whether the i-th cv uses implicit gradient, compute it first.
cv[i_cv]->calc_gradients();
// If the gradient is not implicit, then add the gradients to its atom groups
if ( cv[i_cv]->is_enabled(f_cvc_explicit_gradient) &&
!cv[i_cv]->is_enabled(f_cvc_scalable) &&
!cv[i_cv]->is_enabled(f_cvc_scalable_com)) {
// Temporary variables storing gradients
colvarvalue tmp_cv_grad_v1(cv[i_cv]->value());
colvarvalue tmp_cv_grad_v2(cv[i_cv]->value());
// Compute factors for polynomial combinations
cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
// Loop over all elements of the corresponding colvar value
for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) {
// ds/dz, z = vector of CVs
tmp_cv_grad_v1[j_elem] = -1.0 * sign * 0.5 * dfdv1[i_cv][j_elem] / M;
tmp_cv_grad_v2[j_elem] = sign * 0.5 * dfdv2[i_cv][j_elem] / M;
// Apply the gradients to the atom groups in i-th cv
// Loop over all atom groups
for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
// Loop over all atoms in the k-th atom group
for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) {
// Chain rule
(*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = factor_polynomial * ((*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad * tmp_cv_grad_v1[j_elem] + (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad * tmp_cv_grad_v2[j_elem]);
}
}
}
}
}
}
void colvar::gspathCV::apply_force(colvarvalue const &force) {
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
// If this CV us explicit gradients, then atomic gradients is already calculated
// We can apply the force to atom groups directly
if ( cv[i_cv]->is_enabled(f_cvc_explicit_gradient) &&
!cv[i_cv]->is_enabled(f_cvc_scalable) &&
!cv[i_cv]->is_enabled(f_cvc_scalable_com)
) {
for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
(cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value);
}
} else {
// Temporary variables storing gradients
colvarvalue tmp_cv_grad_v1(cv[i_cv]->value());
colvarvalue tmp_cv_grad_v2(cv[i_cv]->value());
// Compute factors for polynomial combinations
cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) {
// ds/dz, z = vector of CVs
tmp_cv_grad_v1[j_elem] = -1.0 * sign * 0.5 * dfdv1[i_cv][j_elem] / M;
tmp_cv_grad_v2[j_elem] = sign * 0.5 * dfdv2[i_cv][j_elem] / M;
}
colvarvalue cv_force = force.real_value * factor_polynomial * (tmp_cv_grad_v1 + tmp_cv_grad_v2);
cv[i_cv]->apply_force(cv_force);
}
}
}
colvar::gzpathCV::gzpathCV(std::string const &conf): CVBasedPath(conf) {
function_type = "gzpathCV";
cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n"));
// Initialize variables for future calculation
M = cvm::real(total_reference_frames - 1);
m = 1.0;
get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true);
if (use_second_closest_frame == true) {
cvm::log(std::string("Geometric path z(σ) will use the second closest frame to compute s_(m-1)\n"));
} else {
cvm::log(std::string("Geometric path z(σ) will use the neighbouring frame to compute s_(m-1)\n"));
}
get_keyval(conf, "useThirdClosestFrame", use_third_closest_frame, false);
if (use_third_closest_frame == true) {
cvm::log(std::string("Geometric path z(σ) will use the third closest frame to compute s_(m+1)\n"));
} else {
cvm::log(std::string("Geometric path z(σ) will use the neighbouring frame to compute s_(m+1)\n"));
}
bool b_use_z_square = false;
get_keyval(conf, "useZsquare", b_use_z_square, false);
if (b_use_z_square == true) {
cvm::log(std::string("Geometric path z(σ) will use the square of distance from current frame to path compute z\n"));
}
GeometricPathCV::GeometricPathBase<colvarvalue, cvm::real, GeometricPathCV::path_sz::Z>::initialize(cv.size(), ref_cv[0], total_reference_frames, use_second_closest_frame, use_third_closest_frame, b_use_z_square);
x.type(colvarvalue::type_scalar);
use_explicit_gradients = true;
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
if (!cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
use_explicit_gradients = false;
}
}
if (!use_explicit_gradients) {
cvm::log("Geometric path z(σ) will use implicit gradients.\n");
disable(f_cvc_explicit_gradient);
}
}
colvar::gzpathCV::~gzpathCV() {
}
void colvar::gzpathCV::updateReferenceDistances() {
computeReferenceDistance(frame_distances);
}
void colvar::gzpathCV::prepareVectors() {
// Compute v1, v2 and v3
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
// values of sub-cvc are computed in update_distances
// cv[i_cv]->calc_value();
colvarvalue f1_ref_cv_i_value(ref_cv[min_frame_index_1][i_cv]);
colvarvalue f2_ref_cv_i_value(ref_cv[min_frame_index_2][i_cv]);
colvarvalue current_cv_value(cv[i_cv]->value());
// polynomial combination allowed
if (current_cv_value.type() == colvarvalue::type_scalar) {
v1[i_cv] = f1_ref_cv_i_value.real_value - cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np));
v2[i_cv] = cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)) - f2_ref_cv_i_value.real_value;
} else {
v1[i_cv] = f1_ref_cv_i_value - cv[i_cv]->sup_coeff * current_cv_value;
v2[i_cv] = cv[i_cv]->sup_coeff * current_cv_value - f2_ref_cv_i_value;
}
v4[i_cv] = f1_ref_cv_i_value - f2_ref_cv_i_value;
cv[i_cv]->wrap(v1[i_cv]);
cv[i_cv]->wrap(v2[i_cv]);
cv[i_cv]->wrap(v4[i_cv]);
}
if (min_frame_index_3 < 0 || min_frame_index_3 > M) {
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
v3[i_cv] = ref_cv[min_frame_index_1][i_cv] - ref_cv[min_frame_index_2][i_cv];
cv[i_cv]->wrap(v3[i_cv]);
}
} else {
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
v3[i_cv] = ref_cv[min_frame_index_3][i_cv] - ref_cv[min_frame_index_1][i_cv];
cv[i_cv]->wrap(v3[i_cv]);
}
}
}
void colvar::gzpathCV::calc_value() {
computeValue();
x = z;
}
void colvar::gzpathCV::calc_gradients() {
computeDerivatives();
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
// No matter whether the i-th cv uses implicit gradient, compute it first.
cv[i_cv]->calc_gradients();
// If the gradient is not implicit, then add the gradients to its atom groups
if ( cv[i_cv]->is_enabled(f_cvc_explicit_gradient) &&
!cv[i_cv]->is_enabled(f_cvc_scalable) &&
!cv[i_cv]->is_enabled(f_cvc_scalable_com)) {
// Temporary variables storing gradients
colvarvalue tmp_cv_grad_v1 = -1.0 * dzdv1[i_cv];
colvarvalue tmp_cv_grad_v2 = 1.0 * dzdv2[i_cv];
// Compute factors for polynomial combinations
cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) {
// Apply the gradients to the atom groups in i-th cv
// Loop over all atom groups
for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
// Loop over all atoms in the k-th atom group
for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) {
// Chain rule
(*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = factor_polynomial * ((*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad * tmp_cv_grad_v1[j_elem] + (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad * tmp_cv_grad_v2[j_elem]);
}
}
}
}
}
}
void colvar::gzpathCV::apply_force(colvarvalue const &force) {
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
// If this CV us explicit gradients, then atomic gradients is already calculated
// We can apply the force to atom groups directly
if ( cv[i_cv]->is_enabled(f_cvc_explicit_gradient) &&
!cv[i_cv]->is_enabled(f_cvc_scalable) &&
!cv[i_cv]->is_enabled(f_cvc_scalable_com)) {
for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
(cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value);
}
}
else {
colvarvalue tmp_cv_grad_v1 = -1.0 * dzdv1[i_cv];
colvarvalue tmp_cv_grad_v2 = 1.0 * dzdv2[i_cv];
// Temporary variables storing gradients
// Compute factors for polynomial combinations
cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
colvarvalue cv_force = force.real_value * factor_polynomial * (tmp_cv_grad_v1 + tmp_cv_grad_v2);
cv[i_cv]->apply_force(cv_force);
}
}
}
void GeometricPathCV::init_string_cv_map(std::map<std::string, std::function<colvar::cvc* (const std::string& subcv_conf)>>& string_cv_map) {
string_cv_map["distance"] = [](const std::string& conf){return new colvar::distance(conf);};
string_cv_map["dihedral"] = [](const std::string& conf){return new colvar::dihedral(conf);};
string_cv_map["angle"] = [](const std::string& conf){return new colvar::angle(conf);};
string_cv_map["rmsd"] = [](const std::string& conf){return new colvar::rmsd(conf);};
string_cv_map["gyration"] = [](const std::string& conf){return new colvar::gyration(conf);};
string_cv_map["inertia"] = [](const std::string& conf){return new colvar::inertia(conf);};
string_cv_map["inertiaZ"] = [](const std::string& conf){return new colvar::inertia_z(conf);};
string_cv_map["tilt"] = [](const std::string& conf){return new colvar::tilt(conf);};
string_cv_map["distanceZ"] = [](const std::string& conf){return new colvar::distance_z(conf);};
string_cv_map["distanceXY"] = [](const std::string& conf){return new colvar::distance_xy(conf);};
string_cv_map["polarTheta"] = [](const std::string& conf){return new colvar::polar_theta(conf);};
string_cv_map["polarPhi"] = [](const std::string& conf){return new colvar::polar_phi(conf);};
string_cv_map["distanceVec"] = [](const std::string& conf){return new colvar::distance_vec(conf);};
string_cv_map["orientationAngle"] = [](const std::string& conf){return new colvar::orientation_angle(conf);};
string_cv_map["distancePairs"] = [](const std::string& conf){return new colvar::distance_pairs(conf);};
string_cv_map["dipoleMagnitude"] = [](const std::string& conf){return new colvar::dipole_magnitude(conf);};
string_cv_map["coordNum"] = [](const std::string& conf){return new colvar::coordnum(conf);};
string_cv_map["selfCoordNum"] = [](const std::string& conf){return new colvar::selfcoordnum(conf);};
string_cv_map["dipoleAngle"] = [](const std::string& conf){return new colvar::dipole_angle(conf);};
string_cv_map["orientation"] = [](const std::string& conf){return new colvar::orientation(conf);};
string_cv_map["orientationProj"] = [](const std::string& conf){return new colvar::orientation_proj(conf);};
string_cv_map["eigenvector"] = [](const std::string& conf){return new colvar::eigenvector(conf);};
string_cv_map["cartesian"] = [](const std::string& conf){return new colvar::cartesian(conf);};
string_cv_map["alpha"] = [](const std::string& conf){return new colvar::alpha_angles(conf);};
string_cv_map["dihedralPC"] = [](const std::string& conf){return new colvar::dihedPC(conf);};
string_cv_map["linearCombination"] = [](const std::string& conf){return new colvar::linearCombination(conf);};
}
#endif

2
lib/colvars/colvarcomp_protein.cpp
View File

@ -7,6 +7,8 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include <algorithm>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"

4
lib/colvars/colvargrid.cpp
View File

@ -109,7 +109,9 @@ cvm::real colvar_grid_scalar::entropy() const
{
cvm::real sum = 0.0;
for (size_t i = 0; i < nt; i++) {
sum += -1.0 * data[i] * cvm::logn(data[i]);
if (data[i] >0) {
sum += -1.0 * data[i] * cvm::logn(data[i]);
}
}
cvm::real bin_volume = 1.0;
for (size_t id = 0; id < widths.size(); id++) {

6
lib/colvars/colvargrid.h
View File

@ -557,11 +557,11 @@ public:
data[i] *= a;
}
/// \brief Assign all zero elements a scalar constant (fast loop)
inline void remove_zeros(cvm::real const &a)
/// \brief Assign values that are smaller than scalar constant the latter value (fast loop)
inline void remove_small_values(cvm::real const &a)
{
for (size_t i = 0; i < nt; i++)
if(data[i]==0) data[i] = a;
if(data[i]<a) data[i] = a;
}

38
lib/colvars/colvarparse.cpp
View File

@ -7,7 +7,6 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include <sstream>
#include <iostream>
#include <algorithm>
@ -58,6 +57,28 @@ bool colvarparse::get_key_string_value(std::string const &conf,
return b_found_any;
}
bool colvarparse::get_key_string_multi_value(std::string const &conf,
char const *key, std::vector<std::string>& data)
{
bool b_found = false, b_found_any = false;
size_t save_pos = 0, found_count = 0;
data.clear();
do {
std::string data_this = "";
b_found = key_lookup(conf, key, &data_this, &save_pos);
if (b_found) {
if (!b_found_any)
b_found_any = true;
found_count++;
data.push_back(data_this);
}
} while (b_found);
return b_found_any;
}
template<typename TYPE>
void colvarparse::mark_key_set_user(std::string const &key_str,
@ -843,3 +864,18 @@ int colvarparse::check_braces(std::string const &conf,
}
return (brace_count != 0) ? INPUT_ERROR : COLVARS_OK;
}
void colvarparse::split_string(const std::string& data, const std::string& delim, std::vector<std::string>& dest) {
size_t index = 0, new_index = 0;
std::string tmpstr;
while (index != data.length()) {
new_index = data.find(delim, index);
if (new_index != std::string::npos) tmpstr = data.substr(index, new_index - index);
else tmpstr = data.substr(index, data.length());
if (!tmpstr.empty()) {
dest.push_back(tmpstr);
}
if (new_index == std::string::npos) break;
index = new_index + 1;
}
}

10
lib/colvars/colvarparse.h
View File

@ -212,6 +212,10 @@ protected:
bool get_key_string_value(std::string const &conf,
char const *key, std::string &data);
/// Get multiple strings from repeated instances of a same keyword
bool get_key_string_multi_value(std::string const &conf,
char const *key, std::vector<std::string>& data);
/// Template for single-value keyword parsers
template<typename TYPE>
bool _get_keyval_scalar_(std::string const &conf,
@ -322,6 +326,12 @@ public:
/// from this position
static int check_braces(std::string const &conf, size_t const start_pos);
/// \brief Split a string with a specified delimiter into a vector
/// \param data The string to be splitted
/// \param delim A delimiter
/// \param dest A destination vector to store the splitted results
static void split_string(const std::string& data, const std::string& delim, std::vector<std::string>& dest);
protected:
/// \brief List of legal keywords for this object: this is updated

77
lib/colvars/colvarproxy.cpp
View File

@ -7,8 +7,14 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#if !defined(WIN32) || defined(__CYGWIN__)
#include <unistd.h>
#endif
#include <cerrno>
#include <sstream>
#include <cstring>
#include <cstdio>
#if defined(_OPENMP)
#include <omp.h>
@ -704,6 +710,27 @@ std::ostream * colvarproxy_io::output_stream(std::string const &output_name,
if (cvm::debug()) {
cvm::log("Using colvarproxy::output_stream()\n");
}
std::ostream *os = get_output_stream(output_name);
if (os != NULL) return os;
if (!(mode & (std::ios_base::app | std::ios_base::ate))) {
backup_file(output_name);
}
std::ofstream *osf = new std::ofstream(output_name.c_str(), mode);
if (!osf->is_open()) {
cvm::error("Error: cannot write to file/channel \""+output_name+"\".\n",
FILE_ERROR);
return NULL;
}
output_stream_names.push_back(output_name);
output_files.push_back(osf);
return osf;
}
std::ostream *colvarproxy_io::get_output_stream(std::string const &output_name)
{
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
@ -711,21 +738,11 @@ std::ostream * colvarproxy_io::output_stream(std::string const &output_name,
return *osi;
}
}
if (!(mode & (std::ios_base::app | std::ios_base::ate))) {
backup_file(output_name);
}
std::ofstream *os = new std::ofstream(output_name.c_str(), mode);
if (!os->is_open()) {
cvm::error("Error: cannot write to file/channel \""+output_name+"\".\n",
FILE_ERROR);
return NULL;
}
output_stream_names.push_back(output_name);
output_files.push_back(os);
return os;
return NULL;
}
int colvarproxy_io::flush_output_stream(std::ostream *os)
{
std::list<std::ostream *>::iterator osi = output_files.begin();
@ -761,10 +778,45 @@ int colvarproxy_io::close_output_stream(std::string const &output_name)
int colvarproxy_io::backup_file(char const *filename)
{
// TODO implement this using rename_file()
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_io::remove_file(char const *filename)
{
if (std::remove(filename)) {
if (errno != ENOENT) {
return cvm::error("Error: in removing file \""+std::string(filename)+
"\".\n.",
FILE_ERROR);
}
}
return COLVARS_OK;
}
int colvarproxy_io::rename_file(char const *filename, char const *newfilename)
{
int error_code = COLVARS_OK;
#if defined(WIN32) && !defined(__CYGWIN__)
// On straight Windows, must remove the destination before renaming it
error_code |= remove_file(newfilename);
#endif
int rename_exit_code = 0;
while ((rename_exit_code = std::rename(filename, newfilename)) != 0) {
if (errno == EINTR) continue;
// Call log() instead of error to allow the next try
cvm::log("Error: in renaming file \""+std::string(filename)+"\" to \""+
std::string(newfilename)+"\".\n.");
error_code |= FILE_ERROR;
if (errno == EXDEV) continue;
break;
}
return rename_exit_code ? error_code : COLVARS_OK;
}
colvarproxy::colvarproxy()
{
@ -826,4 +878,3 @@ int colvarproxy::get_version_from_string(char const *version_string)
is >> newint;
return newint;
}

24
lib/colvars/colvarproxy.h
View File

@ -571,7 +571,10 @@ public:
/// if this is not open already, then open it
virtual std::ostream *output_stream(std::string const &output_name,
std::ios_base::openmode mode =
std::ios_base::out);
std::ios_base::out);
/// Returns a reference to output_name if it exists, NULL otherwise
virtual std::ostream *get_output_stream(std::string const &output_name);
/// \brief Flushes the given output channel
virtual int flush_output_stream(std::ostream *os);
@ -588,6 +591,25 @@ public:
return backup_file(filename.c_str());
}
/// Remove the given file
int remove_file(char const *filename);
/// Remove the given file
inline int remove_file(std::string const &filename)
{
return remove_file(filename.c_str());
}
/// Rename the given file
int rename_file(char const *filename, char const *newfilename);
/// Rename the given file
inline int rename_file(std::string const &filename,
std::string const &newfilename)
{
return rename_file(filename.c_str(), newfilename.c_str());
}
/// \brief Prefix of the input state file
inline std::string & input_prefix()
{

2
lib/colvars/colvars_version.h
View File

@ -1,5 +1,5 @@
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2019-04-26"
#define COLVARS_VERSION "2019-08-01"
// 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

2
src/USER-COLVARS/colvarproxy_lammps_version.h
View File

@ -1,5 +1,5 @@
#ifndef COLVARPROXY_VERSION
#define COLVARPROXY_VERSION "2019-04-09"
#define COLVARPROXY_VERSION "2019-08-01"
// 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