Merge pull request #381 from seoklab/jnooree/issue-361

feat(algo/crdgen): implement basic DG coordinate generation routine

cmake/NuriKitUtils.cmake

@@ -125,6 +125,26 @@ function(find_or_fetch_eigen)
   endif()
 endfunction()
 
+function(find_or_fetch_spectra)
+  set(BUILD_TESTING OFF)
+
+  find_package(Spectra 1.0 QUIET)
+
+  if(Spectra_FOUND)
+    message(STATUS "Found Spectra ${Spectra_VERSION}")
+  else()
+    include(FetchContent)
+    message(NOTICE "Could not find compatible Spectra. Fetching from github.")
+
+    FetchContent_Declare(
+      spectra
+      GIT_REPOSITORY https://github.com/yixuan/spectra.git
+      GIT_TAG v1.0.1)
+    nuri_make_available_deponly(spectra)
+    add_library(Spectra::Spectra ALIAS Spectra)
+  endif()
+endfunction()
+
 function(find_or_fetch_pybind11)
   set(BUILD_TESTING OFF)
 

include/nuri/algo/crdgen.h

@@ -0,0 +1,47 @@
+//
+// Project NuriKit - Copyright 2024 SNU Compbio Lab.
+// SPDX-License-Identifier: Apache-2.0
+//
+
+#ifndef NURI_ALGO_CRDGEN_H_
+#define NURI_ALGO_CRDGEN_H_
+
+/// @cond
+#include <utility>
+/// @endcond
+
+#include "nuri/eigen_config.h"
+#include "nuri/core/molecule.h"
+
+namespace nuri {
+/**
+ * @brief Generate 3D coordinates for the molecule.
+ *
+ * @param mol The molecule to generate coordinates.
+ * @param conf A matrix to store the generated coordinates.
+ * @param max_trial The maximum number of trials to generate trial distances.
+ * @return Whether the coordinates is successfully generated.
+ */
+extern bool generate_coords(const Molecule &mol, Matrix3Xd &conf,
+                            int max_trial = 10);
+
+/**
+ * @brief Generate 3D coordinates for the molecule.
+ *
+ * @param mol The molecule to generate coordinates.
+ * @return Whether the coordinates is successfully generated.
+ *
+ * @note The generated coordinates is stored as the last conformer of the
+ *       molecule, if and only if the coordinates is successfully generated. The
+ *       molecule is not modified otherwise.
+ */
+inline bool generate_coords(Molecule &mol, int max_trial = 10) {
+  Matrix3Xd conf;
+  bool success = generate_coords(mol, conf, max_trial);
+  if (success)
+    mol.confs().push_back(std::move(conf));
+  return success;
+}
+}  // namespace nuri
+
+#endif /* NURI_ALGO_CRDGEN_H_ */

include/nuri/algo/optim.h

@@ -233,7 +233,7 @@ namespace internal {
  */
 class LBfgsB {
 public:
-  LBfgsB(MutRef<ArrayXd> x, internal::LbfgsbBounds bounds, int m);
+  LBfgsB(MutRef<ArrayXd> x, internal::LbfgsbBounds bounds, int m = 10);
 
   /**
    * @brief Minimize a function using L-BFGS-B algorithm.
@@ -250,8 +250,8 @@ class LBfgsB {
    *         value, and final gradient.
    */
   template <class FuncGrad>
-  LbfgsbResult minimize(FuncGrad fg, double factr, int maxiter, int maxls,
-                        double pgtol);
+  LbfgsbResult minimize(FuncGrad fg, double factr = 1e+7, double pgtol = 1e-5,
+                        int maxiter = 15000, int maxls = 20);
 
   /* State modifiers, only for implementations */
 
@@ -384,8 +384,8 @@ class LBfgsB {
 
 template <class FuncGrad>
 LbfgsbResult LBfgsB::minimize(FuncGrad fg, const double factr,
-                              const int maxiter, const int maxls,
-                              const double pgtol) {
+                              const double pgtol, const int maxiter,
+                              const int maxls) {
   const double tol = factr * internal::kEpsMach;
 
   ArrayXd gx(n());
@@ -545,8 +545,8 @@ LbfgsbResult l_bfgs_b(FuncGrad &&fg, MutRef<ArrayXd> x, const ArrayXi &nbd,
     return { LbfgsbResultCode::kInvalidInput, 0, 0, {} };
 
   LBfgsB lbfgsb(x, { nbd, bounds }, m);
-  return lbfgsb.minimize(std::forward<FuncGrad>(fg), factr, maxiter, maxls,
-                         pgtol);
+  return lbfgsb.minimize(std::forward<FuncGrad>(fg), factr, pgtol, maxiter,
+                         maxls);
 }
 }  // namespace nuri
 

include/nuri/core/geometry.h

@@ -93,6 +93,35 @@ namespace constants {
 
   extern constexpr inline double kTwoPi =
       6.2831853071795864769252867665590057683943387987502116419498891846156328;
+
+  // NOLINTBEGIN(*-identifier-naming)
+  extern constexpr inline double kCos15 =
+      0.9659258262890682867497431997288973676339048390084045504023430763;
+  extern constexpr inline double kCos75 =
+      0.2588190451025207623488988376240483283490689013199305138140032073;
+  extern constexpr inline double kCos100 =
+      -0.173648177666930348851716626769314796000375677184069387236241378;
+  extern constexpr inline double kCos102 =
+      -0.207911690817759337101742284405125166216584760627723836407181973;
+  extern constexpr inline double kCos112 =
+      -0.374606593415912035414963774501195131000158922253676174103440371;
+  extern constexpr inline double kCos115 =
+      -0.422618261740699436186978489647730181563129301194864623444415159;
+  extern constexpr inline double kCos125 =
+      -0.573576436351046096108031912826157864620433371450986351081027118;
+  extern constexpr inline double kCos155 =
+      -0.906307787036649963242552656754316983267712625175864680871298408;
+  extern constexpr inline double kCos175 =
+      -0.996194698091745532295010402473888046183562672645850974525442277;
+  extern constexpr inline double kTan10_2 =
+      0.0874886635259240052220186694349614581194542763681082291452366622;
+  extern constexpr inline double kTan15_2 =
+      0.1316524975873958534715264574097171035928141022232375735535653257;
+  extern constexpr inline double kTan116_2 =
+      1.6003345290410503553267330811833575255040718469227591484115002297;
+  extern constexpr inline double kTan155_2 =
+      4.5107085036620571342899391172547519686713241944553043587162345185;
+  // NOLINTEND(*-identifier-naming)
 }  // namespace constants
 
 template <class DT, std::enable_if_t<std::is_floating_point_v<DT>, int> = 0>
@@ -519,6 +548,38 @@ extern std::pair<Affine3d, double>
 qcp(const Eigen::Ref<const Matrix3Xd> &query,
     const Eigen::Ref<const Matrix3Xd> &templ, AlignMode mode = AlignMode::kBoth,
     double evalprec = 1e-11, double evecprec = 1e-6, int maxiter = 50);
+
+/**
+ * @brief A routine for converting squared pairwise distances to cartesian
+ *        coordinates.
+ * @param pts Destination to which save the generated coordinates (3d).
+ * @param dsqs The squared distances between points. Will be modified in-place.
+ * @return Whether the embedding was successful.
+ *
+ * @note The squared distance matrix must be a N x N symmetric pairwise
+ *       squared-distance matrix, where N is the number of points.
+ *
+ * This implementation is based on the following reference: TF Havel, ID Kuntz,
+ * and GM Crippen. *Bull. Math. Biol.* **1983**, *45* (5), 665-720.
+ * DOI:[10.1007/BF02460044](https://doi.org/10.1007/BF02460044)
+ */
+extern bool embed_distances_3d(Eigen::Ref<Matrix3Xd> pts, MatrixXd &dsqs);
+
+/**
+ * @brief A routine for converting squared pairwise distances to cartesian
+ *        coordinates.
+ * @param pts Destination to which save the generated coordinates (4d).
+ * @param dsqs The squared distances between points. Will be modified in-place.
+ * @return Whether the embedding was successful.
+ *
+ * @note The squared distance matrix must be a N x N symmetric pairwise
+ *       squared-distance matrix, where N is the number of points.
+ *
+ * This implementation is based on the following reference: TF Havel, ID Kuntz,
+ * and GM Crippen. *Bull. Math. Biol.* **1983**, *45* (5), 665-720.
+ * DOI:[10.1007/BF02460044](https://doi.org/10.1007/BF02460044)
+ */
+extern bool embed_distances_4d(Eigen::Ref<Matrix4Xd> pts, MatrixXd &dsqs);
 }  // namespace nuri
 
 #endif /* NURI_CORE_GEOMETRY_H_ */

include/nuri/eigen_config.h

@@ -36,10 +36,14 @@ namespace nuri {
 using Eigen::Array;
 using Eigen::Array3d;
 using Eigen::Array3i;
+using Eigen::Array4d;
 using Eigen::Array4i;
 using Eigen::ArrayX;
 using ArrayXb = Eigen::ArrayX<bool>;
 using Eigen::Array2Xd;
+using Eigen::Array3Xd;
+using Eigen::Array4Xd;
+using Eigen::ArrayX2d;
 using Eigen::ArrayXd;
 using Eigen::ArrayXi;
 using Eigen::ArrayXX;

src/CMakeLists.txt

@@ -6,6 +6,7 @@
 include(NuriKitUtils)
 
 find_or_fetch_eigen()
+find_or_fetch_spectra()
 
 include(GNUInstallDirs)
 
@@ -25,9 +26,13 @@ add_library(nuri_lib SHARED "${NURI_SRCS}")
 set_target_properties(nuri_lib PROPERTIES OUTPUT_NAME nuri)
 target_link_libraries(nuri_lib
   PUBLIC
-  absl::strings absl::flat_hash_map absl::absl_check absl::absl_log
+  absl::strings
+  absl::flat_hash_map
+  absl::random_random
+  absl::absl_check
+  absl::absl_log
 )
-target_system_include_directories(nuri_lib Eigen3::Eigen)
+target_system_include_directories(nuri_lib Eigen3::Eigen Spectra::Spectra)
 handle_boost_dependency(nuri_lib)
 
 if(NURI_LIBRARY_FLAGS)