colvar_geometricpath.h Source File

00001 #ifndef GEOMETRICPATHCV_H
00002 #define GEOMETRICPATHCV_H
00003 // This file is part of the Collective Variables module (Colvars).
00004 // The original version of Colvars and its updates are located at:
00005 // https://github.com/Colvars/colvars
00006 // Please update all Colvars source files before making any changes.
00007 // If you wish to distribute your changes, please submit them to the
00008 // Colvars repository at GitHub.
00009 
00010 
00011 #include "colvarmodule.h"
00012 
00013 #include <vector>
00014 #include <cmath>
00015 #include <algorithm>
00016 #include <string>
00017 
00018 namespace GeometricPathCV {
00019 
00020 enum path_sz {S, Z};
00021 
00022 template <typename element_type, typename scalar_type, path_sz path_type>
00023 class GeometricPathBase {
00024 private:
00025     struct doCompareFrameDistance {
00026         doCompareFrameDistance(const GeometricPathBase& obj): m_obj(obj) {}
00027         const GeometricPathBase& m_obj;
00028         bool operator()(const size_t& i1, const size_t& i2) {
00029             return m_obj.frame_distances[i1] < m_obj.frame_distances[i2];
00030         }
00031     };
00032 protected:
00033     scalar_type v1v1;
00034     scalar_type v2v2;
00035     scalar_type v3v3;
00036     scalar_type v4v4;
00037     scalar_type v1v3;
00038     scalar_type v1v4;
00039     scalar_type f;
00040     scalar_type dx;
00041     scalar_type s;
00042     scalar_type z;
00043     scalar_type zz;
00044     std::vector<element_type> v1;
00045     std::vector<element_type> v2;
00046     std::vector<element_type> v3;
00047     std::vector<element_type> v4;
00048     std::vector<element_type> dfdv1;
00049     std::vector<element_type> dfdv2;
00050     std::vector<element_type> dzdv1;
00051     std::vector<element_type> dzdv2;
00052     std::vector<scalar_type> frame_distances;
00053     std::vector<size_t> frame_index;
00054     bool use_second_closest_frame;
00055     bool use_third_closest_frame;
00056     bool use_z_square;
00057     long min_frame_index_1;
00058     long min_frame_index_2;
00059     long min_frame_index_3;
00060     long sign;
00061     double M;
00062     double m;
00063 public:
00064     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);
00065     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);
00066     GeometricPathBase() {}
00067     virtual ~GeometricPathBase() {}
00068     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);
00069     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);
00070     virtual void prepareVectors() = 0;
00071     virtual void updateDistanceToReferenceFrames() = 0;
00072     virtual void compute();
00073     virtual void determineClosestFrames();
00074     virtual void computeValue();
00075     virtual void computeDerivatives();
00076 };
00077 
00078 template <typename element_type, typename scalar_type, path_sz path_type>
00079 GeometricPathBase<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) {
00080     initialize(vector_size, element, total_frames, p_use_second_closest_frame, p_use_third_closest_frame, p_use_z_square);
00081 }
00082 
00083 template <typename element_type, typename scalar_type, path_sz path_type>
00084 GeometricPathBase<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) {
00085     initialize(vector_size, elements, total_frames, p_use_second_closest_frame, p_use_third_closest_frame, p_use_z_square);
00086 }
00087 
00088 template <typename element_type, typename scalar_type, path_sz path_type>
00089 void 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) {
00090     v1v1 = scalar_type();
00091     v2v2 = scalar_type();
00092     v3v3 = scalar_type();
00093     v4v4 = scalar_type();
00094     v1v3 = scalar_type();
00095     v1v4 = scalar_type();
00096     f = scalar_type();
00097     dx = scalar_type();
00098     z = scalar_type();
00099     zz = scalar_type();
00100     sign = 0;
00101     v1.resize(vector_size, element);
00102     v2.resize(vector_size, element);
00103     v3.resize(vector_size, element);
00104     v4.resize(vector_size, element);
00105     dfdv1.resize(vector_size, element);
00106     dfdv2.resize(vector_size, element);
00107     dzdv1.resize(vector_size, element);
00108     dzdv2.resize(vector_size, element);
00109     frame_distances.resize(total_frames);
00110     frame_index.resize(total_frames);
00111     for (size_t i_frame = 0; i_frame < frame_index.size(); ++i_frame) {
00112         frame_index[i_frame] = i_frame;
00113     }
00114     use_second_closest_frame = p_use_second_closest_frame;
00115     use_third_closest_frame = p_use_third_closest_frame;
00116     use_z_square = p_use_z_square;
00117     M = static_cast<scalar_type>(total_frames - 1);
00118     m = static_cast<scalar_type>(1.0);
00119 }
00120 
00121 template <typename element_type, typename scalar_type, path_sz path_type>
00122 void 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) {
00123     v1v1 = scalar_type();
00124     v2v2 = scalar_type();
00125     v3v3 = scalar_type();
00126     v4v4 = scalar_type();
00127     v1v3 = scalar_type();
00128     v1v4 = scalar_type();
00129     f = scalar_type();
00130     dx = scalar_type();
00131     z = scalar_type();
00132     zz = scalar_type();
00133     sign = 0;
00134     v1 = elements;
00135     v2 = elements;
00136     v3 = elements;
00137     v4 = elements;
00138     dfdv1 = elements;
00139     dfdv2 = elements;
00140     dzdv1 = elements;
00141     dzdv2 = elements;
00142     frame_distances.resize(total_frames);
00143     frame_index.resize(total_frames);
00144     for (size_t i_frame = 0; i_frame < frame_index.size(); ++i_frame) {
00145         frame_index[i_frame] = i_frame;
00146     }
00147     use_second_closest_frame = p_use_second_closest_frame;
00148     use_third_closest_frame = p_use_third_closest_frame;
00149     use_z_square = p_use_z_square;
00150     M = static_cast<scalar_type>(total_frames - 1);
00151     m = static_cast<scalar_type>(1.0);
00152 }
00153 
00154 template <typename element_type, typename scalar_type, path_sz path_type>
00155 void GeometricPathBase<element_type, scalar_type, path_type>::compute() {
00156     computeValue();
00157     computeDerivatives();
00158 }
00159 
00160 template <typename element_type, typename scalar_type, path_sz path_type>
00161 void GeometricPathBase<element_type, scalar_type, path_type>::determineClosestFrames() {
00162     // Find the closest and the second closest frames
00163     std::sort(frame_index.begin(), frame_index.end(), doCompareFrameDistance(*this));
00164     // Determine the sign
00165     sign = static_cast<long>(frame_index[0]) - static_cast<long>(frame_index[1]);
00166     if (sign > 1) {
00167         // sigma(z) is on the left side of the closest frame
00168         sign = 1;
00169     } else if (sign < -1) {
00170         // sigma(z) is on the right side of the closest frame
00171         sign = -1;
00172     }
00173     if (cvm::fabs(static_cast<long>(frame_index[0]) - static_cast<long>(frame_index[1])) > 1) {
00174         std::cout << "Warning: Geometrical pathCV relies on the assumption that the second closest frame is the neighbouring frame\n";
00175         std::cout << "         Please check your configuration or increase restraint on z(σ)\n";
00176         for (size_t i_frame = 0; i_frame < frame_index.size(); ++i_frame) {
00177             std::cout << "Frame index: " << frame_index[i_frame] << " ; optimal RMSD = " << frame_distances[frame_index[i_frame]] << "\n";
00178         }
00179     }
00180     min_frame_index_1 = frame_index[0];                                                         // s_m
00181     min_frame_index_2 = use_second_closest_frame ? frame_index[1] : min_frame_index_1 - sign;   // s_(m-1)
00182     min_frame_index_3 = use_third_closest_frame ? frame_index[2] : min_frame_index_1 + sign;    // s_(m+1)
00183     m = static_cast<double>(frame_index[0]);
00184 }
00185 
00186 template <typename element_type, typename scalar_type, path_sz path_type>
00187 void GeometricPathBase<element_type, scalar_type, path_type>::computeValue() {
00188     updateDistanceToReferenceFrames();
00189     determineClosestFrames();
00190     prepareVectors();
00191     v1v1 = scalar_type();
00192     v2v2 = scalar_type();
00193     v3v3 = scalar_type();
00194     v1v3 = scalar_type();
00195     if (path_type == Z) {
00196         v1v4 = scalar_type();
00197         v4v4 = scalar_type();
00198     }
00199     for (size_t i_elem = 0; i_elem < v1.size(); ++i_elem) {
00200         v1v1 += v1[i_elem] * v1[i_elem];
00201         v2v2 += v2[i_elem] * v2[i_elem];
00202         v3v3 += v3[i_elem] * v3[i_elem];
00203         v1v3 += v1[i_elem] * v3[i_elem];
00204         if (path_type == Z) {
00205             v1v4 += v1[i_elem] * v4[i_elem];
00206             v4v4 += v4[i_elem] * v4[i_elem];
00207         }
00208     }
00209     f = (cvm::sqrt(v1v3 * v1v3 - v3v3 * (v1v1 - v2v2)) - v1v3) / v3v3;
00210     if (path_type == Z) {
00211         dx = 0.5 * (f - 1);
00212         zz = v1v1 + 2 * dx * v1v4 + dx * dx * v4v4;
00213         if (use_z_square) {
00214             z = zz;
00215         } else {
00216             z = cvm::sqrt(cvm::fabs(zz));
00217         }
00218     }
00219     if (path_type == S) {
00220         s = m/M + static_cast<double>(sign) * ((f - 1) / (2 * M));
00221     }
00222 }
00223 
00224 template <typename element_type, typename scalar_type, path_sz path_type>
00225 void GeometricPathBase<element_type, scalar_type, path_type>::computeDerivatives() {
00226     const scalar_type factor1 = 1.0 / (2.0 * v3v3 * cvm::sqrt(v1v3 * v1v3 - v3v3 * (v1v1 - v2v2)));
00227     const scalar_type factor2 = 1.0 / v3v3;
00228     for (size_t i_elem = 0; i_elem < v1.size(); ++i_elem) {
00229         // Compute the derivative of f with vector v1
00230         dfdv1[i_elem] = factor1 * (2.0 * v1v3 * v3[i_elem] - 2.0 * v3v3 * v1[i_elem]) - factor2 * v3[i_elem];
00231         // Compute the derivative of f with respect to vector v2
00232         dfdv2[i_elem] = factor1 * (2.0 * v3v3 * v2[i_elem]);
00233         // dZ(v1(r), v2(r), v3) / dr = ∂Z/∂v1 * dv1/dr + ∂Z/∂v2 * dv2/dr
00234         // dv1/dr = [fitting matrix 1][-1, ..., -1]
00235         // dv2/dr = [fitting matrix 2][1, ..., 1]
00236         // ∂Z/∂v1 = 1/(2*z) * (2v1 + (f-1)v4 + (v1⋅v4)∂f/∂v1 + v4^2 * 1/4 * 2(f-1) * ∂f/∂v1)
00237         // ∂Z/∂v2 = 1/(2*z) * ((v1⋅v4)∂f/∂v2 + v4^2 * 1/4 * 2(f-1) * ∂f/∂v2)
00238         if (path_type == Z) {
00239             if (use_z_square) {
00240                 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];
00241                 dzdv2[i_elem] = v1v4 * dfdv2[i_elem] + v4v4 * 0.25 * 2.0 * (f-1) * dfdv2[i_elem];
00242             } else {
00243                 if (z > static_cast<scalar_type>(0)) {
00244                     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]);
00245                     dzdv2[i_elem] = (1.0 / (2.0 * z)) * (v1v4 * dfdv2[i_elem] + v4v4 * 0.25 * 2.0 * (f-1) * dfdv2[i_elem]);
00246                 } else {
00247                     // workaround at z = 0
00248                     dzdv1[i_elem] = 0;
00249                     dzdv2[i_elem] = 0;
00250                 }
00251             }
00252         }
00253     }
00254 }
00255 
00256 }
00257 
00258 #endif // GEOMETRICPATHCV_H