Constraints included in projector

Works with Eigen, not with Armadillo (singular matrix)
RMeli · Jun 7, 2018 · 8f7f1cb · 8f7f1cb
1 parent 2b31dfc
commit 8f7f1cb
Show file tree

Hide file tree

Showing 3 changed files with 75 additions and 4 deletions.
diff --git a/include/libirc/connectivity.h b/include/libirc/connectivity.h
@@ -971,7 +971,7 @@ std::vector<Dihedral> dihedrals(const Matrix& distance_m,
   // Check if dihedrals are found
   if (n_atoms >= 4 && dih.empty()) {
     std::cerr << "ERROR: Out of plane bending not implemented yet."
-              << std::endl;
+              << std::endl << std::flush;
   }
 
   // Return list of dihedral angles

diff --git a/include/libirc/irc.h b/include/libirc/irc.h
@@ -10,6 +10,8 @@
 #include <algorithm>
 #include <utility>
 
+#include <boost/optional.hpp>
+
 namespace irc {
 
 template<typename Vector3, typename Vector, typename Matrix>
@@ -97,6 +99,10 @@ class IRC {
 
   /// Wilson B matrix
   Matrix B;
+
+  // TODO: Move to std::optional with C++17
+  /// Constraint matrix
+  boost::optional<Matrix> C;
 
   /// Projector
   Matrix P;
@@ -133,6 +139,45 @@ size_t add_without_duplicates(std::vector<T>& v1, const std::vector<T>& v2){
   return n;
 }
 
+// TODO: Switch to std::optional with C++17
+template<typename Matrix>
+boost::optional<Matrix> constraints(const std::vector<connectivity::Bond>& B,
+    const std::vector<connectivity::Angle>& A,
+    const std::vector<connectivity::Dihedral>& D)
+{
+  std::size_t n{B.size() + A.size() + D.size()};
+
+  Matrix C{linalg::zeros<Matrix>(n, n)};
+
+  bool constrained{false};
+
+  std::size_t offset{0};
+  for(std::size_t i{0}; i < B.size(); i++){
+    if(B[i].constraint == connectivity::Constraint::constrained){
+      C(i + offset,i + offset) = 1;
+      constrained=true;
+    }
+  }
+
+  offset = B.size();
+  for(std::size_t i{0}; i < A.size(); i++){
+    if(A[i].constraint == connectivity::Constraint::constrained){
+      C(i + offset,i + offset) = 1;
+      constrained=true;
+    }
+  }
+
+  offset = B.size() + A.size();
+  for(std::size_t i{0}; i < D.size(); i++){
+    if(D[i].constraint == connectivity::Constraint::constrained){
+      C(i + offset,i + offset) = 1;
+      constrained=true;
+    }
+  }
+
+  return constrained ? boost::optional<Matrix>(C) : boost::none;
+}
+
 template<typename Vector3, typename Vector, typename Matrix>
 IRC<Vector3, Vector, Matrix>::IRC(
     const molecule::Molecule<Vector3>& molecule,
@@ -185,9 +230,17 @@ IRC<Vector3, Vector, Matrix>::IRC(
       bonds,
       angles,
       dihedrals);
+
+  // Compute (optional) constraint matrix
+  C = constraints<Matrix>(bonds, angles, dihedrals);
 
   // Compute projector P
-  P = wilson::projector(B);
+  if(C){
+    P = wilson::projector(B, *C);
+  }
+  else{
+    P = wilson::projector(B);
+  }
 }
 
 /// Initial estimate of the inverse Hessian in internal redundant coordinates
@@ -345,7 +398,12 @@ Vector IRC<Vector3, Vector, Matrix>::irc_to_cartesian(const Vector& q_irc_old,
       itc_result.x_c, bonds, angles, dihedrals);
 
   // Update projector P
-  P = wilson::projector(B);
+  if(C){
+    P = wilson::projector(B, *C);
+  }
+  else{
+    P = wilson::projector(B);
+  }
 
   // Return new cartesian coordinates
   return itc_result.x_c;

diff --git a/include/libirc/wilson.h b/include/libirc/wilson.h
@@ -351,10 +351,23 @@ Matrix wilson_matrix_numerical(
 /// \return Projector
 template<typename Matrix>
 Matrix projector(const Matrix& B) {
-  // TODO: Pass iB instead of computing it
   return B * linalg::pseudo_inverse(B);
 }
 
+/// Compute projector from the \param B matrix with constraint \param C
+/// \tparam Matrix
+/// \param B Wilson's B matrix
+/// \param C Constraint matrix
+/// \return Projector
+template<typename Matrix>
+Matrix projector(const Matrix& B, const Matrix& C) {
+  // Standard projector
+  Matrix P{B * linalg::pseudo_inverse(B)};
+
+  // Projector with constraints
+  return P - P * C * linalg::inv<Matrix>(C * P * C) * C * P;
+}
+
 } // namespace wilson
 
 } // namespace irc