ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Update Colvars library to version 2024-06-04

Browse Source
This commit is contained in:
Giacomo Fiorin
2024-08-06 01:07:43 +02:00
parent 278accd9ea
commit 133dee9ac1
74 changed files with 7343 additions and 4676 deletions
Download Patch File Download Diff File Expand all files Collapse all files

442
lib/colvars/colvarcomp_apath.cpp
View File

@ -1,78 +1,312 @@
#if (__cplusplus >= 201103L)
// -*- c++ -*-
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/Colvars/colvars
// Please update all Colvars source files before making any changes.
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include <numeric>
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <numeric>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
#include "colvar_arithmeticpath.h"
colvar::aspathCV::aspathCV(std::string const &conf): CVBasedPath(conf) {
struct ArithmeticPathImpl: public ArithmeticPathCV::ArithmeticPathBase<cvm::real> {
std::vector<std::vector<colvarvalue>> frame_element_distances;
std::vector<std::vector<colvarvalue>> dsdx;
std::vector<std::vector<colvarvalue>> dzdx;
template <typename T>
void updateCartesianDistanceToReferenceFrames(T* obj) {
for (size_t i_frame = 0; i_frame < obj->reference_frames.size(); ++i_frame) {
for (size_t i_atom = 0; i_atom < obj->atoms->size(); ++i_atom) {
frame_element_distances[i_frame][i_atom] = (*(obj->comp_atoms[i_frame]))[i_atom].pos - obj->reference_frames[i_frame][i_atom];
}
}
}
template <typename T>
void updateCVDistanceToReferenceFrames(T* obj) {
for (size_t i_cv = 0; i_cv < obj->cv.size(); ++i_cv) {
obj->cv[i_cv]->calc_value();
}
for (size_t i_frame = 0; i_frame < obj->ref_cv.size(); ++i_frame) {
for (size_t i_cv = 0; i_cv < obj->cv.size(); ++i_cv) {
colvarvalue ref_cv_value(obj->ref_cv[i_frame][i_cv]);
colvarvalue current_cv_value(obj->cv[i_cv]->value());
if (current_cv_value.type() == colvarvalue::type_scalar) {
frame_element_distances[i_frame][i_cv] = 0.5 * obj->cv[i_cv]->dist2_lgrad(obj->cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, obj->cv[i_cv]->sup_np)), ref_cv_value.real_value);
} else {
frame_element_distances[i_frame][i_cv] = 0.5 * obj->cv[i_cv]->dist2_lgrad(obj->cv[i_cv]->sup_coeff * current_cv_value, ref_cv_value);
}
}
}
}
ArithmeticPathImpl(size_t p_num_elements, size_t p_total_frames, cvm::real p_lambda, const std::vector<cvm::real>& p_weights) {
ArithmeticPathCV::ArithmeticPathBase<cvm::real>::initialize(p_num_elements, p_total_frames, p_lambda, p_weights);
frame_element_distances.resize(p_total_frames, std::vector<colvarvalue>(p_num_elements, colvarvalue(colvarvalue::Type::type_notset)));
dsdx.resize(p_total_frames, std::vector<colvarvalue>(p_num_elements, colvarvalue(colvarvalue::Type::type_notset)));
dzdx.resize(p_total_frames, std::vector<colvarvalue>(p_num_elements, colvarvalue(colvarvalue::Type::type_notset)));
}
cvm::real get_lambda() const {return lambda;}
cvm::real compute_s() {
cvm::real s;
computeValue(frame_element_distances, &s, nullptr);
return s;
}
cvm::real compute_z() {
cvm::real z;
computeValue(frame_element_distances, nullptr, &z);
return z;
}
void compute_s_derivatives() {
computeDerivatives<colvarvalue>(frame_element_distances, &dsdx, nullptr);
}
void compute_z_derivatives() {
computeDerivatives<colvarvalue>(frame_element_distances, nullptr, &dzdx);
}
// for debug gradients of implicit sub CVs
template <typename T>
colvarvalue compute_s_analytical_derivative_ij(size_t i, size_t j, cvm::real eps, T* obj) const {
ArithmeticPathImpl tmp_left(*this), tmp_right(*this);
const size_t value_size = frame_element_distances[i][j].size();
colvarvalue result(frame_element_distances[i][j].type());
colvarvalue ref_cv_value(obj->ref_cv[i][j]);
for (size_t k = 0; k < value_size; ++k) {
// get the current CV value
colvarvalue current_cv_value(obj->cv[j]->value());
// save the old values in frame element distance matrices
const auto saved_left = tmp_left.frame_element_distances[i][j][k];
const auto saved_right = tmp_right.frame_element_distances[i][j][k];
// update frame element distance matrices
if (current_cv_value.type() == colvarvalue::type_scalar) {
tmp_left.frame_element_distances[i][j] = 0.5 * obj->cv[j]->dist2_lgrad(obj->cv[j]->sup_coeff * (cvm::pow(current_cv_value.real_value - eps, obj->cv[j]->sup_np)), ref_cv_value.real_value);
tmp_right.frame_element_distances[i][j] = 0.5 * obj->cv[j]->dist2_lgrad(obj->cv[j]->sup_coeff * (cvm::pow(current_cv_value.real_value + eps, obj->cv[j]->sup_np)), ref_cv_value.real_value);
} else {
current_cv_value[k] -= eps;
tmp_left.frame_element_distances[i][j] = 0.5 * obj->cv[j]->dist2_lgrad(obj->cv[j]->sup_coeff * current_cv_value, ref_cv_value);
current_cv_value[k] += eps + eps;
tmp_right.frame_element_distances[i][j] = 0.5 * obj->cv[j]->dist2_lgrad(obj->cv[j]->sup_coeff * current_cv_value, ref_cv_value);
}
const cvm::real s_left = tmp_left.compute_s();
const cvm::real s_right = tmp_right.compute_s();
result[k] = (s_right - s_left) / (2.0 * eps);
tmp_left.frame_element_distances[i][j][k] = saved_left;
tmp_right.frame_element_distances[i][j][k] = saved_right;
}
return result;
}
template <typename T>
colvarvalue compute_z_analytical_derivative_ij(size_t i, size_t j, cvm::real eps, T* obj) const {
ArithmeticPathImpl tmp_left(*this), tmp_right(*this);
const size_t value_size = frame_element_distances[i][j].size();
colvarvalue result(frame_element_distances[i][j].type());
colvarvalue ref_cv_value(obj->ref_cv[i][j]);
for (size_t k = 0; k < value_size; ++k) {
// get the current CV value
colvarvalue current_cv_value(obj->cv[j]->value());
// save the old values in frame element distance matrices
const auto saved_left = tmp_left.frame_element_distances[i][j][k];
const auto saved_right = tmp_right.frame_element_distances[i][j][k];
// update frame element distance matrices
if (current_cv_value.type() == colvarvalue::type_scalar) {
tmp_left.frame_element_distances[i][j] = 0.5 * obj->cv[j]->dist2_lgrad(obj->cv[j]->sup_coeff * (cvm::pow(current_cv_value.real_value - eps, obj->cv[j]->sup_np)), ref_cv_value.real_value);
tmp_right.frame_element_distances[i][j] = 0.5 * obj->cv[j]->dist2_lgrad(obj->cv[j]->sup_coeff * (cvm::pow(current_cv_value.real_value + eps, obj->cv[j]->sup_np)), ref_cv_value.real_value);
} else {
current_cv_value[k] -= eps;
tmp_left.frame_element_distances[i][j] = 0.5 * obj->cv[j]->dist2_lgrad(obj->cv[j]->sup_coeff * current_cv_value, ref_cv_value);
current_cv_value[k] += eps + eps;
tmp_right.frame_element_distances[i][j] = 0.5 * obj->cv[j]->dist2_lgrad(obj->cv[j]->sup_coeff * current_cv_value, ref_cv_value);
}
const cvm::real z_left = tmp_left.compute_z();
const cvm::real z_right = tmp_right.compute_z();
result[k] = (z_right - z_left) / (2.0 * eps);
tmp_left.frame_element_distances[i][j][k] = saved_left;
tmp_right.frame_element_distances[i][j][k] = saved_right;
}
return result;
}
};
colvar::aspath::aspath()
{
set_function_type("aspath");
x.type(colvarvalue::type_scalar);
}
int colvar::aspath::init(std::string const &conf)
{
int error_code = CartesianBasedPath::init(conf);
if (error_code != COLVARS_OK) return error_code;
cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n"));
cvm::real p_lambda;
get_keyval(conf, "lambda", p_lambda, -1.0);
const size_t num_atoms = atoms->size();
std::vector<cvm::real> p_weights(num_atoms, std::sqrt(1.0 / num_atoms));
// ArithmeticPathCV::ArithmeticPathBase<cvm::atom_pos, cvm::real, ArithmeticPathCV::path_sz::S>::initialize(num_atoms, total_reference_frames, p_lambda, reference_frames[0], p_weights);
if (impl_) impl_.reset();
impl_ = std::unique_ptr<ArithmeticPathImpl>(new ArithmeticPathImpl(num_atoms, total_reference_frames, p_lambda, p_weights));
cvm::log(std::string("Lambda is ") + cvm::to_str(impl_->get_lambda()) + std::string("\n"));
return error_code;
}
colvar::aspath::~aspath() {}
void colvar::aspath::calc_value() {
if (impl_->get_lambda() < 0) {
// this implies that the user may not set a valid lambda value
// so recompute it by the suggested value in Parrinello's paper
cvm::log("A non-positive value of lambda is detected, which implies that it may not set in the configuration.\n");
cvm::log("This component (aspath) will recompute a value for lambda following the suggestion in the origin paper.\n");
std::vector<cvm::real> rmsd_between_refs(total_reference_frames - 1, 0.0);
computeDistanceBetweenReferenceFrames(rmsd_between_refs);
impl_->reComputeLambda(rmsd_between_refs);
cvm::log("Ok, the value of lambda is updated to " + cvm::to_str(impl_->get_lambda()));
}
impl_->updateCartesianDistanceToReferenceFrames(this);
x = impl_->compute_s();
}
void colvar::aspath::calc_gradients() {
impl_->compute_s_derivatives();
for (size_t i_frame = 0; i_frame < reference_frames.size(); ++i_frame) {
for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
(*(comp_atoms[i_frame]))[i_atom].grad += impl_->dsdx[i_frame][i_atom];
}
}
}
void colvar::aspath::apply_force(colvarvalue const &force) {
cvm::real const &F = force.real_value;
for (size_t i_frame = 0; i_frame < reference_frames.size(); ++i_frame) {
(*(comp_atoms[i_frame])).apply_colvar_force(F);
}
}
colvar::azpath::azpath()
{
set_function_type("azpath");
x.type(colvarvalue::type_scalar);
}
int colvar::azpath::init(std::string const &conf)
{
int error_code = CartesianBasedPath::init(conf);
if (error_code != COLVARS_OK) return error_code;
cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n"));
x.type(colvarvalue::type_scalar);
cvm::real p_lambda;
get_keyval(conf, "lambda", p_lambda, -1.0);
const size_t num_atoms = atoms->size();
std::vector<cvm::real> p_weights(num_atoms, std::sqrt(1.0 / num_atoms));
if (impl_) impl_.reset();
impl_ = std::unique_ptr<ArithmeticPathImpl>(new ArithmeticPathImpl(num_atoms, total_reference_frames, p_lambda, p_weights));
cvm::log(std::string("Lambda is ") + cvm::to_str(impl_->get_lambda()) + std::string("\n"));
return error_code;
}
colvar::azpath::~azpath() {}
void colvar::azpath::calc_value() {
if (impl_->get_lambda() < 0) {
// this implies that the user may not set a valid lambda value
// so recompute it by the suggested value in Parrinello's paper
cvm::log("A non-positive value of lambda is detected, which implies that it may not set in the configuration.\n");
cvm::log("This component (azpath) will recompute a value for lambda following the suggestion in the origin paper.\n");
std::vector<cvm::real> rmsd_between_refs(total_reference_frames - 1, 0.0);
computeDistanceBetweenReferenceFrames(rmsd_between_refs);
impl_->reComputeLambda(rmsd_between_refs);
cvm::log("Ok, the value of lambda is updated to " + cvm::to_str(impl_->get_lambda()));
}
impl_->updateCartesianDistanceToReferenceFrames(this);
x = impl_->compute_z();
}
void colvar::azpath::calc_gradients() {
impl_->compute_z_derivatives();
for (size_t i_frame = 0; i_frame < reference_frames.size(); ++i_frame) {
for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
(*(comp_atoms[i_frame]))[i_atom].grad += impl_->dzdx[i_frame][i_atom];
}
}
}
void colvar::azpath::apply_force(colvarvalue const &force) {
cvm::real const &F = force.real_value;
for (size_t i_frame = 0; i_frame < reference_frames.size(); ++i_frame) {
(*(comp_atoms[i_frame])).apply_colvar_force(F);
}
}
colvar::aspathCV::aspathCV()
{
set_function_type("aspathCV");
x.type(colvarvalue::type_scalar);
}
int colvar::aspathCV::init(std::string const &conf)
{
int error_code = CVBasedPath::init(conf);
if (error_code != COLVARS_OK) return error_code;
cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n"));
std::vector<cvm::real> p_weights(cv.size(), 1.0);
get_keyval(conf, "weights", p_weights, std::vector<cvm::real>(cv.size(), 1.0));
x.type(colvarvalue::type_scalar);
use_explicit_gradients = true;
cvm::real p_lambda;
get_keyval(conf, "lambda", p_lambda, -1.0);
ArithmeticPathCV::ArithmeticPathBase<colvarvalue, cvm::real, ArithmeticPathCV::path_sz::S>::initialize(cv.size(), total_reference_frames, p_lambda, ref_cv[0], p_weights);
cvm::log(std::string("Lambda is ") + cvm::to_str(lambda) + std::string("\n"));
if (impl_) impl_.reset();
impl_ = std::unique_ptr<ArithmeticPathImpl>(new ArithmeticPathImpl(cv.size(), total_reference_frames, p_lambda, p_weights));
cvm::log(std::string("Lambda is ") + cvm::to_str(impl_->get_lambda()) + std::string("\n"));
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;
}
cvm::log(std::string("The weight of CV ") + cvm::to_str(i_cv) + std::string(" is ") + cvm::to_str(weights[i_cv]) + std::string("\n"));
cvm::log(std::string("The weight of CV ") + cvm::to_str(i_cv) + std::string(" is ") + cvm::to_str(p_weights[i_cv]) + std::string("\n"));
}
return error_code;
}
void colvar::aspathCV::updateDistanceToReferenceFrames() {
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) {
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());
if (current_cv_value.type() == colvarvalue::type_scalar) {
frame_element_distances[i_frame][i_cv] = 0.5 * cv[i_cv]->dist2_lgrad(cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)), ref_cv_value.real_value);
} else {
frame_element_distances[i_frame][i_cv] = 0.5 * cv[i_cv]->dist2_lgrad(cv[i_cv]->sup_coeff * current_cv_value, ref_cv_value);
}
}
}
}
colvar::aspathCV::~aspathCV() {}
void colvar::aspathCV::calc_value() {
if (lambda < 0) {
if (impl_->get_lambda() < 0) {
// this implies that the user may not set a valid lambda value
// so recompute it by the suggested value in Parrinello's paper
cvm::log("A non-positive value of lambda is detected, which implies that it may not set in the configuration.\n");
cvm::log("This component (aspathCV) will recompute a value for lambda following the suggestion in the origin paper.\n");
std::vector<cvm::real> rmsd_between_refs(total_reference_frames - 1, 0.0);
computeDistanceBetweenReferenceFrames(rmsd_between_refs);
reComputeLambda(rmsd_between_refs);
cvm::log("Ok, the value of lambda is updated to " + cvm::to_str(lambda));
impl_->reComputeLambda(rmsd_between_refs);
cvm::log("Ok, the value of lambda is updated to " + cvm::to_str(impl_->get_lambda()));
}
computeValue();
x = s;
impl_->updateCVDistanceToReferenceFrames(this);
x = impl_->compute_s();
}
void colvar::aspathCV::calc_gradients() {
computeDerivatives();
impl_->compute_s_derivatives();
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)) {
cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
// compute the gradient (grad) with respect to the i-th CV
colvarvalue grad(cv[i_cv]->value().type());
// sum up derivatives with respect to all frames
for (size_t m_frame = 0; m_frame < impl_->dsdx.size(); ++m_frame) {
// dsdx is the derivative of s with respect to the m-th frame
grad += impl_->dsdx[m_frame][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 = dsdx[i_cv][j_elem] * factor_polynomial * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad;
(*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = grad[j_elem] * factor_polynomial * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad;
}
}
}
@ -88,97 +322,137 @@ void colvar::aspathCV::apply_force(colvarvalue const &force) {
}
} else {
cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
colvarvalue cv_force = dsdx[i_cv] * force.real_value * factor_polynomial;
cv[i_cv]->apply_force(cv_force);
}
}
}
colvar::aspathCV::~aspathCV() {}
colvar::azpathCV::azpathCV(std::string const &conf): CVBasedPath(conf) {
set_function_type("azpathCV");
cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n"));
std::vector<cvm::real> p_weights(cv.size(), 1.0);
get_keyval(conf, "weights", p_weights, std::vector<cvm::real>(cv.size(), 1.0));
x.type(colvarvalue::type_scalar);
use_explicit_gradients = true;
cvm::real p_lambda;
get_keyval(conf, "lambda", p_lambda, -1.0);
ArithmeticPathCV::ArithmeticPathBase<colvarvalue, cvm::real, ArithmeticPathCV::path_sz::Z>::initialize(cv.size(), total_reference_frames, p_lambda, ref_cv[0], p_weights);
cvm::log(std::string("Lambda is ") + cvm::to_str(lambda) + std::string("\n"));
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;
}
cvm::log(std::string("The weight of CV ") + cvm::to_str(i_cv) + std::string(" is ") + cvm::to_str(weights[i_cv]) + std::string("\n"));
}
}
void colvar::azpathCV::updateDistanceToReferenceFrames() {
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) {
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());
if (current_cv_value.type() == colvarvalue::type_scalar) {
frame_element_distances[i_frame][i_cv] = 0.5 * cv[i_cv]->dist2_lgrad(cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)), ref_cv_value.real_value);
} else {
frame_element_distances[i_frame][i_cv] = 0.5 * cv[i_cv]->dist2_lgrad(cv[i_cv]->sup_coeff * current_cv_value, ref_cv_value);
// compute the gradient (grad) with respect to the i-th CV
colvarvalue grad(cv[i_cv]->value().type());
for (size_t m_frame = 0; m_frame < impl_->dsdx.size(); ++m_frame) {
// dsdx is the derivative of s with respect to the m-th frame
grad += impl_->dsdx[m_frame][i_cv];
}
grad *= factor_polynomial;
cv[i_cv]->apply_force(force.real_value * grad);
// try my best to debug gradients even if the sub-CVs do not have explicit gradients
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging gradients for " + description +
" with respect to sub-CV " + cv[i_cv]->description +
", which has no explicit gradient with respect to its own input(s)");
colvarvalue analytical_grad(cv[i_cv]->value().type());
for (size_t m_frame = 0; m_frame < impl_->dsdx.size(); ++m_frame) {
analytical_grad += impl_->compute_s_analytical_derivative_ij(
m_frame, i_cv, cvm::debug_gradients_step_size, this);
}
cvm::log("dx(actual) = "+cvm::to_str(analytical_grad, 21, 14)+"\n");
cvm::log("dx(interp) = "+cvm::to_str(grad, 21, 14)+"\n");
cvm::log("|dx(actual) - dx(interp)|/|dx(actual)| = "+
cvm::to_str((analytical_grad - grad).norm() /
(analytical_grad).norm(), 12, 5)+"\n");
}
}
}
}
colvar::azpathCV::azpathCV()
{
set_function_type("azpathCV");
x.type(colvarvalue::type_scalar);
}
int colvar::azpathCV::init(std::string const &conf)
{
int error_code = CVBasedPath::init(conf);
if (error_code != COLVARS_OK) return error_code;
cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n"));
std::vector<cvm::real> p_weights(cv.size(), 1.0);
get_keyval(conf, "weights", p_weights, std::vector<cvm::real>(cv.size(), 1.0));
use_explicit_gradients = true;
cvm::real p_lambda;
get_keyval(conf, "lambda", p_lambda, -1.0);
if (impl_) impl_.reset();
impl_ = std::unique_ptr<ArithmeticPathImpl>(new ArithmeticPathImpl(cv.size(), total_reference_frames, p_lambda, p_weights));
cvm::log(std::string("Lambda is ") + cvm::to_str(impl_->get_lambda()) + std::string("\n"));
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;
}
cvm::log(std::string("The weight of CV ") + cvm::to_str(i_cv) + std::string(" is ") + cvm::to_str(p_weights[i_cv]) + std::string("\n"));
}
return error_code;
}
void colvar::azpathCV::calc_value() {
if (lambda < 0) {
if (impl_->get_lambda() < 0) {
// this implies that the user may not set a valid lambda value
// so recompute it by the suggested value in Parrinello's paper
cvm::log("A non-positive value of lambda is detected, which implies that it may not set in the configuration.\n");
cvm::log("This component (azpathCV) will recompute a value for lambda following the suggestion in the origin paper.\n");
std::vector<cvm::real> rmsd_between_refs(total_reference_frames - 1, 0.0);
computeDistanceBetweenReferenceFrames(rmsd_between_refs);
reComputeLambda(rmsd_between_refs);
cvm::log("Ok, the value of lambda is updated to " + cvm::to_str(lambda));
impl_->reComputeLambda(rmsd_between_refs);
cvm::log("Ok, the value of lambda is updated to " + cvm::to_str(impl_->get_lambda()));
}
computeValue();
x = z;
impl_->updateCVDistanceToReferenceFrames(this);
x = impl_->compute_z();
}
void colvar::azpathCV::calc_gradients() {
computeDerivatives();
impl_->compute_z_derivatives();
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)) {
cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
// compute the gradient (grad) with respect to the i-th CV
colvarvalue grad(cv[i_cv]->value().type());
// sum up derivatives with respect to all frames
for (size_t m_frame = 0; m_frame < impl_->dzdx.size(); ++m_frame) {
// dzdx is the derivative of z with respect to the m-th frame
grad += impl_->dzdx[m_frame][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 = dzdx[i_cv][j_elem] * factor_polynomial * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad;
(*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = grad[j_elem] * factor_polynomial * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad;
}
}
}
}
}
}
void colvar::azpathCV::apply_force(colvarvalue const &force) {
// the PCV component itself is a scalar, so force should be scalar
for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
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 {
cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
const colvarvalue cv_force = dzdx[i_cv] * force.real_value * factor_polynomial;
cv[i_cv]->apply_force(cv_force);
const cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
// compute the gradient (grad) with respect to the i-th CV
colvarvalue grad(cv[i_cv]->value().type());
for (size_t m_frame = 0; m_frame < impl_->dzdx.size(); ++m_frame) {
// dzdx is the derivative of z with respect to the m-th frame
grad += impl_->dzdx[m_frame][i_cv];
}
grad *= factor_polynomial;
cv[i_cv]->apply_force(force.real_value * grad);
// try my best to debug gradients even if the sub-CVs do not have explicit gradients
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging gradients for " + description +
" with respect to sub-CV " + cv[i_cv]->description +
", which has no explicit gradient with respect to its own input(s)");
colvarvalue analytical_grad(cv[i_cv]->value().type());
for (size_t m_frame = 0; m_frame < impl_->dzdx.size(); ++m_frame) {
analytical_grad += impl_->compute_z_analytical_derivative_ij(
m_frame, i_cv, cvm::debug_gradients_step_size, this);
}
cvm::log("dx(actual) = "+cvm::to_str(analytical_grad, 21, 14)+"\n");
cvm::log("dx(interp) = "+cvm::to_str(grad, 21, 14)+"\n");
cvm::log("|dx(actual) - dx(interp)|/|dx(actual)| = "+
cvm::to_str((analytical_grad - grad).norm() /
(analytical_grad).norm(), 12, 5)+"\n");
}
}
}
}
colvar::azpathCV::~azpathCV() {}
#endif