ENH: Kill GSL

Co-authored-by: Amrita Goswami <[email protected]>

Author: HaoZeke

environment.yml

@@ -15,7 +15,6 @@ dependencies:
   - liblapack
   - libblas
   - boost-cpp
-  - gsl
   - cmake
   - lua-luafilesystem
   # Pinned

nix/yodaStruct.nix

@@ -21,7 +21,6 @@
 , blas
 , lib
 , boost
-, gsl
 , cmake }:
   clangStdenv.mkDerivation {
   name = "yodaStruct";
@@ -38,7 +37,6 @@
   eigen
   catch2
   boost
-  gsl
   liblapack
   blas
   lua

src/CMakeLists.txt

@@ -17,7 +17,6 @@ set(LIBRARY_OUTPUT_PATH ${CMAKE_BINARY_DIR})
 include(FindBLAS)
 include(FindLAPACK)
 include(FindBoost)
-include(FindGSL)
 add_compile_options(
   # "-Wall" "-Wpedantic" "-Wextra" "-fexceptions"
   "$<$<CONFIG:Debug>:-O0;-g3;-ggdb>"
@@ -34,12 +33,10 @@ include_directories(
   ${LUA_INCLUDE_DIR}
   ${Boost_INCLUDE_DIRS}
   ${Eigen3_INCLUDE_DIRS}
-  ${GSL_INCLUDE_DIRS}
   )
 add_definitions(-DBOOST_ERROR_CODE_HEADER_ONLY)
 find_package(
   "Lua" REQUIRED
-  "GSL" REQUIRED
   "Boost" REQUIRED COMPONENTS system filesystem
   "Eigen3 3.3" REQUIRED
   "yaml-cpp"
@@ -52,8 +49,6 @@ target_link_libraries(
   ${Boost_LIBRARIES}
   ${Eigen3_LIBRARIES}
   ${LUA_LIBRARIES}
-  ${GSL_LIBRARIES}
-  stdc++fs
   fmt
   yaml-cpp
   yodaLib

src/bond.cpp

@@ -292,12 +292,12 @@ bond::populateHbonds(std::string filename,
         } // end of applying PBCs to the O-H and O--O vectors
         //
         // Get the angle between the O--O and O-H vectors
-        double gslAngle = gen::gslVecAngle(ooVec, ohVec);
-        double gslAngleDeg = gen::radDeg(gslAngle);
+        double eigenAngle = gen::eigenVecAngle(ooVec, ohVec);
+        double eigenAngleDeg = gen::radDeg(eigenAngle);
 
         //
         // A hydrogen bond is formed if the angle is less than 30 degrees
-        if (gslAngleDeg > angleCutoff) {
+        if (eigenAngleDeg > angleCutoff) {
           continue;
         } // not a hydrogen bond
 

src/generic.cpp

@@ -140,27 +140,18 @@ int gen::unwrappedCoordShift(
 
 /**
  * @details Function for obtaining the angle between two input vectors
- * (std::vector). Internally, the vectors are converted to GSL vectors. The dot
- * product between the input vectors is used to calculate the angle between
+ * (std::vector). Internally, the vectors are converted to Eigen vectors. The
+ * dot product between the input vectors is used to calculate the angle between
  * them.
  * @param[in] OO The O--O vector (but can be any vector, in general).
  * @param[in] OH The O-H vector (but can be any vector, in general).
  * @return The output angle between the input vectors, in radians
  */
-double gen::gslVecAngle(std::vector<double> OO, std::vector<double> OH) {
-  gsl_vector *gOO = gsl_vector_alloc(3);
-  gsl_vector *gOH = gsl_vector_alloc(3);
-  double norm_gOO, norm_gOH, xDummy, angle;
-  for (int i = 0; i < 3; i++) {
-    gsl_vector_set(gOO, i, OO[i]);
-    gsl_vector_set(gOH, i, OH[i]);
-  }
-  norm_gOO = gsl_blas_dnrm2(gOO);
-  norm_gOH = gsl_blas_dnrm2(gOH);
-  gsl_blas_ddot(gOO, gOH, &xDummy);
-  angle = acos(xDummy / (norm_gOO * norm_gOH));
-  gsl_vector_free(gOO);
-  gsl_vector_free(gOH);
+double gen::eigenVecAngle(std::vector<double> OO, std::vector<double> OH) {
+  Eigen::Vector3d eigOO = Eigen::Map<Eigen::Vector3d>(OO.data(), OO.size());
+  Eigen::Vector3d eigOH = Eigen::Map<Eigen::Vector3d>(OH.data(), OH.size());
+  double angle;
+  angle = acos(eigOO.dot(eigOH) / (eigOO.norm() * eigOH.norm()));
   return angle;
 }
 

src/include/internal/generic.hpp

@@ -18,7 +18,9 @@
 
 // Boost
 #include <boost/math/constants/constants.hpp>
-#include <gsl/gsl_blas.h>
+// Eigen
+#include <eigen3/Eigen/Core>
+#include <eigen3/Eigen/Dense>
 
 /** @file generic.hpp
  *   @brief File for containing generic or common functions.
@@ -54,8 +56,8 @@ const double pi = boost::math::constants::pi<double>();
  */
 inline double radDeg(double angle) { return (angle * 180) / gen::pi; }
 
-//! GSL for getting the angle (in radians) between the O--O and O-H vectors
-double gslVecAngle(std::vector<double> OO, std::vector<double> OH);
+//! Eigen function for getting the angle (in radians) between the O--O and O-H vectors
+double eigenVecAngle(std::vector<double> OO, std::vector<double> OH);
 
 //! Get the average, after excluding the outliers, using quartiles
 double getAverageWithoutOutliers(std::vector<double> inpVec);

src/include/internal/mol_sys.hpp

@@ -20,6 +20,7 @@
 #include <string>
 #include <sys/stat.h>
 #include <vector>
+#include<unordered_map>
 
 // For debugging, instantiate the unordered map [consider removal for
 // production]

tests/CMakeLists.txt

@@ -31,15 +31,13 @@ add_executable(yodaStruct_test
 # include(FindBLAS)
 # include(FindLAPACK)
 # include(FindBoost)
-include(FindGSL)
 
 # Find packages
 find_package(
   "Lua" REQUIRED
   "Boost" REQUIRED
   "Eigen3 3.3" REQUIRED
   "Catch2" REQUIRED
-  "GSL" REQUIRED
   "yaml-cpp"
   "fmt"
   )
@@ -49,7 +47,6 @@ target_link_libraries(yodaStruct_test
   ${Boost_LIBRARIES}
   ${Eigen3_LIBRARIES}
   ${LUA_LIBRARIES}
-  ${GSL_LIBRARIES}
   fmt
   yaml-cpp
   yodaLib
@@ -60,7 +57,6 @@ target_link_libraries(yodaStruct_test
 include_directories(
   ${PROJECT_SOURCE_DIR}/src/include/internal
   ${PROJECT_SOURCE_DIR}/src/include/external
-  ${GSL_INCLUDE_DIRS}
   ${Eigen3_INCLUDE_DIRS}
   )
 # Run unit tests