1#ifndef GEOMETRICPATHCV_H
2#define GEOMETRICPATHCV_H
19namespace GeometricPathCV {
23template <
typename element_type,
typename scalar_type, path_sz path_type>
29 bool operator()(
const size_t& i1,
const size_t& i2) {
30 return m_obj.frame_distances[i1] < m_obj.frame_distances[i2];
45 std::vector<element_type> v1;
46 std::vector<element_type> v2;
47 std::vector<element_type> v3;
48 std::vector<element_type> v4;
49 std::vector<element_type> dfdv1;
50 std::vector<element_type> dfdv2;
51 std::vector<element_type> dzdv1;
52 std::vector<element_type> dzdv2;
53 std::vector<scalar_type> frame_distances;
54 std::vector<size_t> frame_index;
55 bool use_second_closest_frame;
56 bool use_third_closest_frame;
58 long min_frame_index_1;
59 long min_frame_index_2;
60 long min_frame_index_3;
65 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);
66 GeometricPathBase(
size_t vector_size,
const std::vector<element_type>& elements,
size_t total_frames = 1,
bool p_use_second_closest_frame =
true,
bool p_use_third_closest_frame =
false,
bool p_use_z_square =
false);
69 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);
70 virtual void initialize(
size_t vector_size,
const std::vector<element_type>& elements,
size_t total_frames = 1,
bool p_use_second_closest_frame =
true,
bool p_use_third_closest_frame =
false,
bool p_use_z_square =
false);
71 virtual void prepareVectors() = 0;
72 virtual void updateDistanceToReferenceFrames() = 0;
73 virtual void compute();
74 virtual void determineClosestFrames();
75 virtual void computeValue();
76 virtual void computeDerivatives();
79template <
typename element_type,
typename scalar_type, path_sz path_type>
80GeometricPathBase<element_type, scalar_type, path_type>::GeometricPathBase(
size_t vector_size,
const element_type& element,
size_t total_frames,
bool p_use_second_closest_frame,
bool p_use_third_closest_frame,
bool p_use_z_square) {
81 initialize(vector_size, element, total_frames, p_use_second_closest_frame, p_use_third_closest_frame, p_use_z_square);
84template <
typename element_type,
typename scalar_type, path_sz path_type>
85GeometricPathBase<element_type, scalar_type, path_type>::GeometricPathBase(
size_t vector_size,
const std::vector<element_type>& elements,
size_t total_frames,
bool p_use_second_closest_frame,
bool p_use_third_closest_frame,
bool p_use_z_square) {
86 initialize(vector_size, elements, total_frames, p_use_second_closest_frame, p_use_third_closest_frame, p_use_z_square);
89template <
typename element_type,
typename scalar_type, path_sz path_type>
90void GeometricPathBase<element_type, scalar_type, path_type>::initialize(
size_t vector_size,
const element_type& element,
size_t total_frames,
bool p_use_second_closest_frame,
bool p_use_third_closest_frame,
bool p_use_z_square) {
102 v1.resize(vector_size, element);
103 v2.resize(vector_size, element);
104 v3.resize(vector_size, element);
105 v4.resize(vector_size, element);
106 dfdv1.resize(vector_size, element);
107 dfdv2.resize(vector_size, element);
108 dzdv1.resize(vector_size, element);
109 dzdv2.resize(vector_size, element);
110 frame_distances.resize(total_frames);
111 frame_index.resize(total_frames);
112 for (
size_t i_frame = 0; i_frame < frame_index.size(); ++i_frame) {
113 frame_index[i_frame] = i_frame;
115 use_second_closest_frame = p_use_second_closest_frame;
116 use_third_closest_frame = p_use_third_closest_frame;
117 use_z_square = p_use_z_square;
118 M =
static_cast<scalar_type
>(total_frames - 1);
119 m =
static_cast<scalar_type
>(1.0);
122template <
typename element_type,
typename scalar_type, path_sz path_type>
123void GeometricPathBase<element_type, scalar_type, path_type>::initialize(
size_t ,
const std::vector<element_type>& elements,
size_t total_frames,
bool p_use_second_closest_frame,
bool p_use_third_closest_frame,
bool p_use_z_square) {
124 v1v1 = scalar_type();
125 v2v2 = scalar_type();
126 v3v3 = scalar_type();
127 v4v4 = scalar_type();
128 v1v3 = scalar_type();
129 v1v4 = scalar_type();
143 frame_distances.resize(total_frames);
144 frame_index.resize(total_frames);
145 for (
size_t i_frame = 0; i_frame < frame_index.size(); ++i_frame) {
146 frame_index[i_frame] = i_frame;
148 use_second_closest_frame = p_use_second_closest_frame;
149 use_third_closest_frame = p_use_third_closest_frame;
150 use_z_square = p_use_z_square;
151 M =
static_cast<scalar_type
>(total_frames - 1);
152 m =
static_cast<scalar_type
>(1.0);
155template <
typename element_type,
typename scalar_type, path_sz path_type>
156void GeometricPathBase<element_type, scalar_type, path_type>::compute() {
158 computeDerivatives();
161template <
typename element_type,
typename scalar_type, path_sz path_type>
162void GeometricPathBase<element_type, scalar_type, path_type>::determineClosestFrames() {
164 std::sort(frame_index.begin(), frame_index.end(), doCompareFrameDistance(*
this));
166 sign =
static_cast<long>(frame_index[0]) -
static_cast<long>(frame_index[1]);
170 }
else if (sign < -1) {
174 if (
cvm::fabs(
static_cast<long>(frame_index[0]) -
static_cast<long>(frame_index[1])) > 1) {
176 "Warning: Geometrical pathCV relies on the assumption that the second closest frame is "
177 "the neighbouring frame\n"
178 " Please check your configuration or increase restraint on z(σ)\n");
179 for (
size_t i_frame = 0; i_frame < frame_index.size(); ++i_frame) {
180 message +=
"Frame index: " +
cvm::to_str(frame_index[i_frame]) +
181 " ; optimal RMSD = " +
cvm::to_str(frame_distances[frame_index[i_frame]]) +
185 min_frame_index_1 = frame_index[0];
186 min_frame_index_2 = use_second_closest_frame ? frame_index[1] : min_frame_index_1 - sign;
187 min_frame_index_3 = use_third_closest_frame ? frame_index[2] : min_frame_index_1 + sign;
188 m =
static_cast<double>(frame_index[0]);
191template <
typename element_type,
typename scalar_type, path_sz path_type>
192void GeometricPathBase<element_type, scalar_type, path_type>::computeValue() {
193 updateDistanceToReferenceFrames();
194 determineClosestFrames();
196 v1v1 = scalar_type();
197 v2v2 = scalar_type();
198 v3v3 = scalar_type();
199 v1v3 = scalar_type();
200 if (path_type == Z) {
201 v1v4 = scalar_type();
202 v4v4 = scalar_type();
204 for (
size_t i_elem = 0; i_elem < v1.size(); ++i_elem) {
205 v1v1 += v1[i_elem] * v1[i_elem];
206 v2v2 += v2[i_elem] * v2[i_elem];
207 v3v3 += v3[i_elem] * v3[i_elem];
208 v1v3 += v1[i_elem] * v3[i_elem];
209 if (path_type == Z) {
210 v1v4 += v1[i_elem] * v4[i_elem];
211 v4v4 += v4[i_elem] * v4[i_elem];
214 f = (
cvm::sqrt(v1v3 * v1v3 - v3v3 * (v1v1 - v2v2)) - v1v3) / v3v3;
215 if (path_type == Z) {
217 zz = v1v1 + 2 * dx * v1v4 + dx * dx * v4v4;
224 if (path_type == S) {
225 s = m/M +
static_cast<double>(sign) * ((f - 1) / (2 * M));
229template <
typename element_type,
typename scalar_type, path_sz path_type>
230void GeometricPathBase<element_type, scalar_type, path_type>::computeDerivatives() {
231 const scalar_type factor1 = 1.0 / (2.0 * v3v3 *
cvm::sqrt(v1v3 * v1v3 - v3v3 * (v1v1 - v2v2)));
232 const scalar_type factor2 = 1.0 / v3v3;
233 for (
size_t i_elem = 0; i_elem < v1.size(); ++i_elem) {
235 dfdv1[i_elem] = factor1 * (2.0 * v1v3 * v3[i_elem] - 2.0 * v3v3 * v1[i_elem]) - factor2 * v3[i_elem];
237 dfdv2[i_elem] = factor1 * (2.0 * v3v3 * v2[i_elem]);
243 if (path_type == Z) {
245 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];
246 dzdv2[i_elem] = v1v4 * dfdv2[i_elem] + v4v4 * 0.25 * 2.0 * (f-1) * dfdv2[i_elem];
248 if (z >
static_cast<scalar_type
>(0)) {
249 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]);
250 dzdv2[i_elem] = (1.0 / (2.0 * z)) * (v1v4 * dfdv2[i_elem] + v4v4 * 0.25 * 2.0 * (f-1) * dfdv2[i_elem]);
Definition: colvar_geometricpath.h:24
static real sqrt(real const &x)
Reimplemented to work around MS compiler issues.
Definition: colvarmodule.h:133
static real fabs(real const &x)
Reimplemented to work around MS compiler issues.
Definition: colvarmodule.h:127
static std::string to_str(char const *s)
Convert to string for output purposes.
Definition: colvarmodule.cpp:2389
Collective variables main module.
Definition: colvar_geometricpath.h:26