Release 5.0.0

CHANGELOG.rst

@@ -8,6 +8,34 @@ It is intended that only the first two groups (``New Features and Feature Update
 ``Important Technical Changes``) are important for the average user, while
 the last one is mainly aimed at developers and users that link deeply into the code.
 
+Release 5.0.0
+-------------
+
+New Features and Feature Updates
+................................
+- Add second version of Newton Trajectory algorithm (NT2) for elementary step induction
+- Add Gaussian process regression
+- Improve BFGS: guarantees positive-definite approximate Hessian
+- B-Splines: interpolation and compression of trajectories including energy
+- Enable optimization of a geometry together with its unit cell
+  - CP2K:
+    - Add xTB support
+    - Add stress tensor
+- Enable point charge embedding for Turbomole
+- Add linear sum assignment algorithm
+- Add functionality to evaluate the spin contamination of a single determinant
+- Add Python bindings for WavefunctionOutputGenerator and casting utilities
+  from Calculator and CalculatorWithReference.
+
+Important Technical Changes
+...........................
+- Calculators:
+  - Harmonize dispersion correction input
+
+Development Changes
+...................
+- Code deduplication
+
 Release 4.0.0
 -------------
 

CMakeLists.txt

@@ -4,7 +4,7 @@ cmake_minimum_required(VERSION 3.9)
 # tree must then provide a properly namespaced target with the same name as
 # your project.
 project(UtilsOS
-  VERSION 4.0.0
+  VERSION 5.0.0
   DESCRIPTION "Utilities to be used (statically linked) in all other SCINE modules."
 )
 

README.rst

@@ -29,7 +29,7 @@ Dependencies
 
 Required software, minimum required versions in brackets, for this SCINE project are:
 
-- A C++ compiler supporting the C++14 standard
+- A C++ compiler supporting the C++17 standard
 - CMake (3.9)
 - Boost (1.65.0)
 - Eigen3 (3.3.2)
@@ -42,7 +42,7 @@ commands::
 
     git submodule init
     git submodule update
-    mkdir build inst
+    mkdir build
     cd build
     cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=../inst ..
     make -j 4
@@ -72,5 +72,6 @@ SCINE Utilities makes use of the following third-party libraries:
 - `Eigen <http://eigen.tuxfamily.org>`_
 - `Google Test <https://github.com/google/googletest>`_
 - `IRC <https://github.com/rmeli/irc>`_
+- `LBFGSB <https://github.com/yixuan/LBFGSpp>`_
 - `pybind11 <https://github.com/pybind/pybind11>`_
 - `yaml-cpp <https://github.com/jbeder/yaml-cpp>`_

conanfile.py

@@ -4,12 +4,12 @@
 See LICENSE.txt for details.
 """
 
+import sys
 from dev.conan.base import ScineConan
 
-
 class ScineUtilsConan(ScineConan):
     name = "scine_utilities"
-    version = "4.0.0"
+    version = "5.0.0"
     url = "https://github.com/qcscine/utilities"
     description = """
 Functionality which is used in most SCINE modules. It is vital for the correct
@@ -21,16 +21,19 @@ class ScineUtilsConan(ScineConan):
         "python": [True, False],
         "tests": [True, False],
         "coverage": [True, False],
-        "microarch": ["detect", "none"]
+        "microarch": ["detect", "none"],
+        "python_version": "ANY"
     }
+    python_version_string = str(sys.version_info.major) + \
+        "." + str(sys.version_info.minor) 
     default_options = {"shared": True, "python": False,
                        "tests": False, "coverage": False,
-                       "microarch": "none"}
+                       "microarch": "none", "python_version": python_version_string}
     exports = "dev/conan/*.py"
     exports_sources = ["dev/cmake/*", "src/*", "CMakeLists.txt", "README.rst",
                        "LICENSE.txt", "dev/conan/hook.cmake",
                        "dev/conan/glue/*"]
-    requires = ["eigen/[~=3.3.7]", "boost/[>1.65.0]", "yaml-cpp/0.6.3"]
+    requires = ["eigen/[~=3.3.7]", "boost/[>1.65.0]", "yaml-cpp/0.6.3", "lbfgspp/0.1.0"]
     cmake_name = "UtilsOS"
     cmake_definitions = {
         "CMAKE_UNITY_BUILD": "ON",

src/Utils/CMakeLists.txt

@@ -16,6 +16,8 @@ include(ImportYamlCpp)
 import_yamlcpp()
 include(ImportIRC)
 import_irc()
+include(ImportLBFGSB)
+import_lbfgspp()
 include(AddEigen)
 include(ImportGTest)
 import_gtest()
@@ -50,6 +52,7 @@ target_link_libraries(UtilsOS
     $<$<BOOL:${OpenMP_CXX_FOUND}>:OpenMP::OpenMP_CXX>
     yaml-cpp
     irc
+    lbfgspp
 )
 
 target_include_directories(UtilsOS
@@ -74,7 +77,7 @@ set_target_properties(UtilsOSModule PROPERTIES
 # Headers
 install(
   DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/Utils/
-  DESTINATION include/Scine/Utils/
+  DESTINATION include/Scine/Utils
   FILES_MATCHING PATTERN "*.h"
 )
 

src/Utils/Python/AtomCollectionPython.cpp

@@ -52,35 +52,58 @@ void init_atom_collection(pybind11::module& m) {
       pybind11::keep_alive<0, 1>() // Keep object alive while iterator exists
   );
   atom_collection.def("__len__", &AtomCollection::size);
-  atom_collection.def("__getitem__", &AtomCollection::operator[]);
-  atom_collection.def("__getitem__", [&](AtomCollection& coll, pybind11::slice slice) -> AtomCollection {
-    std::size_t start = 0, stop = 0, step = 0, slicelength = 0;
-    if (!slice.compute(coll.size(), &start, &stop, &step, &slicelength)) {
-      throw pybind11::error_already_set();
-    }
-    AtomCollection sliced(slicelength);
-    for (std::size_t i = 0; i < slicelength; ++i, start += step) {
-      sliced.setElement(i, coll.getElement(start));
-      sliced.setPosition(i, coll.getPosition(start));
-    }
-    return sliced;
-  });
-  atom_collection.def("__delitem__", [&](AtomCollection& coll, int i) {
-    const int N = coll.size();
-    if (i >= N) {
-      throw std::out_of_range("Cannot delete out of range index");
-    }
-    AtomCollection result(N - 1);
-    for (int j = 0; j < i; ++j) {
-      result.setElement(j, coll.getElement(j));
-      result.setPosition(j, coll.getPosition(j));
-    }
-    for (int j = i; j < N - 1; ++j) {
-      result.setElement(j, coll.getElement(j + 1));
-      result.setPosition(j, coll.getPosition(j + 1));
-    }
-    coll = result;
-  });
+  atom_collection.def(
+      "__getitem__",
+      [&](AtomCollection& coll, int i) -> Atom {
+        if (i >= coll.size()) {
+          throw std::out_of_range("Given index is out of range");
+        }
+        if (i < 0) {
+          if (-1 * i > coll.size()) {
+            throw std::out_of_range("Given index is out of range");
+          }
+          return coll[coll.size() + i];
+        }
+        return coll[i];
+      },
+      "Access an Atom of the AtomCollection.");
+  atom_collection.def(
+      "__getitem__",
+      [&](AtomCollection& coll, pybind11::slice slice) -> AtomCollection {
+        std::size_t start = 0, stop = 0, step = 0, slicelength = 0;
+        if (!slice.compute(coll.size(), &start, &stop, &step, &slicelength)) {
+          throw pybind11::error_already_set();
+        }
+        AtomCollection sliced(slicelength);
+        for (std::size_t i = 0; i < slicelength; ++i, start += step) {
+          sliced.setElement(i, coll.getElement(start));
+          sliced.setPosition(i, coll.getPosition(start));
+        }
+        return sliced;
+      },
+      "Access a sub-AtomCollection of the AtomCollection based on slicing.");
+  atom_collection.def(
+      "__delitem__",
+      [&](AtomCollection& coll, int i) {
+        const int N = coll.size();
+        if (i >= N || -1 * i > N) {
+          throw std::out_of_range("Cannot delete out of range index");
+        }
+        if (i < 0) {
+          i = N + i;
+        }
+        AtomCollection result(N - 1);
+        for (int j = 0; j < i; ++j) {
+          result.setElement(j, coll.getElement(j));
+          result.setPosition(j, coll.getPosition(j));
+        }
+        for (int j = i; j < N - 1; ++j) {
+          result.setElement(j, coll.getElement(j + 1));
+          result.setPosition(j, coll.getPosition(j + 1));
+        }
+        coll = result;
+      },
+      "Allow the python delete function based on index.");
 
   // Copy support
   atom_collection.def("__copy__", [](const AtomCollection& c) -> AtomCollection { return AtomCollection(c); });

src/Utils/Python/BSplinesPython.cpp

@@ -0,0 +1,105 @@
+/**
+ * @file
+ * @copyright This code is licensed under the 3-clause BSD license.\n
+ *            Copyright ETH Zurich, Laboratory of Physical Chemistry, Reiher Group.\n
+ *            See LICENSE.txt for details.
+ */
+#include <Utils/Geometry/AtomCollection.h>
+#include <Utils/Math/BSplines/ReactionProfileInterpolation.h>
+#include <Utils/Pybind.h>
+#include <pybind11/eigen.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+using namespace Scine;
+
+void init_bspline_functionalities(pybind11::module& m) {
+  auto bsplines = m.def_submodule("bsplines");
+  bsplines.doc() = R"delim(
+    A collection of functions to generate and work with B-Splines of molecular
+    trajectories.
+  )delim";
+
+  pybind11::class_<Utils::BSplines::TrajectorySpline> bspline{bsplines, "TrajectorySpline", R"delim(
+      A class representing a B-Spline fit.
+    )delim"};
+
+  bspline.def(pybind11::init<Utils::ElementTypeCollection, Eigen::VectorXd, Eigen::MatrixXd>());
+
+  bspline.def(pybind11::init<Utils::ElementTypeCollection, Eigen::VectorXd, Eigen::MatrixXd, double>());
+
+  bspline.def("evaluate", &Utils::BSplines::TrajectorySpline::evaluate, pybind11::arg("position"), pybind11::arg("degree") = 0,
+              R"delim(
+      Fit a spline to all currently stored data.
+
+      :param position: The position at which to evaluate the spline. has to be
+                       in the interval [0,1]
+      :param degree: The degree of the polynomial fit.
+      :return: A Tuple of the fitted data: (energy: double , structure: AtomCollection).
+    )delim");
+
+  bspline.def_readonly("elements", &Utils::BSplines::TrajectorySpline::elements,
+                       R"delim(
+      The elements of the atoms that are represented.
+    )delim");
+
+  bspline.def_readonly("data", &Utils::BSplines::TrajectorySpline::data,
+                       R"delim(
+      The spline data to be fitted to.
+    )delim");
+
+  bspline.def_readonly("ts_position", &Utils::BSplines::TrajectorySpline::tsPosition,
+                       R"delim(
+      The TS position in the spline [0.0, 1.0], -1.0 if none is present.
+    )delim");
+
+  bspline.def_readonly("knots", &Utils::BSplines::TrajectorySpline::knots,
+                       R"delim(
+      The knot positions of the spline.
+    )delim");
+
+  pybind11::class_<Utils::BSplines::ReactionProfileInterpolation> rpi{bsplines, "ReactionProfileInterpolation",
+                                                                      R"delim(
+      Factory for B-Splines of reaction paths, including interpolation of an
+      energy associated with the molecular structures.
+    )delim"};
+
+  rpi.def(pybind11::init<>());
+
+  rpi.def("append_structure", &Utils::BSplines::ReactionProfileInterpolation::appendStructure, pybind11::arg("atoms"),
+          pybind11::arg("energy"), pybind11::arg("is_the_transition_state") = false,
+          R"delim(
+      Adds a datapoint to the internal storage.
+      Given atoms have to match those previously given, if there were any.
+
+      :param atoms: An atom collection
+      :param energy: The energy of the given atomic configuration.
+    )delim");
+
+  rpi.def("clear", &Utils::BSplines::ReactionProfileInterpolation::clear,
+          R"delim(
+      Clear all previously stored data.
+    )delim");
+
+  rpi.def("spline", &Utils::BSplines::ReactionProfileInterpolation::spline,
+          pybind11::arg("n_interpolation_points") = 11, pybind11::arg("degree") = 3,
+          R"delim(
+      Fit a spline to all currently stored data.
+
+      :param n_interpolation_points: Number of points to use for interpolation.
+      :param degree: The maximum degree of the fit polynominals.
+      :return: A BSpline object fitted to all stored data.
+    )delim");
+
+  rpi.def("current_ts_position", &Utils::BSplines::ReactionProfileInterpolation::getCurrentTSPosition,
+          R"delim(
+      :return: The current position of the transiton state in the intverval [0.0, 1.0]
+    )delim");
+
+  rpi.def_readwrite("use_quaternion_fit", &Utils::BSplines::ReactionProfileInterpolation::useQuaternionFit,
+                    R"delim(
+      If true, will determine distance along path by fitting each structure
+      to the previous one. This fit is done in mass-weighted coordinates using
+      quaternions.
+    )delim");
+}

src/Utils/Python/Bonds/BondOrderCollectionPython.cpp

@@ -75,6 +75,10 @@ void init_bond_order_collection(pybind11::module& m) {
     )delim");
   bond_order_collection.def("get_order", &getBOOrder, pybind11::arg("i"), pybind11::arg("j"), "Get a bond order");
   bond_order_collection.def("empty", &BondOrderCollection::empty, "Checks whether there are no entries in the BO matrix");
+  bond_order_collection.def(
+      "set_to_absolute_values", &BondOrderCollection::setToAbsoluteValues,
+      "Transfers all values to absolute values. Bonds across periodic boundaries can be signalled by negative "
+      "bond orders. This method essentially removes this information.");
 
   // Operators
   bond_order_collection.def(pybind11::self == pybind11::self);

src/Utils/Python/ConstantsPython.cpp

@@ -14,6 +14,9 @@ void init_constants(pybind11::module& m) {
   m.attr("ELEMENTARY_CHARGE") = Constants::elementaryCharge;
   m.attr("AVOGADRO_NUMBER") = Constants::avogadroNumber;
   m.attr("PI") = Constants::pi;
+  m.attr("BOLTZMANN_CONSTANT") = Constants::boltzmannConstant;
+  m.attr("PLANCK_CONSTANT") = Constants::planckConstant;
+  m.attr("MOLAR_GAS_CONSTANT") = Constants::molarGasConstant;
   m.attr("ATOMIC_MASS_UNIT") = Constants::atomicMassUnit;
   m.attr("ELECTRON_REST_MASS") = Constants::electronRestMass;
   m.attr("INVERSE_FINE_STRUCTURE_CONSTANT") = Constants::inverseFineStructureConstant;