Release 6.0.0

.conan_test_package/conanfile.py

@@ -17,7 +17,7 @@ class TestPackageConan(ConanFile):
     exports_sources = "CMakeLists.txt", "test.cpp"
 
     def _configure(self):
-        cmake = CMake(self)
+        cmake = CMake(self, cmake_program='/usr/bin/cmake')
         cmake.configure()
         return cmake
 

.conan_test_package/test.cpp

@@ -3,7 +3,7 @@
 using namespace Scine::Utils;
 
 int main() {
-  if(ElementInfo::element(1) == ElementType::H) {
+  if (ElementInfo::element(1) == ElementType::H) {
     return 0;
   }
 

.gitignore

@@ -16,9 +16,6 @@ install/*
 # python files
 *__pycache__*
 
-# generated documentation
-resource/Documentation/*
-
 # generated source files
 externalQC_output_location.h
 orca_output_location.h

CHANGELOG.rst

@@ -8,6 +8,30 @@ 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 6.0.0
+-------------
+
+New Features and Feature Updates
+................................
+- Newton Trajectory: Added new extraction options and improved eta bonds in NT2
+- Added more solvents for the Turbomole input creator
+- Added the Pauling electronegativity scale. This is available through the ElementInfo.
+- Added solvate function to place any mixture of solvents around solute.
+
+Important Technical Changes
+...........................
+- PeriodicSystem canonicalizes the given PeriodicBoundaries to ensure read/write stability
+- Calculators:
+  - For the Turbomole calculator, allow the SCF damping value to be specified exactly (instead of predefined
+    settings "default", "low", "medium", and "high")
+  - Add "PBEH-3C" and "B97-3C" as supported method families for the ORCA calculator
+- Settings for placing solvent molecules are summarized in ``SolventPlacementSettings``, all solvate functions
+  take only this structure as input. This change is not backwards compatible.
+
+Development Changes
+...................
+- Code deduplication
+
 Release 5.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 5.0.0
+  VERSION 6.0.0
   DESCRIPTION "Utilities to be used (statically linked) in all other SCINE modules."
 )
 

CONTRIBUTING.rst

@@ -29,3 +29,5 @@ Further Contributors
 
 - Eric Hermes, Sandia National Laboratories
   (https://github.com/qcscine/utilities/pull/1)
+- Sebastian Ehlert (`@awvwgk <https://github.com/awvwgk>`_), University of Bonn
+  (https://github.com/qcscine/utilities/issues/4)

README.rst

@@ -1,5 +1,5 @@
-SCINE - Open Source Utilities
-=============================
+SCINE - Utilities
+=================
 
 Introduction
 ------------
@@ -19,7 +19,7 @@ Installation and Usage
 ----------------------
 
 As a user of one of the SCINE modules (such as Sparrow), you do not need
-to set up Open Source Utilities yourself. Instead, this is done as part of the
+to set up SCINE Utilities yourself. Instead, this is done as part of the
 installation process of the respective SCINE module. Therefore, the following
 instructions are only necessary for developers of SCINE, or those interfacing
 this library directly.

conanfile.py

@@ -9,7 +9,7 @@
 
 class ScineUtilsConan(ScineConan):
     name = "scine_utilities"
-    version = "5.0.0"
+    version = "6.0.0"
     url = "https://github.com/qcscine/utilities"
     description = """
 Functionality which is used in most SCINE modules. It is vital for the correct

src/Utils/CMakeLists.txt

@@ -166,7 +166,10 @@ if(SCINE_BUILD_PYTHON_BINDINGS)
     ${CMAKE_CURRENT_BINARY_DIR}/scine_utilities/pkginit.py
     ${CMAKE_CURRENT_BINARY_DIR}/scine_utilities/__init__.py
   )
-  file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/Python/README.rst DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+  file(
+    COPY ${PROJECT_SOURCE_DIR}/README.rst
+    DESTINATION ${CMAKE_CURRENT_BINARY_DIR}
+  )
   add_custom_command(TARGET UtilsOSModule POST_BUILD
     COMMAND ${CMAKE_COMMAND} -E copy $<TARGET_FILE:UtilsOSModule> ${CMAKE_CURRENT_BINARY_DIR}/scine_utilities
     COMMENT "Copying UtilsOS module into python package directory"

src/Utils/Python/ElementInfoPython.cpp

@@ -73,6 +73,8 @@ void init_element_info(pybind11::module& m) {
     )delim");
   element_info.def_static("vdw_radius", &ElementInfo::vdwRadius, pybind11::arg("element"),
                           "van der Waals radius in atomic units");
+  element_info.def_static("pauling_electronegativity", &ElementInfo::paulingElectronegativity, pybind11::arg("element"),
+                          "Pauling electronegativity");
   element_info.def_static("Z", &ElementInfo::Z, pybind11::arg("element"),
                           R"delim(
       Atomic number of an element

src/Utils/Python/MolecularTrajectoryPython.cpp

@@ -15,7 +15,18 @@ using namespace Scine::Utils;
 void init_molecular_trajectory(pybind11::module& m) {
   pybind11::class_<MolecularTrajectory> molecular_trajectory(m, "MolecularTrajectory");
 
-  molecular_trajectory.def(pybind11::init<>());
+  // constructors
+  molecular_trajectory.def(pybind11::init<>(), "Initialize completely empty trajectory.");
+  molecular_trajectory.def(pybind11::init<const ElementTypeCollection&>(), pybind11::arg("elements"),
+                           "Initialize empty trajectory with given elements.");
+  molecular_trajectory.def(pybind11::init<double>(), pybind11::arg("minimum_rmsd_for_addition"),
+                           "Initialize empty trajectory with a minimum root mean square deviation for a "
+                           "PositionCollection to the previous one to be added.");
+  molecular_trajectory.def(
+      pybind11::init<const ElementTypeCollection&, double>(), pybind11::arg("elements"),
+      pybind11::arg("minimum_rmsd_for_addition"),
+      "Initialize empty trajectory with given elements and a minimum root mean square deviation for a "
+      "PositionCollection to the previous one to be added.");
 
   molecular_trajectory.def_property("elements", &MolecularTrajectory::getElementTypes,
                                     &MolecularTrajectory::setElementTypes, "Element types of the atoms");

src/Utils/Python/PeriodicBoundariesPython.cpp

@@ -58,9 +58,17 @@ void init_periodic_boundaries(pybind11::module& m) {
   periodic_boundaries.def_property_readonly("beta", &PeriodicBoundaries::getBeta, "Unit cell angle beta between a and c");
   periodic_boundaries.def_property_readonly("gamma", &PeriodicBoundaries::getGamma, "Unit cell angle gamma between a and b");
 
+  periodic_boundaries.def_property_readonly("lengths", &PeriodicBoundaries::getLengths, "Lengths of the three unit vectors");
+  periodic_boundaries.def_property_readonly("angles", &PeriodicBoundaries::getAngles,
+                                            "Angles between the three unit vectors in degrees");
+
   periodic_boundaries.def_property("matrix", &PeriodicBoundaries::getCellMatrix, &PeriodicBoundaries::setCellMatrix,
                                    "The underlying matrix governing the periodic boundaries.");
 
+  periodic_boundaries.def_property("periodicity", &PeriodicBoundaries::getPeriodicity,
+                                   pybind11::overload_cast<const std::vector<bool>>(&PeriodicBoundaries::setPeriodicity),
+                                   "The periodicity of the cell.");
+
   periodic_boundaries.def(
       "is_ortho_rhombic", &PeriodicBoundaries::isOrthoRhombic, pybind11::arg("eps") = 1e-2,
       "Returns whether the cell is orthorhombic. The optional parameter gives the tolerance around 90 degrees.");

src/Utils/Python/SoluteSolventComplexPython.cpp

@@ -15,22 +15,43 @@ using namespace Scine::Utils::SoluteSolventComplex;
 void init_solute_solvent_complex(pybind11::module& m) {
   auto solvation_submodule = m.def_submodule("solvation");
 
+  pybind11::class_<SolventPlacementSettings>(solvation_submodule, "placement_settings")
+      .def(pybind11::init<>())
+      .def("__repr__", [](const SolventPlacementSettings&) { return "solvent_placement_settings"; })
+      .def_readwrite("resolution", &SolventPlacementSettings::resolution)
+      .def_readwrite("solvent_offset", &SolventPlacementSettings::solventOffset)
+      .def_readwrite("max_distance", &SolventPlacementSettings::maxDistance)
+      .def_readwrite("step_size", &SolventPlacementSettings::stepSize)
+      .def_readwrite("num_rotamers", &SolventPlacementSettings::numRotamers)
+      .def_readwrite("strategic_solvation", &SolventPlacementSettings::strategicSolv)
+      .def_readwrite("coverage_threshold", &SolventPlacementSettings::coverageThreshold);
+
   solvation_submodule.def(
       "solvate",
-      pybind11::overload_cast<const AtomCollection&, int, const AtomCollection&, int, int, int, double, double, double, int, bool, double>(
-          &solvate),
+      pybind11::overload_cast<const AtomCollection&, int, const AtomCollection&, int, int, SolventPlacementSettings>(&solvate),
       pybind11::arg("solute_complex"), pybind11::arg("solute_size"), pybind11::arg("solvent"),
-      pybind11::arg("number_solvents"), pybind11::arg("seed"), pybind11::arg("resolution") = 32,
-      pybind11::arg("solvent_offset") = 0.0, pybind11::arg("max_distance") = 10.0, pybind11::arg("step_size") = 0.25,
-      pybind11::arg("number_rotamers") = 3, pybind11::arg("strategic_solvation") = false,
-      pybind11::arg("coverage_threshold") = 1.0, "Add systematically a number of solvents to solute.");
+      pybind11::arg("number_solvents"), pybind11::arg("seed"),
+      pybind11::arg("placement_settings") = SolventPlacementSettings(),
+      "Add systematically a number of one type of "
+      "solvent to solute.");
+
+  solvation_submodule.def("solvate_mix", &solvateMix, pybind11::arg("solute_complex"), pybind11::arg("solute_size"),
+                          pybind11::arg("solvents"), pybind11::arg("solvent_ratios"), pybind11::arg("number_solvents"),
+                          pybind11::arg("seed"), pybind11::arg("placement_settings") = SolventPlacementSettings(),
+                          "Add systematically a"
+                          "number of different solvents to solute.");
 
   solvation_submodule.def("solvate_shells", &solvateShells, pybind11::arg("solute_complex"),
                           pybind11::arg("solute_size"), pybind11::arg("solvent"), pybind11::arg("number_shells"),
-                          pybind11::arg("seed"), pybind11::arg("resolution") = 32, pybind11::arg("solvent_offset") = 0.0,
-                          pybind11::arg("max_distance") = 10.0, pybind11::arg("step_size") = 0.25,
-                          pybind11::arg("number_rotamers") = 3, pybind11::arg("strategic_solvation") = false,
-                          pybind11::arg("coverage_threshold") = 1.0, "Add number of solvent shells to solute.");
+                          pybind11::arg("seed"), pybind11::arg("placement_settings") = SolventPlacementSettings(),
+                          "Add number of one type of solvent in shells to solute.");
+
+  solvation_submodule.def("solvate_shells_mix", &solvateShellsMix, pybind11::arg("solute_complex"),
+                          pybind11::arg("solute_size"), pybind11::arg("solvents"), pybind11::arg("solvent_ratios"),
+                          pybind11::arg("number_shells"), pybind11::arg("seed"),
+                          pybind11::arg("placement_settings") = SolventPlacementSettings(),
+                          "Add different"
+                          "solvents in shells to solute.");
 
   solvation_submodule.def("give_solvent_shell_vector", &giveSolventShellVector, pybind11::arg("complex"),
                           pybind11::arg("solute_size"), pybind11::arg("solvent_size_vector"), pybind11::arg("resolution"),
@@ -48,6 +69,10 @@ void init_solute_solvent_complex(pybind11::module& m) {
       "merge_solvent_shell_vector", &mergeSolventShellVector, pybind11::arg("shell_vector"),
       "Merge a vector of a vector of atom collections (solvent shell vector) to one atom collection.");
 
+  solvation_submodule.def("merge_solvent_shell_indices_vector", &mergeSolventShellIndices,
+                          pybind11::arg("shell_indices_vector"),
+                          "Merge a vector of a vector of indices (solvent shell indices vector) to a one flat list.");
+
   solvation_submodule.def("check_distances", &checkDistances, pybind11::arg("molecule_1"), pybind11::arg("molecule_2"),
                           "Check if two atom collections overlap with their VdW spheres.");
 

src/Utils/Python/Tests/test_MolecularTrajectory.py

@@ -27,8 +27,7 @@ def test_access():
     elements = [scine.ElementType.H, scine.ElementType.F]
     pos_a = numpy.array([[0.0, 1.0, 2.0], [0.5, 1.5, 2.5]])
     pos_b = numpy.array([[1.0, 2.0, 3.0], [1.5, 2.5, 3.5]])
-    traj = scine.MolecularTrajectory()
-    traj.elements = elements
+    traj = scine.MolecularTrajectory(elements)
     assert traj.elements[0] == scine.ElementType.H
     assert traj.elements[1] == scine.ElementType.F
     traj.push_back(pos_a)
@@ -56,3 +55,35 @@ def test_access():
         _ = traj[-10]
     assert "out of range" in str(excinfo.value)
 
+
+def test_min_addition():
+    elements = [scine.ElementType.H, scine.ElementType.F]
+    pos_a = numpy.array([[0.0, 1.0, 2.0], [0.5, 1.5, 2.5]])
+    pos_b = numpy.array([[1.0, 2.0, 3.0], [1.5, 2.5, 3.5]])
+    pos_c = numpy.array([[10.0, 2.0, 3.0], [1.5, 2.5, 3.5]])
+    traj = scine.MolecularTrajectory(elements, 1.5)
+    assert traj.elements[0] == scine.ElementType.H
+    assert traj.elements[1] == scine.ElementType.F
+    traj.push_back(pos_a)
+    traj.push_back(pos_b)
+    traj.push_back(pos_b)
+    traj.push_back(pos_b)
+    traj.push_back(pos_a)
+    traj.push_back(pos_c)
+    assert traj.size() == 4
+    stricter_traj = scine.MolecularTrajectory(5.0)
+    stricter_traj.elements = elements
+    assert len(stricter_traj.elements) == 2
+    stricter_traj.push_back(pos_a)
+    stricter_traj.push_back(pos_b)
+    stricter_traj.push_back(pos_b)
+    stricter_traj.push_back(pos_b)
+    stricter_traj.push_back(pos_a)
+    assert stricter_traj.size() == 1
+    stricter_traj.push_back(pos_c)
+    assert stricter_traj.size() == 2
+    stricter_traj.push_back(pos_a)
+    assert stricter_traj.size() == 3
+    stricter_traj.push_back(pos_a)
+    assert stricter_traj.size() == 3
+

src/Utils/Python/setup.py

@@ -57,7 +57,6 @@ def collect_data(pkg_name: str) -> Dict[str, List[str]]:
     author_email="[email protected]",
     description="Utilities used in all SCINE modules",
     long_description=long_description,
-    long_description_content_type="text/markdown",
     url="https://www.scine.ethz.ch",
     packages=["scine_utilities"],
     package_data=collect_data("scine_utilities"),

src/Utils/Tests/CalculatorBasics/CalculationRoutinesTest.cpp

@@ -77,6 +77,11 @@ TEST_F(ACalculationRoutinesTest, DispersionSplitting) {
   ASSERT_THAT(outCorrect.first, std::string{"pbe"});
   ASSERT_THAT(outCorrect.second, std::string{"d3bj"});
 
+  std::string correctException = "dlpno-ccsd(t)";
+  auto outCorrectException = CalculationRoutines::splitIntoMethodAndDispersion(correctException);
+  ASSERT_THAT(outCorrectException.first, std::string{"dlpno-ccsd(t)"});
+  ASSERT_THAT(outCorrectException.second, std::string{""});
+
   std::string noDisp = "pbe";
   auto outNoDisp = CalculationRoutines::splitIntoMethodAndDispersion(noDisp);
   ASSERT_THAT(outNoDisp.first, std::string{"pbe"});

src/Utils/Tests/DataStructures/PeriodicBoundariesTest.cpp

@@ -187,13 +187,14 @@ TEST_F(APeriodicBoundariesTest, ClassIsInitializedCorrectly) {
   ASSERT_TRUE(pbc1 == pbc2);
   ASSERT_DOUBLE_EQ(pbc1.getA().norm(), pbc2.getA().norm());
   ASSERT_DOUBLE_EQ(pbc1.getAlpha(), pbc2.getAlpha());
+  pbc1.setPeriodicity("xy");
   pbc4 = pbc1;
   ASSERT_TRUE(pbc1 == pbc4);
   ASSERT_DOUBLE_EQ(pbc1.getA().norm(), pbc4.getA().norm());
   ASSERT_DOUBLE_EQ(pbc1.getAlpha(), pbc4.getAlpha());
-  PeriodicBoundaries pbc8 = PeriodicBoundaries(Eigen::Matrix3d::Identity());
+  PeriodicBoundaries pbc8 = PeriodicBoundaries(Eigen::Matrix3d::Identity(), pbc1.getPeriodicityString());
   pbc1 = Eigen::Matrix3d::Identity();
-  ASSERT_THAT(pbc1, pbc8);
+  ASSERT_TRUE(pbc1 == pbc8);
 }
 
 TEST_F(APeriodicBoundariesTest, OrthoRhombicWorks) {
@@ -203,6 +204,40 @@ TEST_F(APeriodicBoundariesTest, OrthoRhombicWorks) {
   ASSERT_FALSE(pbc2.isOrthoRhombic());
 }
 
+TEST_F(APeriodicBoundariesTest, CanonicalizationWorks) {
+  Eigen::Matrix3d randomMatrix = Eigen::Matrix3d::Random(3, 3);
+  for (int i = 0; i < 3; ++i) {
+    randomMatrix(i, i) = std::fabs(randomMatrix(i, i));
+  }
+  auto pbc = PeriodicBoundaries(randomMatrix);
+  auto canonic = PeriodicBoundaries(pbc.getPeriodicBoundariesString());
+  auto rot = pbc.getCanonicalizationRotationMatrix();
+  Eigen::Matrix3d manualMatrix = pbc.getCellMatrix() * rot;
+  ASSERT_TRUE(std::fabs(canonic.getCellMatrix()(0, 2) < 1e-12));
+  ASSERT_TRUE(std::fabs(canonic.getCellMatrix()(1, 2) < 1e-12));
+  ASSERT_TRUE(std::fabs(manualMatrix(0, 2) < 1e-12));
+  ASSERT_TRUE(std::fabs(manualMatrix(1, 2) < 1e-12));
+  ASSERT_TRUE(canonic.getCellMatrix().isApprox(manualMatrix, 1e-6));
+  pbc.canonicalize();
+  ASSERT_TRUE(canonic.getCellMatrix().isApprox(pbc.getCellMatrix(), 1e-6));
+  ASSERT_TRUE(std::fabs(pbc.getCellMatrix()(0, 2) < 1e-12));
+  ASSERT_TRUE(std::fabs(pbc.getCellMatrix()(1, 2) < 1e-12));
+  ASSERT_TRUE(pbc.getCanonicalizationRotationMatrix().isApprox(Eigen::Matrix3d::Identity()));
+}
+
+TEST_F(APeriodicBoundariesTest, LengthsAndAngles) {
+  Eigen::Vector3d lengths = lengths2 * Constants::bohr_per_angstrom;
+  auto pbc = PeriodicBoundaries(lengths, angles2);
+  ASSERT_THAT(pbc.getLengths().size(), 3);
+  ASSERT_THAT(pbc.getAngles().size(), 3);
+  ASSERT_DOUBLE_EQ(pbc.getLengths()[0], lengths[0]);
+  ASSERT_DOUBLE_EQ(pbc.getLengths()[1], lengths[1]);
+  ASSERT_DOUBLE_EQ(pbc.getLengths()[2], lengths[2]);
+  ASSERT_DOUBLE_EQ(pbc.getAngles()[0], angles2[0]);
+  ASSERT_DOUBLE_EQ(pbc.getAngles()[1], angles2[1]);
+  ASSERT_DOUBLE_EQ(pbc.getAngles()[2], angles2[2]);
+}
+
 TEST_F(APeriodicBoundariesTest, InverseIsConsistent) {
   PeriodicBoundaries pbc = PeriodicBoundaries(lengths1, angles1, false);
   ASSERT_TRUE(pbc.getInverseCellMatrix().isApprox(pbc.getCellMatrix().inverse()));