Collective Variables Module - Developer Documentation
Loading...
Searching...
No Matches
colvar_geometricpath.h
1#ifndef GEOMETRICPATHCV_H
2#define GEOMETRICPATHCV_H
3// This file is part of the Collective Variables module (Colvars).
4// The original version of Colvars and its updates are located at:
5// https://github.com/Colvars/colvars
6// Please update all Colvars source files before making any changes.
7// If you wish to distribute your changes, please submit them to the
8// Colvars repository at GitHub.
9
10
11#include <algorithm>
12#include <cmath>
13#include <string>
14#include <vector>
15
16#include "colvarmodule.h"
17
18
19namespace GeometricPathCV {
20
21enum path_sz {S, Z};
22
23template <typename element_type, typename scalar_type, path_sz path_type>
24class GeometricPathBase {
25private:
26 struct doCompareFrameDistance {
27 doCompareFrameDistance(const GeometricPathBase& obj): m_obj(obj) {}
28 const GeometricPathBase& m_obj;
29 bool operator()(const size_t& i1, const size_t& i2) {
30 return m_obj.frame_distances[i1] < m_obj.frame_distances[i2];
31 }
32 };
33protected:
34 scalar_type v1v1;
35 scalar_type v2v2;
36 scalar_type v3v3;
37 scalar_type v4v4;
38 scalar_type v1v3;
39 scalar_type v1v4;
40 scalar_type f;
41 scalar_type dx;
42 scalar_type s;
43 scalar_type z;
44 scalar_type zz;
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;
57 bool use_z_square;
58 long min_frame_index_1;
59 long min_frame_index_2;
60 long min_frame_index_3;
61 long sign;
62 double M;
63 double m;
64public:
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);
67 GeometricPathBase() {}
68 virtual ~GeometricPathBase() {}
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();
77};
78
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);
82}
83
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);
87}
88
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) {
91 v1v1 = scalar_type();
92 v2v2 = scalar_type();
93 v3v3 = scalar_type();
94 v4v4 = scalar_type();
95 v1v3 = scalar_type();
96 v1v4 = scalar_type();
97 f = scalar_type();
98 dx = scalar_type();
99 z = scalar_type();
100 zz = scalar_type();
101 sign = 0;
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;
114 }
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);
120}
121
122template <typename element_type, typename scalar_type, path_sz path_type>
123void GeometricPathBase<element_type, scalar_type, path_type>::initialize(size_t /* vector_size */, const std::vector<element_type>& elements, size_t total_frames, bool p_use_second_closest_frame, bool p_use_third_closest_frame, bool p_use_z_square) {
124 v1v1 = scalar_type();
125 v2v2 = scalar_type();
126 v3v3 = scalar_type();
127 v4v4 = scalar_type();
128 v1v3 = scalar_type();
129 v1v4 = scalar_type();
130 f = scalar_type();
131 dx = scalar_type();
132 z = scalar_type();
133 zz = scalar_type();
134 sign = 0;
135 v1 = elements;
136 v2 = elements;
137 v3 = elements;
138 v4 = elements;
139 dfdv1 = elements;
140 dfdv2 = elements;
141 dzdv1 = elements;
142 dzdv2 = elements;
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;
147 }
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);
153}
154
155template <typename element_type, typename scalar_type, path_sz path_type>
156void GeometricPathBase<element_type, scalar_type, path_type>::compute() {
157 computeValue();
158 computeDerivatives();
159}
160
161template <typename element_type, typename scalar_type, path_sz path_type>
162void GeometricPathBase<element_type, scalar_type, path_type>::determineClosestFrames() {
163 // Find the closest and the second closest frames
164 std::sort(frame_index.begin(), frame_index.end(), doCompareFrameDistance(*this));
165 // Determine the sign
166 sign = static_cast<long>(frame_index[0]) - static_cast<long>(frame_index[1]);
167 if (sign > 1) {
168 // sigma(z) is on the left side of the closest frame
169 sign = 1;
170 } else if (sign < -1) {
171 // sigma(z) is on the right side of the closest frame
172 sign = -1;
173 }
174 if (cvm::fabs(static_cast<long>(frame_index[0]) - static_cast<long>(frame_index[1])) > 1) {
175 std::string message(
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]]) +
182 "\n";
183 }
184 }
185 min_frame_index_1 = frame_index[0]; // s_m
186 min_frame_index_2 = use_second_closest_frame ? frame_index[1] : min_frame_index_1 - sign; // s_(m-1)
187 min_frame_index_3 = use_third_closest_frame ? frame_index[2] : min_frame_index_1 + sign; // s_(m+1)
188 m = static_cast<double>(frame_index[0]);
189}
190
191template <typename element_type, typename scalar_type, path_sz path_type>
192void GeometricPathBase<element_type, scalar_type, path_type>::computeValue() {
193 updateDistanceToReferenceFrames();
194 determineClosestFrames();
195 prepareVectors();
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();
203 }
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];
212 }
213 }
214 f = (cvm::sqrt(v1v3 * v1v3 - v3v3 * (v1v1 - v2v2)) - v1v3) / v3v3;
215 if (path_type == Z) {
216 dx = 0.5 * (f - 1);
217 zz = v1v1 + 2 * dx * v1v4 + dx * dx * v4v4;
218 if (use_z_square) {
219 z = zz;
220 } else {
221 z = cvm::sqrt(cvm::fabs(zz));
222 }
223 }
224 if (path_type == S) {
225 s = m/M + static_cast<double>(sign) * ((f - 1) / (2 * M));
226 }
227}
228
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) {
234 // Compute the derivative of f with vector v1
235 dfdv1[i_elem] = factor1 * (2.0 * v1v3 * v3[i_elem] - 2.0 * v3v3 * v1[i_elem]) - factor2 * v3[i_elem];
236 // Compute the derivative of f with respect to vector v2
237 dfdv2[i_elem] = factor1 * (2.0 * v3v3 * v2[i_elem]);
238 // dZ(v1(r), v2(r), v3) / dr = ∂Z/∂v1 * dv1/dr + ∂Z/∂v2 * dv2/dr
239 // dv1/dr = [fitting matrix 1][-1, ..., -1]
240 // dv2/dr = [fitting matrix 2][1, ..., 1]
241 // ∂Z/∂v1 = 1/(2*z) * (2v1 + (f-1)v4 + (v1⋅v4)∂f/∂v1 + v4^2 * 1/4 * 2(f-1) * ∂f/∂v1)
242 // ∂Z/∂v2 = 1/(2*z) * ((v1⋅v4)∂f/∂v2 + v4^2 * 1/4 * 2(f-1) * ∂f/∂v2)
243 if (path_type == Z) {
244 if (use_z_square) {
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];
247 } else {
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]);
251 } else {
252 // workaround at z = 0
253 dzdv1[i_elem] = 0;
254 dzdv2[i_elem] = 0;
255 }
256 }
257 }
258 }
259}
260
261}
262
263#endif // GEOMETRICPATHCV_H
GeometricPathCV::GeometricPathBase
Definition: colvar_geometricpath.h:24
colvarmodule::sqrt
static real sqrt(real const &x)
Reimplemented to work around MS compiler issues.
Definition: colvarmodule.h:133
colvarmodule::fabs
static real fabs(real const &x)
Reimplemented to work around MS compiler issues.
Definition: colvarmodule.h:127
colvarmodule::to_str
static std::string to_str(char const *s)
Convert to string for output purposes.
Definition: colvarmodule.cpp:2392
colvarmodule.h
Collective variables main module.
Generated on Mon Nov 18 2024 18:38:26 for Collective Variables Module - Developer Documentation by doxygen 1.9.5