diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf index 6194424212..14deceeb87 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/Makefile.common b/lib/colvars/Makefile.common index e3fa4662e6..5958765077 100644 --- a/lib/colvars/Makefile.common +++ b/lib/colvars/Makefile.common @@ -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 \ diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index 723a54ad39..e3676084ac 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -791,6 +791,11 @@ int colvar::init_components(std::string const &conf) "inertia", "inertia"); error_code |= init_components_type(conf, "moment of inertia around an axis", "inertiaZ"); error_code |= init_components_type(conf, "eigenvector", "eigenvector"); + error_code |= init_components_type(conf, "geometrical path collective variables (s)", "gspath"); + error_code |= init_components_type(conf, "geometrical path collective variables (z)", "gzpath"); + error_code |= init_components_type(conf, "linear combination of other collective variables", "subColvar"); + error_code |= init_components_type(conf, "geometrical path collective variables (s) for other CVs", "gspathCV"); + error_code |= init_components_type(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)) { diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index 74f7fdee51..e2ab0f3c1d 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -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; diff --git a/lib/colvars/colvar_geometricpath.h b/lib/colvars/colvar_geometricpath.h new file mode 100644 index 0000000000..a7ef7f337c --- /dev/null +++ b/lib/colvars/colvar_geometricpath.h @@ -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 +#include +#include +#include +#include +#include + +namespace GeometricPathCV { + +enum path_sz {S, Z}; + +template +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 v1; + std::vector v2; + std::vector v3; + std::vector v4; + std::vector dfdv1; + std::vector dfdv2; + std::vector dzdv1; + std::vector dzdv2; + std::vector frame_distances; + std::vector 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& 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& 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 +GeometricPathBase::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 +GeometricPathBase::GeometricPathBase(size_t vector_size, const std::vector& 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 +void GeometricPathBase::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(total_frames - 1); + m = static_cast(1.0); +} + +template +void GeometricPathBase::initialize(size_t vector_size, const std::vector& 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(total_frames - 1); + m = static_cast(1.0); +} + +template +void GeometricPathBase::prepareVectors() { + std::cout << "Warning: you should not call the prepareVectors() in base class!\n"; + std::cout << std::flush; +} + +template +void GeometricPathBase::updateReferenceDistances() { + std::cout << "Warning: you should not call the updateReferenceDistances() in base class!\n"; + std::cout << std::flush; +} + +template +void GeometricPathBase::compute() { + computeValue(); + computeDerivatives(); +} + +template +void GeometricPathBase::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(frame_index[0]) - static_cast(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(frame_index[0]) - static_cast(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(frame_index[0]); +} + +template +void GeometricPathBase::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(sign) * ((f - 1) / (2 * M)); + } +} + +template +void GeometricPathBase::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(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 diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index 27781ec733..8540e4a945 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -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::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::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; } - - diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index 0ba2bef1c3..8e98274b2d 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -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; diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp index 90588f5a1f..ab02820cf0 100644 --- a/lib/colvars/colvarbias_restraint.cpp +++ b/lib/colvars/colvarbias_restraint.cpp @@ -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(); } diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index 3075ed82ca..431884a877 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -7,6 +7,8 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. +#include + #include "colvarmodule.h" #include "colvarvalue.h" #include "colvar.h" diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index f615680ba6..5e7a91e872 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -24,6 +24,13 @@ #include "colvar.h" #include "colvaratoms.h" +#if (__cplusplus >= 201103L) +#include "colvar_geometricpath.h" +#include +#endif // C++11 checking + +#include + /// \brief Colvar component (base class for collective variables) /// @@ -688,9 +695,6 @@ protected: /// Reference coordinates std::vector 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 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& result); + /// Selected atoms + cvm::atom_group *atoms; + /// Fitting options + bool has_user_defined_fitting; + /// Reference frames + std::vector> reference_frames; + std::vector> reference_fitting_frames; + /// Atom groups for RMSD calculation together with reference frames + std::vector 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 +{ +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 +{ +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> string_cv_map; + /// Sub-colvar components + std::vector 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> string_cv_map; + /// Sub-colvar components + std::vector cv; + /// Refernce colvar values from path + std::vector> 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& 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 +{ +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 +{ +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 diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index d9cd9d55e4..7cb4c30c31 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -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; } diff --git a/lib/colvars/colvarcomp_gpath.cpp b/lib/colvars/colvarcomp_gpath.cpp new file mode 100644 index 0000000000..530aaf3284 --- /dev/null +++ b/lib/colvars/colvarcomp_gpath.cpp @@ -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 +#include +#include +#include +#include + +#include "colvarmodule.h" +#include "colvarvalue.h" +#include "colvarparse.h" +#include "colvar.h" +#include "colvarcomp.h" + +namespace GeometricPathCV { +void init_string_cv_map(std::map>& 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 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 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& 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::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 tmp_reference_frame_1(reference_frames[min_frame_index_1].size()); + std::vector 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 tmp_reference_fitting_frame_1(reference_fitting_frames[min_frame_index_1].size()); + std::vector 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 tmp_reference_frame_1(reference_frames[min_frame_index_1].size()); + std::vector 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 tmp_reference_fitting_frame_1(reference_fitting_frames[min_frame_index_1].size()); + std::vector 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::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 tmp_reference_frame_1(reference_frames[min_frame_index_1].size()); + std::vector 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 tmp_reference_fitting_frame_1; + std::vector 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 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 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 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 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 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 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& 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::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::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>& 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 diff --git a/lib/colvars/colvarcomp_protein.cpp b/lib/colvars/colvarcomp_protein.cpp index b9f9c60cdb..aa61cdf1dc 100644 --- a/lib/colvars/colvarcomp_protein.cpp +++ b/lib/colvars/colvarcomp_protein.cpp @@ -7,6 +7,8 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. +#include + #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp index dc1a709edb..bc5ac3c57c 100644 --- a/lib/colvars/colvargrid.cpp +++ b/lib/colvars/colvargrid.cpp @@ -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++) { diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h index 2ba0566e49..7792b7d0df 100644 --- a/lib/colvars/colvargrid.h +++ b/lib/colvars/colvargrid.h @@ -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] #include #include @@ -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& 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 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& 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; + } +} diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h index 8501ee8c14..f4ead26601 100644 --- a/lib/colvars/colvarparse.h +++ b/lib/colvars/colvarparse.h @@ -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& data); + /// Template for single-value keyword parsers template 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& dest); + protected: /// \brief List of legal keywords for this object: this is updated diff --git a/lib/colvars/colvarproxy.cpp b/lib/colvars/colvarproxy.cpp index 5f8e82d30a..f077984556 100644 --- a/lib/colvars/colvarproxy.cpp +++ b/lib/colvars/colvarproxy.cpp @@ -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 +#endif +#include + #include #include +#include #if defined(_OPENMP) #include @@ -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::iterator osi = output_files.begin(); std::list::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::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; } - diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index 3bbdfe522e..21944106f1 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -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() { diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h index 2521fdf872..0f1ee9a0a0 100644 --- a/lib/colvars/colvars_version.h +++ b/lib/colvars/colvars_version.h @@ -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 diff --git a/src/USER-COLVARS/colvarproxy_lammps_version.h b/src/USER-COLVARS/colvarproxy_lammps_version.h index 0a4f9fdf4f..d976f40510 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 "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