From d49bed3c27ea4ed2a375b413aa9d83bd8e1b2edd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=B4me=20H=C3=A9nin?= Date: Mon, 16 Dec 2024 16:36:51 +0100 Subject: [PATCH] Fixes for the Colvars library - 759 min_image fix; Saves long runs from crashes; https://github.com/Colvars/colvars/pull/759 (@PolyachenkoYA) - 728 Fix undefined behavior when getting the current working directory from std::filesystem https://github.com/Colvars/colvars/pull/728 (@giacomofiorin) - 724 Fix gradients and metric functions of distanceDir https://github.com/Colvars/colvars/pull/724 (@giacomofiorin) - 715 Add missing rotation in orientation component https://github.com/Colvars/colvars/pull/715 (@giacomofiorin) - 713 fix: try to solve #87 for non-scala components https://github.com/Colvars/colvars/pull/713 (@HanatoK) - 706 BUGFIX for Segmentation fault in colvarbias_meta::calc_energy() with useGrids off https://github.com/Colvars/colvars/pull/706 (@alphataubio) - 701 Do not try accessing LAMMPS proxy object before allocating it https://github.com/Colvars/colvars/pull/701 (@giacomofiorin) Authors: @alphataubio, @giacomofiorin, @HanatoK, @PolyachenkoYA --- lib/colvars/colvaratoms.cpp | 125 +++++++++++++++++++----- lib/colvars/colvaratoms.h | 55 ++++++++++- lib/colvars/colvarbias_meta.cpp | 85 ++++++++--------- lib/colvars/colvarcomp_distances.cpp | 29 +++--- lib/colvars/colvarcomp_rotations.cpp | 9 +- lib/colvars/colvargrid.cpp | 2 +- lib/colvars/colvargrid.h | 2 +- lib/colvars/colvarmodule.h | 2 +- lib/colvars/colvarproxy_io.cpp | 56 ++++++++++- lib/colvars/colvarproxy_io.h | 6 ++ lib/colvars/colvarvalue.cpp | 137 ++++++++++++++++++++------- lib/colvars/colvarvalue.h | 86 ----------------- 12 files changed, 383 insertions(+), 211 deletions(-) diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index e15b9301a12..745cb032816 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -1217,23 +1217,30 @@ void cvm::atom_group::calc_fit_gradients() if (cvm::debug()) cvm::log("Calculating fit gradients.\n"); + cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this; + + auto accessor_main = [this](size_t i){return atoms[i].grad;}; + auto accessor_fitting = [&group_for_fit](size_t j, const cvm::rvector& grad){group_for_fit->fit_gradients[j] = grad;}; if (is_enabled(f_ag_center) && is_enabled(f_ag_rotate)) - calc_fit_gradients_impl(); + calc_fit_forces_impl(accessor_main, accessor_fitting); if (is_enabled(f_ag_center) && !is_enabled(f_ag_rotate)) - calc_fit_gradients_impl(); + calc_fit_forces_impl(accessor_main, accessor_fitting); if (!is_enabled(f_ag_center) && is_enabled(f_ag_rotate)) - calc_fit_gradients_impl(); + calc_fit_forces_impl(accessor_main, accessor_fitting); if (!is_enabled(f_ag_center) && !is_enabled(f_ag_rotate)) - calc_fit_gradients_impl(); + calc_fit_forces_impl(accessor_main, accessor_fitting); if (cvm::debug()) cvm::log("Done calculating fit gradients.\n"); } -template -void cvm::atom_group::calc_fit_gradients_impl() { - cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this; +template +void cvm::atom_group::calc_fit_forces_impl( + main_force_accessor_T accessor_main, + fitting_force_accessor_T accessor_fitting) const { + const cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this; // the center of geometry contribution to the gradients cvm::rvector atom_grad; // the rotation matrix contribution to the gradients @@ -1245,15 +1252,15 @@ void cvm::atom_group::calc_fit_gradients_impl() { for (size_t i = 0; i < size(); i++) { cvm::atom_pos pos_orig; if (B_ag_center) { - atom_grad += atoms[i].grad; + atom_grad += accessor_main(i); if (B_ag_rotate) pos_orig = rot_inv * (atoms[i].pos - ref_pos_cog); } else { - if (B_ag_rotate) pos_orig = atoms[i].pos; + if (B_ag_rotate) pos_orig = rot_inv * atoms[i].pos; } if (B_ag_rotate) { // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i cvm::quaternion const dxdq = - rot.q.position_derivative_inner(pos_orig, atoms[i].grad); + rot.q.position_derivative_inner(pos_orig, accessor_main(i)); sum_dxdq[0] += dxdq[0]; sum_dxdq[1] += dxdq[1]; sum_dxdq[2] += dxdq[2]; @@ -1261,26 +1268,45 @@ void cvm::atom_group::calc_fit_gradients_impl() { } } if (B_ag_center) { - if (B_ag_rotate) atom_grad = rot.inverse().matrix() * atom_grad; + if (B_ag_rotate) atom_grad = rot_inv * atom_grad; atom_grad *= (-1.0)/(cvm::real(group_for_fit->size())); } // loop 2: iterate over the fitting group if (B_ag_rotate) rot_deriv->prepare_derivative(rotation_derivative_dldq::use_dq); for (size_t j = 0; j < group_for_fit->size(); j++) { + cvm::rvector fitting_force_grad{0, 0, 0}; if (B_ag_center) { - group_for_fit->fit_gradients[j] = atom_grad; + fitting_force_grad += atom_grad; } if (B_ag_rotate) { rot_deriv->calc_derivative_wrt_group1(j, nullptr, &dq0_1); // multiply by {\partial q}/\partial\vec{x}_j and add it to the fit gradients - group_for_fit->fit_gradients[j] += sum_dxdq[0] * dq0_1[0] + - sum_dxdq[1] * dq0_1[1] + - sum_dxdq[2] * dq0_1[2] + - sum_dxdq[3] * dq0_1[3]; + fitting_force_grad += sum_dxdq[0] * dq0_1[0] + + sum_dxdq[1] * dq0_1[1] + + sum_dxdq[2] * dq0_1[2] + + sum_dxdq[3] * dq0_1[3]; + } + if (cvm::debug()) { + cvm::log(cvm::to_str(fitting_force_grad)); } + accessor_fitting(j, fitting_force_grad); } } +template +void cvm::atom_group::calc_fit_forces( + main_force_accessor_T accessor_main, + fitting_force_accessor_T accessor_fitting) const { + if (is_enabled(f_ag_center) && is_enabled(f_ag_rotate)) + calc_fit_forces_impl(accessor_main, accessor_fitting); + if (is_enabled(f_ag_center) && !is_enabled(f_ag_rotate)) + calc_fit_forces_impl(accessor_main, accessor_fitting); + if (!is_enabled(f_ag_center) && is_enabled(f_ag_rotate)) + calc_fit_forces_impl(accessor_main, accessor_fitting); + if (!is_enabled(f_ag_center) && !is_enabled(f_ag_rotate)) + calc_fit_forces_impl(accessor_main, accessor_fitting); +} + std::vector cvm::atom_group::positions() const { @@ -1452,17 +1478,72 @@ void cvm::atom_group::apply_force(cvm::rvector const &force) return; } - if (is_enabled(f_ag_rotate)) { + auto ag_force = get_group_force_object(); + for (size_t i = 0; i < size(); ++i) { + ag_force.add_atom_force(i, atoms[i].mass / total_mass * force); + } +} - const auto rot_inv = rot.inverse().matrix(); - for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { - ai->apply_force(rot_inv * ((ai->mass/total_mass) * force)); +cvm::atom_group::group_force_object cvm::atom_group::get_group_force_object() { + return cvm::atom_group::group_force_object(this); +} + +cvm::atom_group::group_force_object::group_force_object(cvm::atom_group* ag): +m_ag(ag), m_group_for_fit(m_ag->fitting_group ? m_ag->fitting_group : m_ag), +m_has_fitting_force(m_ag->is_enabled(f_ag_center) || m_ag->is_enabled(f_ag_rotate)) { + if (m_has_fitting_force) { + if (m_ag->group_forces.size() != m_ag->size()) { + m_ag->group_forces.assign(m_ag->size(), 0); + } else { + std::fill(m_ag->group_forces.begin(), + m_ag->group_forces.end(), 0); } + } +} +cvm::atom_group::group_force_object::~group_force_object() { + if (m_has_fitting_force) { + apply_force_with_fitting_group(); + } +} + +void cvm::atom_group::group_force_object::add_atom_force(size_t i, const cvm::rvector& force) { + if (m_has_fitting_force) { + m_ag->group_forces[i] += force; } else { + // Apply the force directly if we don't use fitting + (*m_ag)[i].apply_force(force); + } +} - for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { - ai->apply_force((ai->mass/total_mass) * force); +void cvm::atom_group::group_force_object::apply_force_with_fitting_group() { + const cvm::rmatrix rot_inv = m_ag->rot.inverse().matrix(); + if (cvm::debug()) { + cvm::log("Applying force on main group " + m_ag->name + ":\n"); + } + for (size_t ia = 0; ia < m_ag->size(); ++ia) { + const cvm::rvector f_ia = rot_inv * m_ag->group_forces[ia]; + (*m_ag)[ia].apply_force(f_ia); + if (cvm::debug()) { + cvm::log(cvm::to_str(f_ia)); + } + } + // Gradients are only available with scalar components, so for a scalar component, + // if f_ag_fit_gradients is disabled, then the forces on the fitting group is not + // computed. For a vector component, we can only know the forces on the fitting + // group, but checking this flag can mimic results that the users expect (if + // "enableFitGradients no" then there is no force on the fitting group). + if (m_ag->is_enabled(f_ag_fit_gradients)) { + auto accessor_main = [this](size_t i){return m_ag->group_forces[i];}; + auto accessor_fitting = [this](size_t j, const cvm::rvector& fitting_force){ + (*(m_group_for_fit))[j].apply_force(fitting_force); + }; + if (cvm::debug()) { + cvm::log("Applying force on the fitting group of main group" + m_ag->name + ":\n"); + } + m_ag->calc_fit_forces(accessor_main, accessor_fitting); + if (cvm::debug()) { + cvm::log("Done applying force on the fitting group of main group" + m_ag->name + ":\n"); } } } diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index d16ca7bd567..1493046e033 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -257,8 +257,27 @@ class colvarmodule::atom_group /// \brief Index in the colvarproxy arrays (if the group is scalable) int index; + /// \brief The temporary forces acting on the main group atoms. + /// Currently this is only used for calculating the fitting group forces for + /// non-scalar components. + std::vector group_forces; + public: + class group_force_object { + public: + group_force_object(cvm::atom_group* ag); + ~group_force_object(); + void add_atom_force(size_t i, const cvm::rvector& force); + private: + cvm::atom_group* m_ag; + cvm::atom_group* m_group_for_fit; + bool m_has_fitting_force; + void apply_force_with_fitting_group(); + }; + + group_force_object get_group_force_object(); + inline cvm::atom & operator [] (size_t const i) { return atoms[i]; @@ -497,15 +516,47 @@ class colvarmodule::atom_group /// \brief Calculate the derivatives of the fitting transformation void calc_fit_gradients(); -/*! @brief Actual implementation of `calc_fit_gradients`. The template is +/*! @brief Actual implementation of `calc_fit_gradients` and + * `calc_fit_forces`. The template is * used to avoid branching inside the loops in case that the CPU * branch prediction is broken (or further migration to GPU code). * @tparam B_ag_center Centered the reference to origin? This should follow * the value of `is_enabled(f_ag_center)`. * @tparam B_ag_rotate Calculate the optimal rotation? This should follow * the value of `is_enabled(f_ag_rotate)`. + * @tparam main_force_accessor_T The type of accessor of the main + * group forces or gradients. + * @tparam fitting_force_accessor_T The type of accessor of the fitting group + * forces or gradients. + * @param accessor_main The accessor of the main group forces or gradients. + * accessor_main(i) should return the i-th force or gradient of the + * main group. + * @param accessor_fitting The accessor of the fitting group forces or gradients. + * accessor_fitting(j, v) should store/apply the j-th atom gradient or + * force in the fitting group. + */ + template + void calc_fit_forces_impl( + main_force_accessor_T accessor_main, + fitting_force_accessor_T accessor_fitting) const; + +/*! @brief Calculate or apply the fitting group forces from the main group forces. + * @tparam main_force_accessor_T The type of accessor of the main + * group forces or gradients. + * @tparam fitting_force_accessor_T The type of accessor of the fitting group + * forces or gradients. + * @param accessor_main The accessor of the main group forces or gradients. + * accessor_main(i) should return the i-th force or gradient of the + * main group. + * @param accessor_fitting The accessor of the fitting group forces or gradients. + * accessor_fitting(j, v) should store/apply the j-th atom gradient or + * force in the fitting group. */ - template void calc_fit_gradients_impl(); + template + void calc_fit_forces( + main_force_accessor_T accessor_main, + fitting_force_accessor_T accessor_fitting) const; /// \brief Derivatives of the fitting transformation std::vector fit_gradients; diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index 905cd17883b..004dcbb2e63 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -11,24 +11,6 @@ #include #include -// Define function to get the absolute path of a replica file -#if defined(_WIN32) && !defined(__CYGWIN__) -#include -#define GETCWD(BUF, SIZE) ::_getcwd(BUF, SIZE) -#define PATHSEP "\\" -#else -#include -#define GETCWD(BUF, SIZE) ::getcwd(BUF, SIZE) -#define PATHSEP "/" -#endif - -#ifdef __cpp_lib_filesystem -// When std::filesystem is available, use it -#include -#undef GETCWD -#define GETCWD(BUF, SIZE) (std::filesystem::current_path().string().c_str()) -#endif - #include "colvarmodule.h" #include "colvarproxy.h" #include "colvar.h" @@ -451,8 +433,11 @@ int colvarbias_meta::update() error_code |= update_grid_params(); // add new biasing energy/forces error_code |= update_bias(); - // update grid content to reflect new bias - error_code |= update_grid_data(); + + if (use_grids) { + // update grid content to reflect new bias + error_code |= update_grid_data(); + } if (comm != single_replica && (cvm::step_absolute() % replica_update_freq) == 0) { @@ -670,11 +655,20 @@ int colvarbias_meta::calc_energy(std::vector const *values) replicas[ir]->bias_energy = 0.0; } - std::vector const curr_bin = values ? - hills_energy->get_colvars_index(*values) : - hills_energy->get_colvars_index(); + bool index_ok = false; + std::vector curr_bin; + + if (use_grids) { + + curr_bin = values ? + hills_energy->get_colvars_index(*values) : + hills_energy->get_colvars_index(); + + index_ok = hills_energy->index_ok(curr_bin); + + } - if (hills_energy->index_ok(curr_bin)) { + if ( index_ok ) { // index is within the grid: get the energy from there for (ir = 0; ir < replicas.size(); ir++) { @@ -723,11 +717,20 @@ int colvarbias_meta::calc_forces(std::vector const *values) } } - std::vector const curr_bin = values ? - hills_energy->get_colvars_index(*values) : - hills_energy->get_colvars_index(); + bool index_ok = false; + std::vector curr_bin; - if (hills_energy->index_ok(curr_bin)) { + if (use_grids) { + + curr_bin = values ? + hills_energy->get_colvars_index(*values) : + hills_energy->get_colvars_index(); + + index_ok = hills_energy->index_ok(curr_bin); + + } + + if ( index_ok ) { for (ir = 0; ir < replicas.size(); ir++) { cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value(curr_bin)); for (ic = 0; ic < num_variables(); ic++) { @@ -1718,29 +1721,17 @@ int colvarbias_meta::setup_output() if (comm == multiple_replicas) { - // TODO: one may want to specify the path manually for intricated filesystems? - char *pwd = new char[3001]; - if (GETCWD(pwd, 3000) == nullptr) { - if (pwd != nullptr) { // - delete[] pwd; - } - return cvm::error("Error: cannot get the path of the current working directory.\n", - COLVARS_BUG_ERROR); - } - + auto const pwd = cvm::main()->proxy->get_current_work_dir(); replica_list_file = - (std::string(pwd)+std::string(PATHSEP)+ - this->name+"."+replica_id+".files.txt"); + cvm::main()->proxy->join_paths(pwd, this->name + "." + replica_id + ".files.txt"); // replica_hills_file and replica_state_file are those written // by the current replica; within the mirror biases, they are // those by another replica - replica_hills_file = - (std::string(pwd)+std::string(PATHSEP)+ - cvm::output_prefix()+".colvars."+this->name+"."+replica_id+".hills"); - replica_state_file = - (std::string(pwd)+std::string(PATHSEP)+ - cvm::output_prefix()+".colvars."+this->name+"."+replica_id+".state"); - delete[] pwd; + replica_hills_file = cvm::main()->proxy->join_paths( + pwd, cvm::output_prefix() + ".colvars." + this->name + "." + replica_id + ".hills"); + + replica_state_file = cvm::main()->proxy->join_paths( + pwd, cvm::output_prefix() + ".colvars." + this->name + "." + replica_id + ".state"); // now register this replica diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index 319190c385a..1c7267f1e5e 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -384,32 +384,30 @@ void colvar::distance_dir::apply_force(colvarvalue const &force) cvm::real const iprod = force.rvector_value * x.rvector_value; cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value; - if (!group1->noforce) - group1->apply_force(-1.0 * force_tang); - - if (!group2->noforce) - group2->apply_force( force_tang); + if (!group1->noforce) { + group1->apply_force(-1.0 / dist_v.norm() * force_tang); + } + if (!group2->noforce) { + group2->apply_force( 1.0 / dist_v.norm() * force_tang); + } } -cvm::real colvar::distance_dir::dist2(colvarvalue const &x1, - colvarvalue const &x2) const +cvm::real colvar::distance_dir::dist2(colvarvalue const &x1, colvarvalue const &x2) const { - return (x1.rvector_value - x2.rvector_value).norm2(); + return x1.dist2(x2); } -colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const +colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const { - return colvarvalue((x1.rvector_value - x2.rvector_value), colvarvalue::type_unit3vectorderiv); + return x1.dist2_grad(x2); } -colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const +colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const { - return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vectorderiv); + return x2.dist2_grad(x1); } @@ -1403,11 +1401,12 @@ void colvar::cartesian::apply_force(colvarvalue const &force) size_t ia, j; if (!atoms->noforce) { cvm::rvector f; + auto ag_force = atoms->get_group_force_object(); for (ia = 0; ia < atoms->size(); ia++) { for (j = 0; j < dim; j++) { f[axes[j]] = force.vector1d_value[dim*ia + j]; } - (*atoms)[ia].apply_force(f); + ag_force.add_atom_force(ia, f); } } } diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp index a04ace851a0..12dd90e33be 100644 --- a/lib/colvars/colvarcomp_rotations.cpp +++ b/lib/colvars/colvarcomp_rotations.cpp @@ -137,11 +137,14 @@ void colvar::orientation::apply_force(colvarvalue const &force) if (!atoms->noforce) { rot_deriv_impl->prepare_derivative(rotation_derivative_dldq::use_dq); cvm::vector1d dq0_2; + auto ag_force = atoms->get_group_force_object(); for (size_t ia = 0; ia < atoms->size(); ia++) { rot_deriv_impl->calc_derivative_wrt_group2(ia, nullptr, &dq0_2); - for (size_t i = 0; i < 4; i++) { - (*atoms)[ia].apply_force(FQ[i] * dq0_2[i]); - } + const auto f_ia = FQ[0] * dq0_2[0] + + FQ[1] * dq0_2[1] + + FQ[2] * dq0_2[2] + + FQ[3] * dq0_2[3]; + ag_force.add_atom_force(ia, f_ia); } } } diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp index 11693a75874..36f83b31ae3 100644 --- a/lib/colvars/colvargrid.cpp +++ b/lib/colvars/colvargrid.cpp @@ -617,7 +617,7 @@ integrate_potential::integrate_potential(std::vector &colvars, std::sh } -integrate_potential::integrate_potential(std::shared_ptr gradients) +integrate_potential::integrate_potential(colvar_grid_gradient * gradients) : b_smoothed(false), gradients(gradients) { diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h index 4cbbb109614..0f168af07fb 100644 --- a/lib/colvars/colvargrid.h +++ b/lib/colvars/colvargrid.h @@ -1832,7 +1832,7 @@ class integrate_potential : public colvar_grid_scalar integrate_potential(std::vector &colvars, std::shared_ptr gradients); /// Constructor from a gradient grid (for processing grid files without a Colvars config) - integrate_potential(std::shared_ptr gradients); + integrate_potential(colvar_grid_gradient * gradients); /// \brief Calculate potential from divergence (in 2D); return number of steps int integrate(const int itmax, const cvm::real & tol, cvm::real & err, bool verbose = true); diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index fa84b1ad75d..cf551e3d7c3 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -84,7 +84,7 @@ class colvarmodule { int version_int = 0; /// Patch version number (non-zero in patch releases of other packages) - int patch_version_int = 0; + int patch_version_int = 4; public: diff --git a/lib/colvars/colvarproxy_io.cpp b/lib/colvars/colvarproxy_io.cpp index 4cfdfeec26f..08666d138e6 100644 --- a/lib/colvars/colvarproxy_io.cpp +++ b/lib/colvars/colvarproxy_io.cpp @@ -8,13 +8,20 @@ // Colvars repository at GitHub. // Using access() to check if a file exists (until we can assume C++14/17) -#if !defined(_WIN32) || defined(__CYGWIN__) +#if defined(_WIN32) && !defined(__CYGWIN__) +#include +#else #include #endif + #if defined(_WIN32) #include #endif +#ifdef __cpp_lib_filesystem +#include +#endif + #include #include @@ -64,6 +71,53 @@ int colvarproxy_io::set_frame(long int) } +std::string colvarproxy_io::get_current_work_dir() const +{ +#ifdef __cpp_lib_filesystem + + return std::filesystem::current_path().string(); + +#else + + // Legacy code + size_t constexpr buf_size = 3001; + char buf[buf_size]; + +#if defined(_WIN32) && !defined(__CYGWIN__) + char *getcwd_result = ::_getcwd(buf, buf_size); +#else + char *getcwd_result = ::getcwd(buf, buf_size); +#endif + + if (getcwd_result == nullptr) { + cvm::error("Error: cannot read the current working directory.\n", COLVARS_INPUT_ERROR); + return std::string(""); + } + + return std::string(getcwd_result); +#endif +} + + +std::string colvarproxy_io::join_paths(std::string const &path1, std::string const &path2) const +{ +#ifdef __cpp_lib_filesystem + + return (std::filesystem::path(path1) / std::filesystem::path(path2)).string(); + +#else + + // Legacy code +#if defined(_WIN32) && !defined(__CYGWIN__) + return (path1 + "\\" + path2); +#else + return (path1 + "/" + path2); +#endif + +#endif +} + + int colvarproxy_io::backup_file(char const *filename) { // Simplified version of NAMD_file_exists() diff --git a/lib/colvars/colvarproxy_io.h b/lib/colvars/colvarproxy_io.h index 726f915c970..eaf750366d6 100644 --- a/lib/colvars/colvarproxy_io.h +++ b/lib/colvars/colvarproxy_io.h @@ -38,6 +38,12 @@ class colvarproxy_io { // Returns error code virtual int set_frame(long int); + /// Get the current working directory of this process + std::string get_current_work_dir() const; + + /// Join two paths using the operating system's path separation + std::string join_paths(std::string const &path1, std::string const &path2) const; + /// \brief Rename the given file, before overwriting it virtual int backup_file(char const *filename); diff --git a/lib/colvars/colvarvalue.cpp b/lib/colvars/colvarvalue.cpp index 3b8077d2e77..0a7223d9e77 100644 --- a/lib/colvars/colvarvalue.cpp +++ b/lib/colvars/colvarvalue.cpp @@ -153,29 +153,6 @@ std::string const colvarvalue::type_keyword(Type t) } -size_t colvarvalue::num_df(Type t) -{ - switch (t) { - case colvarvalue::type_notset: - default: - return 0; break; - case colvarvalue::type_scalar: - return 1; break; - case colvarvalue::type_3vector: - return 3; break; - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return 2; break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return 3; break; - case colvarvalue::type_vector: - // the size of a vector is unknown without its object - return 0; break; - } -} - - size_t colvarvalue::num_dimensions(Type t) { switch (t) { @@ -591,6 +568,97 @@ cvm::real operator * (colvarvalue const &x1, } +cvm::real colvarvalue::norm2() const +{ + switch (value_type) { + case colvarvalue::type_scalar: + return (this->real_value)*(this->real_value); + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + return (this->rvector_value).norm2(); + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + return (this->quaternion_value).norm2(); + case colvarvalue::type_vector: + if (elem_types.size() > 0) { + // if we have information about non-scalar types, use it + cvm::real result = 0.0; + size_t i; + for (i = 0; i < elem_types.size(); i++) { + result += (this->get_elem(i)).norm2(); + } + return result; + } else { + return vector1d_value.norm2(); + } + break; + case colvarvalue::type_notset: + default: + return 0.0; + } +} + + +cvm::real colvarvalue::sum() const +{ + switch (value_type) { + case colvarvalue::type_scalar: + return (this->real_value); + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + return (this->rvector_value).x + (this->rvector_value).y + + (this->rvector_value).z; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + return (this->quaternion_value).q0 + (this->quaternion_value).q1 + + (this->quaternion_value).q2 + (this->quaternion_value).q3; + case colvarvalue::type_vector: + return (this->vector1d_value).sum(); + case colvarvalue::type_notset: + default: + return 0.0; + } +} + + +cvm::real colvarvalue::dist2(colvarvalue const &x2) const +{ + colvarvalue::check_types(*this, x2); + + switch (this->type()) { + case colvarvalue::type_scalar: + return (this->real_value - x2.real_value) * (this->real_value - x2.real_value); + case colvarvalue::type_3vector: + return (this->rvector_value - x2.rvector_value).norm2(); + case colvarvalue::type_unit3vector: { + cvm::rvector const &v1 = this->rvector_value; + cvm::rvector const &v2 = x2.rvector_value; + cvm::real const theta = cvm::acos(v1 * v2); + return theta * theta; + } + case colvarvalue::type_quaternion: + // angle between (*this) and x2 is the distance, the quaternion + // object has it implemented internally + return this->quaternion_value.dist2(x2.quaternion_value); + case colvarvalue::type_vector: + return (this->vector1d_value - x2.vector1d_value).norm2(); + case colvarvalue::type_unit3vectorderiv: + case colvarvalue::type_quaternionderiv: + cvm::error("Error: computing a squared-distance between two variables of type \"" + + type_desc(this->type()) + "\", for which it is not defined.\n", + COLVARS_BUG_ERROR); + case colvarvalue::type_notset: + default: + this->undef_op(); + return 0.0; + }; + + return 0.0; +} + + colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const { colvarvalue::check_types(*this, x2); @@ -600,25 +668,30 @@ colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const return 2.0 * (this->real_value - x2.real_value); case colvarvalue::type_3vector: return 2.0 * (this->rvector_value - x2.rvector_value); - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - { - cvm::rvector const &v1 = this->rvector_value; - cvm::rvector const &v2 = x2.rvector_value; - cvm::real const cos_t = v1 * v2; - return colvarvalue(2.0 * (cos_t * v1 - v2), colvarvalue::type_unit3vectorderiv); - } + case colvarvalue::type_unit3vector: { + cvm::rvector const &v1 = this->rvector_value; + cvm::rvector const &v2 = x2.rvector_value; + cvm::real const cos_t = v1 * v2; + return colvarvalue(2.0 * std::acos(cos_t) * -1.0 / cvm::sqrt(1.0 - cos_t * cos_t) * v2, + colvarvalue::type_unit3vectorderiv); + } case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: return this->quaternion_value.dist2_grad(x2.quaternion_value); case colvarvalue::type_vector: return colvarvalue(2.0 * (this->vector1d_value - x2.vector1d_value), colvarvalue::type_vector); break; + case colvarvalue::type_unit3vectorderiv: + case colvarvalue::type_quaternionderiv: + cvm::error("Error: computing a squared-distance gradient between two variables of type \"" + + type_desc(this->type()) + "\", for which it is not defined.\n", + COLVARS_BUG_ERROR); case colvarvalue::type_notset: default: this->undef_op(); return colvarvalue(colvarvalue::type_notset); }; + + return colvarvalue(colvarvalue::type_notset); } diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h index e8a6a849d35..61f1bf718b0 100644 --- a/lib/colvars/colvarvalue.h +++ b/lib/colvars/colvarvalue.h @@ -109,9 +109,6 @@ class colvarvalue { /// User keywords for specifying value types in the configuration static std::string const type_keyword(Type t); - /// Number of degrees of freedom for each supported type - static size_t num_df(Type t); - /// Number of dimensions for each supported type (used to allocate vector1d_value) static size_t num_dimensions(Type t); @@ -671,87 +668,4 @@ inline cvm::vector1d const colvarvalue::as_vector() const } -inline cvm::real colvarvalue::norm2() const -{ - switch (value_type) { - case colvarvalue::type_scalar: - return (this->real_value)*(this->real_value); - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return (this->rvector_value).norm2(); - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return (this->quaternion_value).norm2(); - case colvarvalue::type_vector: - if (elem_types.size() > 0) { - // if we have information about non-scalar types, use it - cvm::real result = 0.0; - size_t i; - for (i = 0; i < elem_types.size(); i++) { - result += (this->get_elem(i)).norm2(); - } - return result; - } else { - return vector1d_value.norm2(); - } - break; - case colvarvalue::type_notset: - default: - return 0.0; - } -} - - -inline cvm::real colvarvalue::sum() const -{ - switch (value_type) { - case colvarvalue::type_scalar: - return (this->real_value); - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return (this->rvector_value).x + (this->rvector_value).y + - (this->rvector_value).z; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return (this->quaternion_value).q0 + (this->quaternion_value).q1 + - (this->quaternion_value).q2 + (this->quaternion_value).q3; - case colvarvalue::type_vector: - return (this->vector1d_value).sum(); - case colvarvalue::type_notset: - default: - return 0.0; - } -} - - -inline cvm::real colvarvalue::dist2(colvarvalue const &x2) const -{ - colvarvalue::check_types(*this, x2); - - switch (this->type()) { - case colvarvalue::type_scalar: - return (this->real_value - x2.real_value)*(this->real_value - x2.real_value); - case colvarvalue::type_3vector: - return (this->rvector_value - x2.rvector_value).norm2(); - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - // angle between (*this) and x2 is the distance - return cvm::acos(this->rvector_value * x2.rvector_value) * cvm::acos(this->rvector_value * x2.rvector_value); - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - // angle between (*this) and x2 is the distance, the quaternion - // object has it implemented internally - return this->quaternion_value.dist2(x2.quaternion_value); - case colvarvalue::type_vector: - return (this->vector1d_value - x2.vector1d_value).norm2(); - case colvarvalue::type_notset: - default: - this->undef_op(); - return 0.0; - }; -} - - #endif