diff --git a/.conan_test_package/conanfile.py b/.conan_test_package/conanfile.py index bed3673..9bcdca4 100644 --- a/.conan_test_package/conanfile.py +++ b/.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 diff --git a/.conan_test_package/test.cpp b/.conan_test_package/test.cpp index 2f38876..2cfde64 100644 --- a/.conan_test_package/test.cpp +++ b/.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; } diff --git a/.gitignore b/.gitignore index 2b67c52..24d0fac 100644 --- a/.gitignore +++ b/.gitignore @@ -16,9 +16,6 @@ install/* # python files *__pycache__* -# generated documentation -resource/Documentation/* - # generated source files externalQC_output_location.h orca_output_location.h diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 58ea310..0548519 100644 --- a/CHANGELOG.rst +++ b/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 ------------- diff --git a/CMakeLists.txt b/CMakeLists.txt index 3989161..dd5fac0 100644 --- a/CMakeLists.txt +++ b/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." ) diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst index eb615a4..5f5993f 100644 --- a/CONTRIBUTING.rst +++ b/CONTRIBUTING.rst @@ -29,3 +29,5 @@ Further Contributors - Eric Hermes, Sandia National Laboratories (https://github.com/qcscine/utilities/pull/1) +- Sebastian Ehlert (`@awvwgk `_), University of Bonn + (https://github.com/qcscine/utilities/issues/4) diff --git a/README.rst b/README.rst index 9779be9..3e9ad60 100644 --- a/README.rst +++ b/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. diff --git a/conanfile.py b/conanfile.py index 84c396e..d1bcaec 100644 --- a/conanfile.py +++ b/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 diff --git a/dev b/dev index abb0acc..0ffba74 160000 --- a/dev +++ b/dev @@ -1 +1 @@ -Subproject commit abb0acc9e1bb05da2a8dd4c5804d415b7f596e1a +Subproject commit 0ffba74b61c40c6df9ce4289e6585c4fe766b3a8 diff --git a/src/Utils/CMakeLists.txt b/src/Utils/CMakeLists.txt index 5056d5a..85bc41c 100644 --- a/src/Utils/CMakeLists.txt +++ b/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 $ ${CMAKE_CURRENT_BINARY_DIR}/scine_utilities COMMENT "Copying UtilsOS module into python package directory" diff --git a/src/Utils/Python/ElementInfoPython.cpp b/src/Utils/Python/ElementInfoPython.cpp index b6c08f2..ebd9888 100644 --- a/src/Utils/Python/ElementInfoPython.cpp +++ b/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 diff --git a/src/Utils/Python/MolecularTrajectoryPython.cpp b/src/Utils/Python/MolecularTrajectoryPython.cpp index f030b6a..e154bf8 100644 --- a/src/Utils/Python/MolecularTrajectoryPython.cpp +++ b/src/Utils/Python/MolecularTrajectoryPython.cpp @@ -15,7 +15,18 @@ using namespace Scine::Utils; void init_molecular_trajectory(pybind11::module& m) { pybind11::class_ molecular_trajectory(m, "MolecularTrajectory"); - molecular_trajectory.def(pybind11::init<>()); + // constructors + molecular_trajectory.def(pybind11::init<>(), "Initialize completely empty trajectory."); + molecular_trajectory.def(pybind11::init(), pybind11::arg("elements"), + "Initialize empty trajectory with given elements."); + molecular_trajectory.def(pybind11::init(), 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(), 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"); diff --git a/src/Utils/Python/PeriodicBoundariesPython.cpp b/src/Utils/Python/PeriodicBoundariesPython.cpp index c405f44..b42d5fb 100644 --- a/src/Utils/Python/PeriodicBoundariesPython.cpp +++ b/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>(&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."); diff --git a/src/Utils/Python/README.rst b/src/Utils/Python/README.rst deleted file mode 100644 index 4ef299f..0000000 --- a/src/Utils/Python/README.rst +++ /dev/null @@ -1,24 +0,0 @@ -SCINE - Utilities Python Bindings -================================= - -Introduction ------------- - -This repository contains functionality which is used in most SCINE modules. -It is vital for the correct functioning of all SCINE modules, but it does not -directly provide any functionality for end users. - - -Installation and Usage ----------------------- - -We currently offer only a ``manylinux2014`` build for Linux platforms. This may -be expanded in the future. If you are on a different platform, you can build -the Python bindings from source (https://github.com/qcscine/utilities/). - - -Support and Contact -------------------- - -In case you should encounter problems or bugs, please write a short message -to scine@phys.chem.ethz.ch. diff --git a/src/Utils/Python/SoluteSolventComplexPython.cpp b/src/Utils/Python/SoluteSolventComplexPython.cpp index d57fc33..f5f77c8 100644 --- a/src/Utils/Python/SoluteSolventComplexPython.cpp +++ b/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_(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( - &solvate), + pybind11::overload_cast(&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."); diff --git a/src/Utils/Python/Tests/test_MolecularTrajectory.py b/src/Utils/Python/Tests/test_MolecularTrajectory.py index bdb988a..809dde9 100644 --- a/src/Utils/Python/Tests/test_MolecularTrajectory.py +++ b/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 + diff --git a/src/Utils/Python/setup.py b/src/Utils/Python/setup.py index 7d16b12..8909537 100644 --- a/src/Utils/Python/setup.py +++ b/src/Utils/Python/setup.py @@ -57,7 +57,6 @@ def collect_data(pkg_name: str) -> Dict[str, List[str]]: author_email="scine@phys.chem.ethz.ch", 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"), diff --git a/src/Utils/Tests/CalculatorBasics/CalculationRoutinesTest.cpp b/src/Utils/Tests/CalculatorBasics/CalculationRoutinesTest.cpp index 4cc6a58..c023446 100644 --- a/src/Utils/Tests/CalculatorBasics/CalculationRoutinesTest.cpp +++ b/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"}); diff --git a/src/Utils/Tests/DataStructures/PeriodicBoundariesTest.cpp b/src/Utils/Tests/DataStructures/PeriodicBoundariesTest.cpp index d468e87..5329c20 100644 --- a/src/Utils/Tests/DataStructures/PeriodicBoundariesTest.cpp +++ b/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())); diff --git a/src/Utils/Tests/Geometry/PeriodicSystemTest.cpp b/src/Utils/Tests/Geometry/PeriodicSystemTest.cpp index b2c2f11..661fc6a 100644 --- a/src/Utils/Tests/Geometry/PeriodicSystemTest.cpp +++ b/src/Utils/Tests/Geometry/PeriodicSystemTest.cpp @@ -7,6 +7,7 @@ #include "Utils/Bonds/BondDetector.h" #include "Utils/Geometry/GeometryUtilities.h" #include "Utils/Geometry/PeriodicSystem.h" +#include "Utils/Math/QuaternionFit.h" #include #include @@ -52,7 +53,8 @@ class PeriodicSystemTest : public Test { TEST_F(PeriodicSystemTest, CanBeConstructed) { auto ps0 = PeriodicSystem(pbc); ASSERT_THAT(ps0.atoms.size(), 0); - ASSERT_THAT(ps0.pbc, pbc); + pbc.canonicalize(); + ASSERT_TRUE(ps0.pbc == pbc); auto ps1 = PeriodicSystem(pbc, 1); ASSERT_THAT(ps1.atoms.size(), 1); std::unordered_set unimportant = {1}; @@ -62,6 +64,26 @@ TEST_F(PeriodicSystemTest, CanBeConstructed) { ASSERT_THAT(ps2.atoms.getElements(), randomElements); } +TEST_F(PeriodicSystemTest, CanonicalizationWorks) { + auto atoms = AtomCollection(randomElements, randomPositions); + auto ps = PeriodicSystem(pbc, atoms); + auto copyPbc = PeriodicBoundaries(pbc.getPeriodicBoundariesString()); + if (std::fabs(pbc.getCellMatrix()(0, 2)) > 1e-6 || std::fabs(pbc.getCellMatrix()(1, 2)) > 1e-6) { + ASSERT_FALSE(pbc == copyPbc); + } + else { + ASSERT_TRUE(pbc == copyPbc); + } + auto rotation = pbc.getCanonicalizationRotationMatrix(); + PositionCollection oldPositions = atoms.getPositions(); + atoms.setPositions(atoms.getPositions() * rotation); + auto fit = QuaternionFit(oldPositions, atoms.getPositions()); + ASSERT_LT(fit.getRMSD(), 1e-5); + auto copyPs = PeriodicSystem(copyPbc, atoms); + ASSERT_TRUE(ps.pbc.getCellMatrix().isApprox(copyPs.pbc.getCellMatrix(), 1e-5)); + ASSERT_TRUE(ps.atoms.getPositions().isApprox(copyPs.atoms.getPositions(), 1e-5)); +} + TEST_F(PeriodicSystemTest, GetsCorrectImageAtoms) { /* * | H H H | @@ -79,7 +101,7 @@ TEST_F(PeriodicSystemTest, GetsCorrectImageAtoms) { auto ps = PeriodicSystem(pbc, h3); BondOrderCollection bo = BondDetector::detectBonds(h3, pbc, true); - ASSERT_THAT(bo, ps.constructBondOrders()); + ASSERT_TRUE(bo == ps.constructBondOrders()); ps.constructImageAtoms(bo); ps.constructImageAtoms(); auto atoms = ps.getAtomCollectionWithImages(); @@ -508,9 +530,9 @@ TEST_F(PeriodicSystemTest, MultiplicationWorksSurface) { srand(42); int factor = std::rand() % 10 + 1; auto otherPs = ps * factor; - ASSERT_THAT(otherPs, otherPs); + ASSERT_TRUE(otherPs == otherPs); ASSERT_THAT(otherPs.atoms.size(), structure.size() * std::pow(factor, 3)); - ASSERT_THAT(otherPs.pbc.getCellMatrix(), originalCell * factor); + ASSERT_TRUE(otherPs.pbc.getCellMatrix().isApprox(originalCell * factor)); ASSERT_THAT(otherPs.solidStateAtomIndices.size(), solidStateIndices.size() * std::pow(factor, 3)); std::vector randoms; for (int i = 0; i < 3; ++i) { @@ -545,9 +567,9 @@ TEST_F(PeriodicSystemTest, MultiplicationWorksMolecular) { int factor = std::rand() % 10 + 1; auto otherPs = ps * factor; ps *= factor; - ASSERT_THAT(ps, otherPs); + ASSERT_TRUE(ps == otherPs); ASSERT_THAT(ps.atoms.size(), structure.size() * std::pow(factor, 3)); - ASSERT_THAT(ps.pbc.getCellMatrix(), Eigen::Matrix3d::Identity() * factor * basisLength); + ASSERT_TRUE(ps.pbc.getCellMatrix().isApprox(Eigen::Matrix3d::Identity() * factor * basisLength)); ASSERT_TRUE(ps.solidStateAtomIndices.empty()); auto rm = Eigen::Matrix3d::Random().cwiseAbs(); Eigen::Matrix3d randomMatrix = rm; diff --git a/src/Utils/Tests/GeometryOptimizer/NtOptimizer2Test.cpp b/src/Utils/Tests/GeometryOptimizer/NtOptimizer2Test.cpp index 0d841ae..f90ff42 100644 --- a/src/Utils/Tests/GeometryOptimizer/NtOptimizer2Test.cpp +++ b/src/Utils/Tests/GeometryOptimizer/NtOptimizer2Test.cpp @@ -12,11 +12,13 @@ #include "Utils/CalculatorBasics/PropertyList.h" #include "Utils/CalculatorBasics/Results.h" #include "Utils/CalculatorBasics/TestCalculator.h" +#include "Utils/IO/MolecularTrajectoryIO.h" #include "Utils/Math/QuaternionFit.h" #include "Utils/Settings.h" #include "Utils/Technical/CloneInterface.h" /* External */ #include +#include using namespace testing; namespace Scine { @@ -120,10 +122,7 @@ TEST(NtOptimizer2Tests, DefaultEmptyLists) { ASSERT_THROW(nt.optimize(atoms, logger), std::logic_error); } -TEST(NtOptimizer2Tests, SN2) { - Core::Log logger = Core::Log::silent(); - TestCalculator testCalc; - NtOptimizer2 nt(testCalc); +inline void runSN2Test(NtOptimizer2& nt, Core::Log& logger) { nt.coordinateSystem = CoordinateSystem::Cartesian; nt.associationList = {0, 1}; nt.dissociationList = {1, 5}; @@ -170,6 +169,29 @@ TEST(NtOptimizer2Tests, SN2) { EXPECT_NEAR(referenceBonds.row(3).value(), (positions.row(1) - positions.row(5)).norm(), 1e-3); // C-Br } +TEST(NtOptimizer2Tests, SN2) { + Core::Log logger = Core::Log::silent(); + TestCalculator testCalc; + NtOptimizer2 nt(testCalc); + runSN2Test(nt, logger); +} + +TEST(NtOptimizer2Tests, SN2FirstMaximum) { + Core::Log logger = Core::Log::silent(); + TestCalculator testCalc; + NtOptimizer2 nt(testCalc); + nt.extractionCriterion = NtOptimizer2::ntExtractFirst; + runSN2Test(nt, logger); +} + +TEST(NtOptimizer2Tests, SN2HighestMaximum) { + Core::Log logger = Core::Log::silent(); + TestCalculator testCalc; + NtOptimizer2 nt(testCalc); + nt.extractionCriterion = NtOptimizer2::ntExtractHighest; + runSN2Test(nt, logger); +} + TEST(NtOptimizer2Tests, SingleStepFailure) { Core::Log logger = Core::Log::silent(); NtOpt2MockCalculator mockCalculator; @@ -183,10 +205,76 @@ TEST(NtOptimizer2Tests, SingleStepFailure) { // Enforce a scan that stops after the initial evaluation. positions(0, 1) = 0.1; AtomCollection atoms(elements, positions); - // Expect failure du to missing TS (trajectory has only one point) + // Expect failure due to missing TS (trajectory has only one point) ASSERT_THROW(nt.optimize(atoms, logger), std::runtime_error); } +TEST(NtOptimizer2Tests, ConnectedNucleiWorks) { + Core::Log logger = Core::Log::silent(); + NtOpt2MockCalculator mockCalculator; + mockCalculator.setLog(logger); + auto pathToResources = boost::dll::program_location().parent_path(); + pathToResources /= "Resources"; + const auto trajectoryFile = (pathToResources / "dissociation.trj.xyz").string(); + auto trajectory = Utils::MolecularTrajectoryIO::read(MolecularTrajectoryIO::format::xyz, trajectoryFile); + // testing connected state + auto atomsFirst = Utils::AtomCollection(trajectory.getElementTypes(), trajectory.front()); + mockCalculator.setStructure(atomsFirst); + auto result = mockCalculator.calculate(); + ASSERT_TRUE(result.has()); + const auto boFirst = result.get(); + std::vector indices = {0, 2, 4}; + auto molecules = NtUtils::connectedNuclei(indices, boFirst); + ASSERT_THAT(molecules.size(), 1); + ASSERT_THAT(molecules[0], indices); + // testing unconnected state + auto atomsLast = Utils::AtomCollection(trajectory.getElementTypes(), trajectory.back()); + mockCalculator.setStructure(atomsLast); + result = mockCalculator.calculate(); + ASSERT_TRUE(result.has()); + const auto boLast = result.get(); + indices = {1, 9}; + molecules = NtUtils::connectedNuclei(indices, boLast); + ASSERT_THAT(molecules.size(), 2); + ASSERT_TRUE((molecules[0] == std::vector{indices[0]} && molecules[1] == std::vector{indices[1]}) || + (molecules[1] == std::vector{indices[0]} && molecules[0] == std::vector{indices[1]})); +} + +TEST(NtOptimizer2Tests, Center2CenterWorks) { + PositionCollection pos = PositionCollection::Zero(4, 3); + // clang-format off + pos << -1.0, 0.0, 0.0, + 1.0, 2.0, 10.0, + 1.0, -2.0, 10.0, + 4.2, 3.1, -0.9; + // clang-format on + Displacement expected; + expected << -2.0, 0.0, -10.0; + std::vector lhs = {0}; + std::vector rhs = {1, 2}; + auto c2c = NtUtils::centerToCenterVector(pos, lhs, rhs); + ASSERT_TRUE(expected.isApprox(c2c)); +} + +TEST(NtOptimizer2Tests, SmallestCovRadiusWorks) { + auto elements = ElementTypeCollection{ElementType::Cl, ElementType::C, ElementType::H, + ElementType::H, ElementType::H, ElementType::Br}; + PositionCollection startPositions = PositionCollection::Zero(6, 3); + // clang-format off + startPositions << 3.7376961460e+00, 2.1020866350e-04, 4.5337439168e-02, + -3.8767703481e+00, -2.4803422157e-05, -1.2049608882e-01, + -2.3620148614e+00, 1.3238308540e+00, 1.0376490681e-01, + -2.3809041075e+00, -8.2773666259e-01, 9.6331578315e-01, + -2.3309449521e+00, -4.9652606314e-01, -1.3293307598e+00, + -7.4798903722e+00, 2.6536371103e-04, -1.9897114399e-01; + // clang-format on + AtomCollection atoms(elements, startPositions); + std::vector i = {0, 1, 5}; + std::vector j = {1, 2, 3}; + ASSERT_DOUBLE_EQ(NtUtils::smallestCovalentRadius(atoms, i), ElementInfo::covalentRadius(ElementType::C)); + ASSERT_DOUBLE_EQ(NtUtils::smallestCovalentRadius(atoms, j), ElementInfo::covalentRadius(ElementType::H)); +} + } /* namespace Tests */ } /* namespace Utils */ } /* namespace Scine */ diff --git a/src/Utils/Tests/Module/TurbomoleTest.cpp b/src/Utils/Tests/Module/TurbomoleTest.cpp index 34d5e53..7146af0 100644 --- a/src/Utils/Tests/Module/TurbomoleTest.cpp +++ b/src/Utils/Tests/Module/TurbomoleTest.cpp @@ -18,7 +18,6 @@ #include #include #include -#include #include using namespace testing; @@ -584,7 +583,7 @@ TEST_F(ATurbomoleTest, ImprovedScfConvergenceSettingsAreAppliedCorrectly) { const char* envVariablePtr = std::getenv("TURBODIR"); if (envVariablePtr) { calculator.settings().modifyDouble(Utils::ExternalQC::SettingsNames::scfOrbitalShift, 0.5); - calculator.settings().modifyString(Utils::SettingsNames::scfDamping, "high"); + calculator.settings().modifyBool(Utils::SettingsNames::scfDamping, true); calculator.settings().modifyInt(Utils::SettingsNames::maxScfIterations, 250); calculator.settings().modifyDouble(Utils::SettingsNames::selfConsistenceCriterion, 1e-4); @@ -618,7 +617,7 @@ TEST_F(ATurbomoleTest, ImprovedScfConvergenceSettingsAreAppliedCorrectly) { input.close(); std::regex scfIterRegex(R"((scfiterlimit)+ +([0-9]+))"); - std::regex scfDamp(R"(scfdamp start=8.500 step=0.10 min=0.50)"); + std::regex scfDamp(R"(scfdamp start=0.500 step=0.05 min=0.10)"); std::regex scfOrbitalShift(R"((scforbitalshift closedshell=)++([-\.0-9]+))"); std::regex convThreshold(R"((scfconv)+ +([0-9]+))"); diff --git a/src/Utils/Tests/MolecularTrajectoryTest.cpp b/src/Utils/Tests/MolecularTrajectoryTest.cpp new file mode 100644 index 0000000..55f898d --- /dev/null +++ b/src/Utils/Tests/MolecularTrajectoryTest.cpp @@ -0,0 +1,74 @@ +/** + * @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/MolecularTrajectory.h" +#include "Utils/Typenames.h" +#include "gmock/gmock.h" + +using namespace testing; +namespace Scine { +namespace Utils { +namespace Tests { + +/** + * @class MolecularTrajectoryTest MolecularTrajectoryTest.cpp + * @brief Comprises separate tests for the MolecularTrajectory class, + * which is mostly tested in the MolecularDynamicsTest + * @test + */ +class MolecularTrajectoryTest : public Test { + public: + ElementTypeCollection elements; + PositionCollection p1, p2, p3; + + protected: + void SetUp() override { + elements = {ElementType::H, ElementType::F}; + p1 = Eigen::MatrixX3d::Zero(2, 3); + p2 = Eigen::MatrixX3d::Zero(2, 3); + p3 = Eigen::MatrixX3d::Zero(2, 3); + // clang-format off + p1 << 0.0, 1.0, 2.0, + 0.5, 1.5, 2.5; + p2 << 1.0, 2.0, 3.0, + 1.5, 2.5, 3.5; + p3 << 10.0, 2.0, 3.0, + 1.5, 2.5, 3.5; + // clang-format on + } +}; + +TEST_F(MolecularTrajectoryTest, RmsdLimitWorks) { + auto traj = MolecularTrajectory(elements, 1.5); + ASSERT_TRUE(traj.getElementTypes().size() == 2); + ASSERT_TRUE(traj.getElementTypes()[0] == ElementType::H); + ASSERT_TRUE(traj.getElementTypes()[1] == ElementType::F); + traj.push_back(p1); + traj.push_back(p2); + traj.push_back(p2); + traj.push_back(p2); + traj.push_back(p1); + traj.push_back(p3); + ASSERT_TRUE(traj.size() == 4); + auto stricter_traj = MolecularTrajectory(5.0); + stricter_traj.setElementTypes(elements); + stricter_traj.push_back(p1); + stricter_traj.push_back(p2); + stricter_traj.push_back(p2); + stricter_traj.push_back(p2); + stricter_traj.push_back(p1); + ASSERT_TRUE(stricter_traj.size() == 1); + stricter_traj.push_back(p3); + ASSERT_TRUE(stricter_traj.size() == 2); + stricter_traj.push_back(p1); + ASSERT_TRUE(stricter_traj.size() == 3); + stricter_traj.push_back(p1); + ASSERT_TRUE(stricter_traj.size() == 3); +} + +} /* namespace Tests */ +} /* namespace Utils */ +} /* namespace Scine */ diff --git a/src/Utils/Tests/Resources/dissociation.trj.xyz b/src/Utils/Tests/Resources/dissociation.trj.xyz index 265de5f..6e64332 100644 --- a/src/Utils/Tests/Resources/dissociation.trj.xyz +++ b/src/Utils/Tests/Resources/dissociation.trj.xyz @@ -614,4 +614,3 @@ H -0.1121260277 0.2586275150 1.7132975753 H -1.6267430639 -0.6671056597 1.5832306617 O -2.0817216425 0.1235493398 -0.8338138002 H -2.6899663147 0.8589072273 -0.9844087902 -12 \ No newline at end of file diff --git a/src/Utils/Tests/Solvation/SoluteSolventComplexTest.cpp b/src/Utils/Tests/Solvation/SoluteSolventComplexTest.cpp index afb43de..ee9b0ce 100644 --- a/src/Utils/Tests/Solvation/SoluteSolventComplexTest.cpp +++ b/src/Utils/Tests/Solvation/SoluteSolventComplexTest.cpp @@ -4,14 +4,10 @@ * Copyright ETH Zurich, Laboratory of Physical Chemistry, Reiher Group.\n * See LICENSE.txt for details. */ -#include "Utils/IO/ChemicalFileFormats/ChemicalFileHandler.h" -#include "Utils/IO/ChemicalFileFormats/XyzStreamHandler.h" #include "Utils/Solvation/RandomIndexGenerator.h" #include "Utils/Solvation/SoluteSolventComplex.h" #include #include -#include -#include #include using namespace testing; @@ -32,6 +28,17 @@ static void assertPositionCollectionsEqual(const PositionCollection& p1, const P } } +TEST(SoluteSolventComplexTest, CheckDefaultSettings) { + SoluteSolventComplex::SolventPlacementSettings settings; + ASSERT_THAT(settings.resolution, 18); + ASSERT_THAT(settings.solventOffset, 0.0); + ASSERT_THAT(settings.maxDistance, 5.0); + ASSERT_THAT(settings.stepSize, 1.0); + ASSERT_THAT(settings.numRotamers, 3); + ASSERT_THAT(settings.strategicSolv, true); + ASSERT_THAT(settings.coverageThreshold, 0.85); +} + TEST(SoluteSolventComplexTest, DoesCheckDistanceWork) { Atom testAtom(ElementType::Xe); @@ -183,6 +190,20 @@ TEST(SoluteSolventComplexTest, DoesStrategyDeciderWork) { ASSERT_THAT(SoluteSolventComplex::solvationStrategy(0), 1); } +TEST(SoluteSolventComplexTest, EmptyRandomIndexGeneratorAndSetter) { + int size = 100; + int seed = 40; + SoluteSolventComplex::RandomIndexGenerator empty; + empty.setSeed(seed); + empty.setSize(size); + int i = 1; + while (i < 101) { + int randInt = empty.next(); + ASSERT_THAT(true, randInt >= 0 && randInt < 100); + i += 1; + } +} + // Test if Random Index Generator produces same number with same seed and if the number lies between 0 and size - 1 TEST(SoluteSolventComplexTest, DoesRandomIndexGeneratorWork) { int size = 100; @@ -201,6 +222,28 @@ TEST(SoluteSolventComplexTest, DoesRandomIndexGeneratorWork) { } } +TEST(SoluteSolventComplexTest, DoesDifferentListSizeThrowError) { + try { + SoluteSolventComplex::getSolventIndices(10, {4, 5, 3}, 2); + } + catch (const std::runtime_error& re) { + ASSERT_STREQ(re.what(), "The same number of ratios and solvents has to be given."); + } +} + +TEST(SoluteSolventComplexTest, GetCorrectIndexList) { + std::vector testOneRound = SoluteSolventComplex::getSolventIndices(6, {3, 2, 1}, 3); + std::vector test_match = {0, 0, 0, 1, 1, 2}; + for (unsigned long int i = 0; i < test_match.size(); i++) { + ASSERT_THAT(testOneRound.at(i), test_match.at(i)); + } + std::vector testTwoRounds = SoluteSolventComplex::getSolventIndices(8, {3, 2, 1}, 3); + std::vector test_match2 = {0, 0, 0, 1, 1, 2, 0, 0}; + for (unsigned long int i = 0; i < test_match2.size(); i++) { + ASSERT_THAT(testTwoRounds.at(i), test_match2.at(i)); + } +} + TEST(SoluteSolventComplexTest, DoesMergeAtomCollectionVectorWork) { AtomCollection test1; AtomCollection test2; @@ -264,6 +307,17 @@ TEST(SoluteSolventComplexTest, DoesMergeSolventVectorWork) { } } +TEST(SoluteSolventComplexTest, DoesMergeSolventIndicesWork) { + std::vector> test = {{0, 0, 1}, {2, 0, 1, 1}, {0, 0, 1, 1, 2}}; + std::vector control = {0, 0, 1, 2, 0, 1, 1, 0, 0, 1, 1, 2}; + + std::vector flatTest = SoluteSolventComplex::mergeSolventShellIndices(test); + + for (int i = 0; i < int(control.size()); i++) { + ASSERT_THAT(flatTest[i], control[i]); + } +} + TEST(SoluteSolventComplexTest, DoesSolvateAddsCorrectNumberOfSolvents) { ElementTypeCollection waterEC(3); waterEC.at(0) = ElementType::O; @@ -280,11 +334,18 @@ TEST(SoluteSolventComplexTest, DoesSolvateAddsCorrectNumberOfSolvents) { AtomCollection solute; solute.push_back(Atom(ElementType::Na, Position(0, 0, 0))); - int sampleRes = 18; + SoluteSolventComplex::SolventPlacementSettings testSettings; + testSettings.resolution = 18; + testSettings.solventOffset = 0.0; + testSettings.maxDistance = 7.0; + testSettings.stepSize = 1.0; + testSettings.numRotamers = 3; + testSettings.strategicSolv = false; + testSettings.coverageThreshold = 1.0; + int numSolvents = 24; - auto solventComplex = - SoluteSolventComplex::solvate(solute, solute.size(), water, numSolvents, 5, sampleRes, 0, 7, 1, 3, false); + auto solventComplex = SoluteSolventComplex::solvate(solute, solute.size(), water, numSolvents, 5, testSettings); ASSERT_THAT(numSolvents, SoluteSolventComplex::mergeSolventShellVector(solventComplex).size() / water.size()); @@ -293,18 +354,19 @@ TEST(SoluteSolventComplexTest, DoesSolvateAddsCorrectNumberOfSolvents) { auto iterComplex = solute; while (i < 18) { - auto newComplex = SoluteSolventComplex::solvate(iterComplex, solute.size(), water, 1, i + 5, sampleRes, 0, 7, 1, 3); + auto newComplex = SoluteSolventComplex::solvate(iterComplex, solute.size(), water, 1, i + 5, testSettings); iterComplex += SoluteSolventComplex::mergeSolventShellVector(newComplex); i++; } - auto iterComplexSurface = MolecularSurface::getVisibleMolecularSurface(iterComplex, 0, solute.size(), sampleRes); + auto iterComplexSurface = + MolecularSurface::getVisibleMolecularSurface(iterComplex, 0, solute.size(), testSettings.resolution); ASSERT_THAT((iterComplex.size() - solute.size()) / water.size(), 18); ASSERT_THAT(iterComplexSurface.size(), 0); } -// Run this test only in release builds +//// Run this test only in release builds #ifdef NDEBUG TEST(SoluteSolventComplexTest, DoesSolvateFinishesSolventShell) { ElementTypeCollection waterEC(3); @@ -322,19 +384,26 @@ TEST(SoluteSolventComplexTest, DoesSolvateFinishesSolventShell) { AtomCollection solute; solute.push_back(Atom(ElementType::Na, Position(0, 0, 0))); - int sampleRes = 18; + SoluteSolventComplex::SolventPlacementSettings testSettings1; + testSettings1.resolution = 18; + testSettings1.solventOffset = 0.0; + testSettings1.maxDistance = 7.0; + testSettings1.stepSize = 1.0; + testSettings1.numRotamers = 3; + testSettings1.strategicSolv = false; + testSettings1.coverageThreshold = 1.0; + int numShells = 2; - auto solventComplex = - SoluteSolventComplex::solvateShells(solute, solute.size(), water, numShells, 5, sampleRes, 0, 7, 1, 3, false); + auto solventComplex = SoluteSolventComplex::solvateShells(solute, solute.size(), water, numShells, 5, testSettings1); auto soluteComplex = solute + SoluteSolventComplex::mergeSolventShellVector(solventComplex); - // calculate solvent vector with reduced coverage - auto reducedSolventComplex = - SoluteSolventComplex::solvateShells(solute, solute.size(), water, numShells, 5, sampleRes, 0, 7, 1, 3, false, 0.75); + SoluteSolventComplex::SolventPlacementSettings testSettings2 = testSettings1; + testSettings2.coverageThreshold = 0.75; - auto soluteComplexSurface = - MolecularSurface::getVisibleMolecularSurface(soluteComplex, solute.size(), solute.size() + 18 * 3, sampleRes); + auto reducedSolventComplex = SoluteSolventComplex::solvateShells(solute, solute.size(), water, numShells, 5, testSettings2); + auto soluteComplexSurface = MolecularSurface::getVisibleMolecularSurface( + soluteComplex, solute.size(), solute.size() + 18 * 3, testSettings1.resolution); // first solvent shell is full after addition of 18 water molecules ASSERT_THAT(soluteComplexSurface.size(), 0); @@ -363,11 +432,18 @@ TEST(SoluteSolventComplexTest, DoesGiveSolventVectorGiveCorrectStructure) { AtomCollection solute; solute.push_back(Atom(ElementType::Na, Position(0, 0, 0))); - int sampleRes = 18; + SoluteSolventComplex::SolventPlacementSettings testSettings; + testSettings.resolution = 18; + testSettings.solventOffset = 0.0; + testSettings.maxDistance = 7.0; + testSettings.stepSize = 1.0; + testSettings.numRotamers = 3; + testSettings.strategicSolv = false; + testSettings.coverageThreshold = 1.0; + int numShells = 1; - auto solventComplex = - SoluteSolventComplex::solvateShells(solute, solute.size(), water, numShells, 5, sampleRes, 0, 7, 1, 3, false); + auto solventComplex = SoluteSolventComplex::solvateShells(solute, solute.size(), water, numShells, 5, testSettings); auto soluteComplex = solute + SoluteSolventComplex::mergeSolventShellVector(solventComplex); auto solventSizeVector = SoluteSolventComplex::transferSolventShellVector(solventComplex); @@ -380,8 +456,8 @@ TEST(SoluteSolventComplexTest, DoesGiveSolventVectorGiveCorrectStructure) { Core::Log log = Core::Log::silent(); // give solvent shell vector from soluteComplex with threshold of 100% - auto givenSolventShellVector = - SoluteSolventComplex::giveSolventShellVector(soluteComplex, solute.size(), solventSizeVector, sampleRes, log); + auto givenSolventShellVector = SoluteSolventComplex::giveSolventShellVector( + soluteComplex, solute.size(), solventSizeVector, testSettings.resolution, log); // solvent vector obtained by giveSolventShellVector function has to be the same as solventComplex ASSERT_THAT(givenSolventShellVector.size(), solventComplex.size() + 1); @@ -390,7 +466,7 @@ TEST(SoluteSolventComplexTest, DoesGiveSolventVectorGiveCorrectStructure) { assertPositionCollectionsEqual(givenSolventShellVector[0][17].getPositions(), solventComplex[0][17].getPositions()); auto reducedSolventShellVector = SoluteSolventComplex::giveSolventShellVector( - soluteComplex, solute.size(), solventSizeVector, sampleRes, log, true, 0.90); + soluteComplex, solute.size(), solventSizeVector, testSettings.resolution, log, true, 0.90); // reduced solvent shell vector due to threshold ASSERT_THAT(reducedSolventShellVector.size(), solventComplex.size() + 1); @@ -400,6 +476,107 @@ TEST(SoluteSolventComplexTest, DoesGiveSolventVectorGiveCorrectStructure) { assertPositionCollectionsEqual(reducedSolventShellVector[1][0].getPositions(), solventComplex[0][17].getPositions()); } +TEST(SoluteSolventComplexTest, DoesAdditionOfMixedSolventsWork) { + // Formaldehyde + PositionCollection methanalPC(4, 3); + methanalPC.row(0) = Position(2, 0, 0); + methanalPC.row(1) = Position(0, 0, 0); + methanalPC.row(2) = Position(-1.8, 1.5, 0); + methanalPC.row(3) = Position(-1.8, -1.5, 0); + + ElementTypeCollection methanalEC(4); + methanalEC.at(0) = ElementType::O; + methanalEC.at(1) = ElementType::C; + methanalEC.at(2) = ElementType::H; + methanalEC.at(3) = ElementType::H; + + AtomCollection methanal(methanalEC, methanalPC); + // Water + ElementTypeCollection waterEC(3); + waterEC.at(0) = ElementType::O; + waterEC.at(1) = ElementType::H; + waterEC.at(2) = ElementType::H; + + PositionCollection waterPC0(3, 3); + waterPC0.row(0) = Position(0, 0, 0); + waterPC0.row(1) = Position(-1.8, 1.25, 0); + waterPC0.row(2) = Position(-1.8, -1.25, 0); + + AtomCollection water(waterEC, waterPC0); + + // Argon + AtomCollection argon; + argon.push_back(Atom(ElementType::Ar, Position(0, 0, 0))); + + std::vector solventList = {methanal, water, argon}; + + AtomCollection solute; + solute.push_back(Atom(ElementType::Na, Position(0, 0, 0))); + + SoluteSolventComplex::SolventPlacementSettings testSettings; + + testSettings.maxDistance = 2.5; + testSettings.stepSize = 1.5; + + // Test number of solvents + int numSolvents = 19; + auto solventComplexTuple = SoluteSolventComplex::solvateMix(solute, solute.size(), solventList, + std::vector{2, 6, 11}, numSolvents, 42, testSettings); + auto solventComplex = std::get<0>(solventComplexTuple); + auto solventIndices = std::get<1>(solventComplexTuple); + + auto soluteComplex = solute + SoluteSolventComplex::mergeSolventShellVector(solventComplex); + auto flatIndices = SoluteSolventComplex::mergeSolventShellIndices(solventIndices); + + // Check for correct size of cluster + ASSERT_THAT(soluteComplex.size(), 1 + 2 * 4 + 6 * 3 + 11 * 1); + ASSERT_THAT(flatIndices.size(), numSolvents); + // Matching sizes of shell and solvent sizes + std::vector matchShellSize = {17, 2}; + // Sizes of Indices {4, 3, 1} + std::vector> matchSolventSize = {{3, 3, 1, 1, 1, 1, 4, 1, 1, 1, 3, 4, 1, 3, 3, 1, 1}, {3, 1}}; + std::vector> matchSolventIndices = {{1, 1, 2, 2, 2, 2, 0, 2, 2, 2, 1, 0, 2, 1, 1, 2, 2}, {1, 2}}; + for (int i = 0; i < int(solventComplex.size()); i++) { + ASSERT_THAT(solventComplex[i].size(), matchShellSize[i]); + ASSERT_THAT(solventIndices[i].size(), matchShellSize[i]); + for (int j = 0; j < int(solventComplex[i].size()); j++) { + ASSERT_THAT(solventComplex[i][j].size(), matchSolventSize[i][j]); + ASSERT_THAT(solventIndices[i][j], matchSolventIndices[i][j]); + } + } + + // Test number of shells + int numShells = 2; + SoluteSolventComplex::SolventPlacementSettings testSettingsShell = testSettings; + testSettingsShell.coverageThreshold = 0.40; + + auto solventComplexShellTuple = SoluteSolventComplex::solvateShellsMix( + solute, solute.size(), solventList, std::vector{11, 6, 2}, numShells, 42, testSettingsShell); + + auto solventComplexShell = std::get<0>(solventComplexShellTuple); + auto solventComplexShellIndices = std::get<1>(solventComplexShellTuple); + + auto soluteComplexShell = solute + SoluteSolventComplex::mergeSolventShellVector(solventComplexShell); + + ASSERT_THAT(solventComplexShell.size(), numShells + 1); + ASSERT_THAT(solventComplexShellIndices.size(), numShells + 1); + // Matching sizes of shell and solvent sizes + std::vector matchShellSize2 = {9, 22}; + std::vector> matchSolventSize2 = { + {4, 4, 3, 4, 4, 1, 3, 4, 4}, {3, 4, 3, 3, 3, 3, 4, 4, 3, 3, 3, 4, 4, 4, 3, 4, 3, 3, 3, 4, 4, 3}}; + std::vector> matchSolventIndices2 = { + {0, 0, 1, 0, 0, 2, 1, 0, 0}, {1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1}}; + + for (int i = 0; i < int(solventComplexShell.size()) - 1; i++) { + ASSERT_THAT(solventComplexShellIndices[i].size(), matchShellSize2[i]); + ASSERT_THAT(solventComplexShell[i].size(), matchShellSize2[i]); + for (int j = 0; j < int(solventComplexShell[i].size()); j++) { + ASSERT_THAT(solventComplexShell[i][j].size(), matchSolventSize2[i][j]); + ASSERT_THAT(solventComplexShellIndices[i][j], matchSolventIndices2[i][j]); + } + } +} + } /* namespace Tests */ } /* namespace Utils */ } /* namespace Scine */ diff --git a/src/Utils/Tests/UniversalSettings/Settings.cpp b/src/Utils/Tests/UniversalSettings/Settings.cpp index 6f1433c..5c071b6 100644 --- a/src/Utils/Tests/UniversalSettings/Settings.cpp +++ b/src/Utils/Tests/UniversalSettings/Settings.cpp @@ -34,6 +34,54 @@ TEST(Settings, InvalidValueExceptions) { ASSERT_THROW(settings.throwIncorrectSettings(), InvalidSettingsException); } +TEST(Settings, DivergingKeysAndEquality) { + using namespace UniversalSettings; + + DescriptorCollection descriptors; + IntDescriptor multiplicity("spin multiplicity"); + IntDescriptor charge("charge"); + StringDescriptor method("method"); + + descriptors.push_back("multiplicity", multiplicity); + descriptors.push_back("charge", charge); + descriptors.push_back("method", method); + + ValueCollection values; + values.addInt("multiplicity", 1); + values.addInt("charge", 0); + values.addString("method", "bla"); + + ValueCollection otherValues; + otherValues.addInt("multiplicity", 3); + otherValues.addInt("charge", 0); + otherValues.addString("method", "blubb"); + + ValueCollection duplicateValues; + duplicateValues.addInt("multiplicity", 1); + duplicateValues.addInt("charge", 0); + duplicateValues.addString("method", "bla"); + + Settings settings{values, descriptors}; + Settings otherSettings{otherValues, descriptors}; + Settings duplicateSettings{duplicateValues, descriptors}; + + ASSERT_TRUE(settings.valid()); + ASSERT_TRUE(otherSettings.valid()); + + ASSERT_FALSE(values == otherValues); + ASSERT_FALSE(duplicateValues == otherValues); + ASSERT_TRUE(values == duplicateValues); + + auto otherDiverging = values.getDivergingKeys(otherValues); + ASSERT_TRUE(otherDiverging.size() == 2); + ASSERT_TRUE(values.getDivergingKeys(otherValues, true).size() == 1); + ASSERT_TRUE(std::find(otherDiverging.begin(), otherDiverging.end(), "multiplicity") != otherDiverging.end()); + ASSERT_TRUE(std::find(otherDiverging.begin(), otherDiverging.end(), "method") != otherDiverging.end()); + ASSERT_TRUE(std::find(otherDiverging.begin(), otherDiverging.end(), "charge") == otherDiverging.end()); + + ASSERT_TRUE(values.getDivergingKeys(duplicateValues).empty()); +} + TEST(Settings, CoerceStringCases) { using namespace UniversalSettings; diff --git a/src/Utils/Utils/Bonds/BondOrderCollection.h b/src/Utils/Utils/Bonds/BondOrderCollection.h index cc2c461..b68ba3c 100644 --- a/src/Utils/Utils/Bonds/BondOrderCollection.h +++ b/src/Utils/Utils/Bonds/BondOrderCollection.h @@ -121,6 +121,24 @@ class BondOrderCollection { rangeCheck(i, j); return bondOrderMatrix_.coeff(i, j); } + /** + * @brief Get bond partner of an index based on a minimum bond threshold + * @tparam Index + * @param i The index of whom we want the bond partner. + * @param threshold The threshold above which a bond exists. + * @return std::vector The list of bond partners. + */ + template + std::vector getBondPartners(Index i, double threshold = 0.5) const { + rangeCheck(i, i); + std::vector result; + for (Index j = 0; j < getSystemSize(); ++j) { + if (bondOrderMatrix_.coeff(i, j) > threshold) { + result.push_back(j); + } + } + return result; + } /** * @brief Checks whether two bond order collections are approximately equal diff --git a/src/Utils/Utils/CalculatorBasics/CalculationRoutines.h b/src/Utils/Utils/CalculatorBasics/CalculationRoutines.h index bb1a333..67da914 100644 --- a/src/Utils/Utils/CalculatorBasics/CalculationRoutines.h +++ b/src/Utils/Utils/CalculatorBasics/CalculationRoutines.h @@ -139,6 +139,21 @@ inline std::pair splitIntoMethodAndDispersion(const st if (input.empty()) { return std::make_pair("", ""); } + // check for exceptions + // if input contains one of these as a substring + // we return the input as method and empty dispersion + std::vector exceptions = { + "PNO-CC", + "HF-3C", + "PBEH-3C", + "B97-3C", + }; + std::string inputCopy(input.size(), 0); + std::transform(input.begin(), input.end(), inputCopy.begin(), [](unsigned char c) { return std::toupper(c); }); + if (std::any_of(exceptions.begin(), exceptions.end(), + [&, inputCopy](const std::string& e) { return inputCopy.find(e) != std::string::npos; })) + return std::make_pair(input, ""); + std::string segment; std::vector segments; std::stringstream ss(input); diff --git a/src/Utils/Utils/DataStructures/PeriodicBoundaries.cpp b/src/Utils/Utils/DataStructures/PeriodicBoundaries.cpp index 0b21130..62876f5 100644 --- a/src/Utils/Utils/DataStructures/PeriodicBoundaries.cpp +++ b/src/Utils/Utils/DataStructures/PeriodicBoundaries.cpp @@ -15,8 +15,8 @@ PeriodicBoundaries::PeriodicBoundaries(double cubeLength, const std::string& per : PeriodicBoundaries(Eigen::Matrix3d::Identity() * cubeLength, periodicity) { } -PeriodicBoundaries::PeriodicBoundaries(const PeriodicBoundaries& rhs, const std::string& periodicity) - : PeriodicBoundaries(rhs.getCellMatrix(), periodicity) { +PeriodicBoundaries::PeriodicBoundaries(const PeriodicBoundaries& rhs) + : PeriodicBoundaries(rhs.getCellMatrix(), rhs.getPeriodicityString()) { } PeriodicBoundaries::PeriodicBoundaries(Eigen::Matrix3d matrix, const std::string& periodicity) @@ -301,6 +301,21 @@ bool PeriodicBoundaries::isWithinCell(const PositionCollection& positions) const return true; } +void PeriodicBoundaries::canonicalize() { + auto canonicPbc = PeriodicBoundaries(getLengths(), getAngles(), true, true, getPeriodicityString()); + setCellMatrix(canonicPbc.getCellMatrix()); +} + +Eigen::Matrix3d PeriodicBoundaries::getCanonicalizationRotationMatrix() const { + auto canonicPbc = PeriodicBoundaries(getLengths(), getAngles(), true, true, getPeriodicityString()); + if (canonicPbc == *this) { + // we are already canonic, avoid numeric instability, just return I + return Eigen::Matrix3d::Identity(); + } + Eigen::Matrix3d rotationMatrix = this->_invMatrix * canonicPbc.getCellMatrix(); + return rotationMatrix; +} + PeriodicBoundaries& PeriodicBoundaries::operator=(const PeriodicBoundaries& rhs) { this->setCellMatrix(rhs._matrix); this->setPeriodicity(rhs._periodicity); diff --git a/src/Utils/Utils/DataStructures/PeriodicBoundaries.h b/src/Utils/Utils/DataStructures/PeriodicBoundaries.h index 5b5dac4..e7213c1 100644 --- a/src/Utils/Utils/DataStructures/PeriodicBoundaries.h +++ b/src/Utils/Utils/DataStructures/PeriodicBoundaries.h @@ -41,7 +41,7 @@ namespace Utils { class PeriodicBoundaries { public: explicit PeriodicBoundaries(double cubeLength = 1.0, const std::string& periodicity = "xyz"); - PeriodicBoundaries(const PeriodicBoundaries& rhs, const std::string& periodicity = "xyz"); + PeriodicBoundaries(const PeriodicBoundaries& rhs); explicit PeriodicBoundaries(Eigen::Matrix3d matrix, const std::string& periodicity = "xyz"); explicit PeriodicBoundaries(const Eigen::Vector3d& lengths, const Eigen::Vector3d& angles, bool isBohr = true, @@ -213,10 +213,26 @@ class PeriodicBoundaries { */ bool isWithinCell(const PositionCollection& positions) const; + void canonicalize(); + + Eigen::Matrix3d getCanonicalizationRotationMatrix() const; + inline bool isOrthoRhombic(double eps = 1e-2) const { return std::fabs(_alpha - 90.0) < eps && std::fabs(_beta - 90.0) < eps && std::fabs(_gamma - 90.0) < eps; }; + inline Eigen::Vector3d getLengths() const { + Eigen::Vector3d result; + result << _aNorm, _bNorm, _cNorm; + return result; + }; + + inline Eigen::Vector3d getAngles() const { + Eigen::Vector3d result; + result << _alpha, _beta, _gamma; + return result; + }; + inline std::string getPeriodicBoundariesString(const std::string& delimiter = ",") const { return std::to_string(_aNorm) + delimiter + std::to_string(_bNorm) + delimiter + std::to_string(_cNorm) + delimiter + std::to_string(_alpha) + delimiter + std::to_string(_beta) + delimiter + std::to_string(_gamma) + diff --git a/src/Utils/Utils/ExternalQC/ExternalProgram.cpp b/src/Utils/Utils/ExternalQC/ExternalProgram.cpp index cd34e16..d98cf3b 100644 --- a/src/Utils/Utils/ExternalQC/ExternalProgram.cpp +++ b/src/Utils/Utils/ExternalQC/ExternalProgram.cpp @@ -34,6 +34,10 @@ void ExternalProgram::createWorkingDirectory() { } } +void ExternalProgram::setErrorOutFile(const std::string& filename) { + _errorOut = filename; +} + void ExternalProgram::executeCommand(const std::string& command) const { executeCommand(command, "", ""); } @@ -55,20 +59,30 @@ int ExternalProgram::executeCommandImpl(const std::string& command, const std::s const std::string& outputFile) const { bool hasInput = !inputFile.empty(); bool hasOutput = !outputFile.empty(); + bool hasError = !_errorOut.empty(); auto workingDirectory = bp::start_dir(workingDirectory_); if (workingDirectory_.empty()) { workingDirectory = bp::start_dir(FilesystemHelpers::currentDirectory()); } + if (hasInput && hasOutput && hasError) { + return bp::system(command, bp::std_out > outputFile, bp::std_err > _errorOut, bp::std_in < inputFile, workingDirectory); + } if (hasInput && hasOutput) { return bp::system(command, bp::std_out > outputFile, bp::std_in < inputFile, workingDirectory); } + if (hasInput && hasError) { + return bp::system(command, bp::std_err > _errorOut, bp::std_in < inputFile, workingDirectory); + } + if (hasOutput && hasError) { + return bp::system(command, bp::std_out > outputFile, bp::std_err > _errorOut, workingDirectory); + } if (hasInput) { return bp::system(command, bp::std_in < inputFile, workingDirectory); } - if (hasOutput) { - return bp::system(command, bp::std_out > outputFile); + if (hasError) { + return bp::system(command, bp::std_err > _errorOut, workingDirectory); } if (hasOutput) { return bp::system(command, bp::std_out > outputFile, workingDirectory); diff --git a/src/Utils/Utils/ExternalQC/ExternalProgram.h b/src/Utils/Utils/ExternalQC/ExternalProgram.h index 6d439c6..fbd360d 100644 --- a/src/Utils/Utils/ExternalQC/ExternalProgram.h +++ b/src/Utils/Utils/ExternalQC/ExternalProgram.h @@ -51,10 +51,15 @@ class ExternalProgram { * @brief Generate filename by prepending the working directory. */ std::string generateFullFilename(const std::string& filename) const; + /** + * @brief Setter to write error out to specific file. + */ + void setErrorOutFile(const std::string& filename); private: int executeCommandImpl(const std::string& command, const std::string& inputFile, const std::string& outputFile) const; std::string workingDirectory_{}; + std::string _errorOut; }; } // namespace ExternalQC diff --git a/src/Utils/Utils/ExternalQC/Orca/OrcaCalculator.h b/src/Utils/Utils/ExternalQC/Orca/OrcaCalculator.h index d88cd9b..7ee9f55 100644 --- a/src/Utils/Utils/ExternalQC/Orca/OrcaCalculator.h +++ b/src/Utils/Utils/ExternalQC/Orca/OrcaCalculator.h @@ -162,10 +162,11 @@ class OrcaCalculator final : public CloneInterface availableSolvationModels_ = std::vector{"cpcm", "smd"}; - const std::vector availableMethodFamilies_ = std::vector{"DFT", "HF", "CC", "HF-3C", "PBEH-3C"}; + const std::vector availableMethodFamilies_ = + std::vector{"DFT", "HF", "CC", "HF-3C", "PBEH-3C", "B97-3C"}; // incomplete list of methods with no analytical Hessian, add if you find a new one const std::vector numFreqMethods_ = - std::vector{"M06", "DLPNO-CCSD(T)", "DLPNO-CCSD", "HF-3C", "PBEH-3C"}; + std::vector{"M06", "DLPNO-CCSD(T)", "DLPNO-CCSD", "HF-3C", "PBEH-3C", "B97-3C"}; }; } // namespace ExternalQC diff --git a/src/Utils/Utils/ExternalQC/SettingsNames.h b/src/Utils/Utils/ExternalQC/SettingsNames.h index ae7e286..21e7ef4 100644 --- a/src/Utils/Utils/ExternalQC/SettingsNames.h +++ b/src/Utils/Utils/ExternalQC/SettingsNames.h @@ -40,6 +40,7 @@ static constexpr const char* additionalOutputFile = "additional_output_file"; static constexpr const char* gaussianFilenameBase = "gaussian_filename_base"; // turbomole static constexpr const char* scfOrbitalShift = "scf_orbitalshift"; +static constexpr const char* scfDampingValue = "scf_damping_value"; } // namespace SettingsNames } // namespace ExternalQC diff --git a/src/Utils/Utils/ExternalQC/Turbomole/TurbomoleCalculatorSettings.h b/src/Utils/Utils/ExternalQC/Turbomole/TurbomoleCalculatorSettings.h index c7c155f..ae012bb 100644 --- a/src/Utils/Utils/ExternalQC/Turbomole/TurbomoleCalculatorSettings.h +++ b/src/Utils/Utils/ExternalQC/Turbomole/TurbomoleCalculatorSettings.h @@ -36,6 +36,7 @@ class TurbomoleCalculatorSettings : public Scine::Utils::Settings { void addTemperature(UniversalSettings::DescriptorCollection& settings); void addElectronicTemperature(UniversalSettings::DescriptorCollection& settings); void addScfDamping(UniversalSettings::DescriptorCollection& settings); + void addScfDampingValue(UniversalSettings::DescriptorCollection& settings); void addScfOrbitalShift(UniversalSettings::DescriptorCollection& settings); void addSolvent(UniversalSettings::DescriptorCollection& settings); void addSolvation(UniversalSettings::DescriptorCollection& settings); @@ -57,6 +58,7 @@ class TurbomoleCalculatorSettings : public Scine::Utils::Settings { addBaseWorkingDirectory(_fields); addTemperature(_fields); addScfDamping(_fields); + addScfDampingValue(_fields); addScfOrbitalShift(_fields); addElectronicTemperature(_fields); addSolvent(_fields); @@ -140,15 +142,17 @@ inline void TurbomoleCalculatorSettings::addTemperature(UniversalSettings::Descr } inline void TurbomoleCalculatorSettings::addScfDamping(UniversalSettings::DescriptorCollection& settings) { - Utils::UniversalSettings::OptionListDescriptor scfDamping("Specify SCF damping (low/medium/high)."); - scfDamping.addOption("default"); - scfDamping.addOption("low"); - scfDamping.addOption("medium"); - scfDamping.addOption("high"); - scfDamping.setDefaultOption("default"); + Utils::UniversalSettings::BoolDescriptor scfDamping("Enable stronger SCF damping (true/false)."); + scfDamping.setDefaultValue(false); settings.push_back(Utils::SettingsNames::scfDamping, std::move(scfDamping)); } +inline void TurbomoleCalculatorSettings::addScfDampingValue(UniversalSettings::DescriptorCollection& settings) { + Utils::UniversalSettings::DoubleDescriptor scfDampingValue("Specify exact SCF damping value to be used."); + scfDampingValue.setDefaultValue(0.5); + settings.push_back(SettingsNames::scfDampingValue, std::move(scfDampingValue)); +} + inline void TurbomoleCalculatorSettings::addScfOrbitalShift(UniversalSettings::DescriptorCollection& settings) { Utils::UniversalSettings::DoubleDescriptor scfOrbitalShift( "Shift closed shells to lower energies to aid convergence."); diff --git a/src/Utils/Utils/ExternalQC/Turbomole/TurbomoleInputFileCreator.cpp b/src/Utils/Utils/ExternalQC/Turbomole/TurbomoleInputFileCreator.cpp index 1e59cd7..53f3bfa 100644 --- a/src/Utils/Utils/ExternalQC/Turbomole/TurbomoleInputFileCreator.cpp +++ b/src/Utils/Utils/ExternalQC/Turbomole/TurbomoleInputFileCreator.cpp @@ -109,13 +109,16 @@ void TurbomoleInputFileCreator::prepareDefineSession(const Settings& settings, c if (!methodInput.second.empty()) { // make sure that dispersion is uppercase only std::for_each(methodInput.second.begin(), methodInput.second.end(), [](char& c) { c = ::toupper(c); }); - auto it = std::find(availableD3Params_.begin(), availableD3Params_.end(), methodInput.second); - if (it - availableD3Params_.begin() == 0) { + auto it = std::find(availableDispersionParams_.begin(), availableDispersionParams_.end(), methodInput.second); + if (it - availableDispersionParams_.begin() == 0) { out << "dsp\non\n\n"; } - else if (it - availableD3Params_.begin() == 1) { + else if (it - availableDispersionParams_.begin() == 1) { out << "dsp\nbj\n\n"; } + else if (it - availableDispersionParams_.begin() == 2) { + out << "dsp\nd4\n\n"; + } else { throw std::runtime_error("Invalid dispersion correction!"); } @@ -139,7 +142,7 @@ void TurbomoleInputFileCreator::checkAndUpdateControlFile(const Settings& settin bool functionalIsCorrect = false; double scfConvCriterion = settings.getDouble(Utils::SettingsNames::selfConsistenceCriterion); - auto scfDamping = settings.getString(Utils::SettingsNames::scfDamping); + bool scfDamping = settings.getBool(Utils::SettingsNames::scfDamping); auto scfOrbitalShift = settings.getDouble(SettingsNames::scfOrbitalShift); auto solvent = settings.getString(Utils::SettingsNames::solvent); @@ -166,10 +169,9 @@ void TurbomoleInputFileCreator::checkAndUpdateControlFile(const Settings& settin helper.mapDftFunctionalToTurbomoleStringRepresentation(method); while (std::getline(in, line)) { - auto it = availableDampingParameter_.find(scfDamping); - if ((line.find("$scfdamp") != std::string::npos) && (it != availableDampingParameter_.end())) { - out << "$scfdamp start=" << std::fixed << std::setprecision(3) << it->second[0] << " step=" << std::fixed - << std::setprecision(2) << it->second[1] << " min=" << std::fixed << std::setprecision(2) << it->second[2] << "\n"; + if ((line.find("$scfdamp") != std::string::npos) && scfDamping) { + double dampingValue = settings.getDouble(SettingsNames::scfDampingValue); + out << "$scfdamp start=" << std::fixed << std::setprecision(3) << dampingValue << " step=0.05 min=0.10\n"; } else if ((line.find("scforbitalshift") != std::string::npos)) { out << "$scforbitalshift closedshell=" << scfOrbitalShift << "\n"; diff --git a/src/Utils/Utils/ExternalQC/Turbomole/TurbomoleInputFileCreator.h b/src/Utils/Utils/ExternalQC/Turbomole/TurbomoleInputFileCreator.h index c25df7c..8da4e4f 100644 --- a/src/Utils/Utils/ExternalQC/Turbomole/TurbomoleInputFileCreator.h +++ b/src/Utils/Utils/ExternalQC/Turbomole/TurbomoleInputFileCreator.h @@ -52,22 +52,45 @@ class TurbomoleInputFileCreator { std::string defineExecutableBase_ = "define"; // maps the available solvents to the corresponding dielectric constants and radii of the solvent molecules in // Angstrom. Constants have been taken from the ADF documentation (https://www.scm.com/doc/ADF/Input/COSMO.html, last - // visited Sept 06, 2021). + // visited March 23, 2022). std::unordered_map> availableSolventModels_ = { - {"acetone", std::make_pair(21.7, 3.08)}, {"ammonia", std::make_pair(16.9, 2.24)}, - {"benzene", std::make_pair(2.3, 3.28)}, {"chloroform", std::make_pair(4.8, 3.17)}, - {"dmso", std::make_pair(46.7, 3.04)}, {"ethanol", std::make_pair(24.5, 2.85)}, - {"hexane", std::make_pair(1.88, 3.74)}, {"h2o", std::make_pair(78.4, 1.93)}, - {"methanol", std::make_pair(32.6, 2.53)}, {"nitrobenzene", std::make_pair(34.8, 3.44)}, - {"thf", std::make_pair(7.58, 3.18)}, {"toluene", std::make_pair(2.38, 3.48)}, - {"water", std::make_pair(78.4, 1.93)}, {"isopropanol", std::make_pair(19.9, 3.12)}}; + {"aceticacid", std::make_pair(6.19, 2.83)}, + {"acetonitrile", std::make_pair(37.5, 2.76)}, + {"aniline", std::make_pair(6.8, 3.31)}, + {"benzylalcohol", std::make_pair(13.1, 3.45)}, + {"bromoform", std::make_pair(4.3, 3.26)}, + {"butanol", std::make_pair(17.5, 3.31)}, + {"isobutanol", std::make_pair(17.9, 3.33)}, + {"tertbutanol", std::make_pair(12.4, 3.35)}, + {"carbondisulfide", std::make_pair(2.6, 2.88)}, + {"carbontetrachloride", std::make_pair(2.2, 3.37)}, + {"cyclohexane", std::make_pair(2.0, 3.5)}, + {"cyclohexanone", std::make_pair(15.0, 3.46)}, + {"dichlorobenzene", std::make_pair(9.8, 3.54)}, + {"diethylether", std::make_pair(4.34, 3.46)}, + {"dioxane", std::make_pair(2.2, 3.24)}, + {"dmfa", std::make_pair(37.0, 3.13)}, + {"ethylacetate", std::make_pair(6.02, 3.39)}, + {"dichloroethane", std::make_pair(10.66, 3.15)}, + {"ethyleneglycol", std::make_pair(37.7, 2.81)}, + {"formicacid", std::make_pair(58.5, 2.47)}, + {"acetone", std::make_pair(20.7, 3.08)}, + {"ammonia", std::make_pair(16.9, 2.24)}, + {"benzene", std::make_pair(2.3, 3.28)}, + {"chloroform", std::make_pair(4.8, 3.17)}, + {"dmso", std::make_pair(46.7, 3.04)}, + {"ethanol", std::make_pair(24.55, 2.85)}, + {"hexane", std::make_pair(1.88, 3.74)}, + {"h2o", std::make_pair(78.39, 1.93)}, + {"methanol", std::make_pair(32.6, 2.53)}, + {"nitrobenzene", std::make_pair(34.8, 3.44)}, + {"thf", std::make_pair(7.58, 3.18)}, + {"toluene", std::make_pair(2.38, 3.48)}, + {"water", std::make_pair(78.39, 1.93)}, + {"isopropanol", std::make_pair(19.9, 3.12)}}; + + const std::vector availableDispersionParams_ = {"D3", "D3BJ", "D4"}; - const std::vector availableD3Params_ = {"D3", "D3BJ"}; - // Sets of damping parameters to aid SCF convergence. The old - // Fock operator is added to the current one with a specific weight (first parameter). This weight is reduced in - // specific steps (second parameter) until a minimum is reached (third parameter). - std::unordered_map> availableDampingParameter_{ - {"low", {1.5, 0.05, 0.1}}, {"medium", {5.0, 0.1, 0.5}}, {"high", {8.5, 0.1, 0.5}}}; TurbomoleFiles files_; }; diff --git a/src/Utils/Utils/Geometry/ElementData.cpp b/src/Utils/Utils/Geometry/ElementData.cpp index 531e3d3..fb4d96b 100644 --- a/src/Utils/Utils/Geometry/ElementData.cpp +++ b/src/Utils/Utils/Geometry/ElementData.cpp @@ -51,122 +51,129 @@ const ElementDataSingleton::ElementData& ElementDataSingleton::lookup(const Elem return data().at(e); } +double ElementDataSingleton::ElementData::paulingElectronegativity() const { + if (d_paulingElectronegativity > -1) { + return d_paulingElectronegativity; + } + throw DataNotAvailable(); +} + const std::unordered_map& ElementDataSingleton::data() { // clang-format off static const std::unordered_map map = { - {ElementType::none, ElementData("None", 0, 0.0, 0.0, 0.0, -1, 0, 0, 0, 0)}, - {ElementType::H, ElementData("H", 1, 1.0079, 32, 109, -1, 1, 0, 0, 0)}, - {ElementType::He, ElementData("He", 2, 4.0026, 37, 140, -1, 2, 0, 0, 0)}, - {ElementType::Li, ElementData("Li", 3, 6.941, 130, 182, -1, 1, 0, 0, 0)}, - {ElementType::Be, ElementData("Be", 4, 9.0122, 99, 153, -1, 2, 0, 0, 0)}, - {ElementType::B, ElementData("B", 5, 10.811, 84, 192, -1, 2, 1, 0, 0)}, - {ElementType::C, ElementData("C", 6, 12.011, 75, 170, -1, 2, 2, 0, 0)}, - {ElementType::N, ElementData("N", 7, 14.007, 71, 155, -1, 2, 3, 0, 0)}, - {ElementType::O, ElementData("O", 8, 15.999, 64, 152, -1, 2, 4, 0, 0)}, - {ElementType::F, ElementData("F", 9, 18.988, 60, 147, -1, 2, 5, 0, 0)}, - {ElementType::Ne, ElementData("Ne", 10, 20.180, 62, 154, -1, 2, 6, 0, 0)}, - {ElementType::Na, ElementData("Na", 11, 22.990, 160, 227, -1, 1, 0, 0, 0)}, - {ElementType::Mg, ElementData("Mg", 12, 24.305, 140, 173, -1, 2, 0, 0, 0)}, - {ElementType::Al, ElementData("Al", 13, 26.982, 124, 184, -1, 2, 1, 0, 0)}, - {ElementType::Si, ElementData("Si", 14, 28.086, 114, 210, -1, 2, 2, 0, 0)}, - {ElementType::P, ElementData("P", 15, 30.974, 109, 180, -1, 2, 3, 0, 0)}, - {ElementType::S, ElementData("S", 16, 32.065, 104, 180, -1, 2, 4, 0, 0)}, - {ElementType::Cl, ElementData("Cl", 17, 35.453, 100, 175, -1, 2, 5, 0, 0)}, - {ElementType::Ar, ElementData("Ar", 18, 39.948, 101, 188, -1, 2, 6, 0, 0)}, - {ElementType::K, ElementData("K", 19, 39.098, 200, 275, -1, 1, 0, 0, 0)}, - {ElementType::Ca, ElementData("Ca", 20, 40.078, 174, 231, -1, 2, 0, 0, 0)}, - {ElementType::Sc, ElementData("Sc", 21, 44.956, 159, 215, -1, 2, 0, 1, 0)}, - {ElementType::Ti, ElementData("Ti", 22, 47.867, 148, 211, -1, 2, 0, 2, 0)}, - {ElementType::V, ElementData("V", 23, 50.942, 144, 207, -1, 2, 0, 3, 0)}, - {ElementType::Cr, ElementData("Cr", 24, 51.996, 130, 206, -1, 2, 0, 4, 0)}, - {ElementType::Mn, ElementData("Mn", 25, 54.938, 129, 205, -1, 2, 0, 5, 0)}, - {ElementType::Fe, ElementData("Fe", 26, 55.938, 124, 204, -1, 2, 0, 6, 0)}, - {ElementType::Co, ElementData("Co", 27, 58.933, 118, 200, -1, 2, 0, 7, 0)}, - {ElementType::Ni, ElementData("Ni", 28, 58.693, 117, 197, -1, 2, 0, 8, 0)}, - {ElementType::Cu, ElementData("Cu", 29, 63.546, 122, 196, -1, 2, 0, 9, 0)}, - {ElementType::Zn, ElementData("Zn", 30, 65.38, 120, 201, -1, 2, 0, 10, 0)}, - {ElementType::Ga, ElementData("Ga", 31, 69.723, 123, 187, -1, 2, 1, 10, 0)}, - {ElementType::Ge, ElementData("Ge", 32, 72.64, 120, 211, -1, 2, 2, 10, 0)}, - {ElementType::As, ElementData("As", 33, 74.922, 120, 185, -1, 2, 3, 10, 0)}, - {ElementType::Se, ElementData("Se", 34, 78.96, 118, 190, -1, 2, 4, 10, 0)}, - {ElementType::Br, ElementData("Br", 35, 79.904, 117, 185, -1, 2, 5, 10, 0)}, - {ElementType::Kr, ElementData("Kr", 36, 83.798, 116, 202, -1, 2, 6, 10, 0)}, - {ElementType::Rb, ElementData("Rb", 37, 83.468, 215, 303, -1, 1, 0, 0, 0)}, - {ElementType::Sr, ElementData("Sr", 38, 87.62, 190, 249, -1, 2, 0, 0, 0)}, - {ElementType::Y, ElementData("Y", 39, 88.906, 176, 232, -1, 2, 0, 1, 0)}, - {ElementType::Zr, ElementData("Zr", 40, 91.224, 164, 223, -1, 2, 0, 2, 0)}, - {ElementType::Nb, ElementData("Nb", 41, 92.906, 156, 218, -1, 2, 0, 3, 0)}, - {ElementType::Mo, ElementData("Mo", 42, 95.96, 146, 217, -1, 2, 0, 4, 0)}, - {ElementType::Tc, ElementData("Tc", 43, 98.91, 138, 216, -1, 2, 0, 5, 0)}, - {ElementType::Ru, ElementData("Ru", 44, 101.07, 136, 213, -1, 2, 0, 6, 0)}, - {ElementType::Rh, ElementData("Rh", 45, 102.91, 134, 210, -1, 2, 0, 7, 0)}, - {ElementType::Pd, ElementData("Pd", 46, 106.42, 130, 210, -1, 2, 0, 8, 0)}, - {ElementType::Ag, ElementData("Ag", 47, 107.87, 136, 211, -1, 2, 0, 9, 0)}, - {ElementType::Cd, ElementData("Cd", 48, 112.41, 140, 218, -1, 2, 0, 10, 0)}, - {ElementType::In, ElementData("In", 49, 114.82, 142, 193, -1, 2, 1, 10, 0)}, - {ElementType::Sn, ElementData("Sn", 50, 118.71, 140, 217, -1, 2, 2, 10, 0)}, - {ElementType::Sb, ElementData("Sb", 51, 121.76, 140, 206, -1, 2, 3, 10, 0)}, - {ElementType::Te, ElementData("Te", 52, 127.60, 137, 206, -1, 2, 4, 10, 0)}, - {ElementType::I, ElementData("I", 53, 126.90, 136, 198, -1, 2, 5, 10, 0)}, - {ElementType::Xe, ElementData("Xe", 54, 131.29, 136, 216, -1, 2, 6, 10, 0)}, - {ElementType::Cs, ElementData("Cs", 55, 132.91, 238, 343, -1, 1, 0, 0, 0)}, - {ElementType::Ba, ElementData("Ba", 56, 137.33, 206, 268, -1, 2, 0, 0, 0)}, - {ElementType::La, ElementData("La", 57, 138.91, 194, 243, -1, 2, 0, 0, 1)}, - {ElementType::Ce, ElementData("Ce", 58, 140.12, 184, 242, -1, 2, 0, 0, 2)}, - {ElementType::Pr, ElementData("Pr", 59, 140.91, 190, 240, -1, 2, 0, 0, 3)}, - {ElementType::Nd, ElementData("Nd", 60, 144.24, 188, 239, -1, 2, 0, 0, 4)}, - {ElementType::Pm, ElementData("Pm", 61, 146.90, 186, 238, -1, 2, 0, 0, 5)}, - {ElementType::Sm, ElementData("Sm", 62, 150.36, 185, 236, -1, 2, 0, 0, 6)}, - {ElementType::Eu, ElementData("Eu", 63, 151.96, 183, 235, -1, 2, 0, 0, 7)}, - {ElementType::Gd, ElementData("Gd", 64, 157.25, 182, 234, -1, 2, 0, 0, 8)}, - {ElementType::Tb, ElementData("Tb", 65, 158.93, 181, 233, -1, 2, 0, 0, 9)}, - {ElementType::Dy, ElementData("Dy", 66, 162.50, 180, 231, -1, 2, 0, 0, 10)}, - {ElementType::Ho, ElementData("Ho", 67, 164.93, 179, 230, -1, 2, 0, 0, 11)}, - {ElementType::Er, ElementData("Er", 68, 167.26, 177, 229, -1, 2, 0, 0, 12)}, - {ElementType::Tm, ElementData("Tm", 69, 168.93, 177, 227, -1, 2, 0, 0, 13)}, - {ElementType::Yb, ElementData("Yb", 70, 173.05, 178, 226, -1, 2, 0, 0, 14)}, - {ElementType::Lu, ElementData("Lu", 71, 174.97, 174, 224, -1, 2, 0, 1, 14)}, - {ElementType::Hf, ElementData("Hf", 72, 178.49, 164, 223, -1, 2, 0, 2, 14)}, - {ElementType::Ta, ElementData("Ta", 73, 180.95, 158, 222, -1, 2, 0, 3, 14)}, - {ElementType::W, ElementData("W", 74, 183.84, 150, 218, -1, 2, 0, 4, 14)}, - {ElementType::Re, ElementData("Re", 75, 186.21, 141, 216, -1, 2, 0, 5, 14)}, - {ElementType::Os, ElementData("Os", 76, 190.23, 136, 216, -1, 2, 0, 6, 14)}, - {ElementType::Ir, ElementData("Ir", 77, 192.22, 132, 213, -1, 2, 0, 7, 14)}, - {ElementType::Pt, ElementData("Pt", 78, 195.08, 130, 213, -1, 2, 0, 8, 14)}, - {ElementType::Au, ElementData("Au", 79, 196.97, 130, 214, -1, 2, 0, 9, 14)}, - {ElementType::Hg, ElementData("Hg", 80, 200.59, 132, 223, -1, 2, 0, 10, 14)}, - {ElementType::Tl, ElementData("Tl", 81, 204.38, 144, 196, -1, 2, 1, 10, 14)}, - {ElementType::Pb, ElementData("Pb", 82, 207.2, 145, 202, -1, 2, 2, 10, 14)}, - {ElementType::Bi, ElementData("Bi", 83, 208.98, 150, 207, -1, 2, 3, 10, 14)}, - {ElementType::Po, ElementData("Po", 84, 209.98, 142, 197, -1, 2, 4, 10, 14)}, - {ElementType::At, ElementData("At", 85, 210, 148, 202, -1, 2, 5, 10, 14)}, - {ElementType::Rn, ElementData("Rn", 86, 222, 146, 220, -1, 2, 6, 10, 14)}, - {ElementType::Fr, ElementData("Fr", 87, 223, 242, 348, -1, 1, 0, 0, 0)}, - {ElementType::Ra, ElementData("Ra", 88, 226.03, 211, 283, -1, 2, 0, 0, 0)}, - {ElementType::Ac, ElementData("Ac", 89, 227, 201, 247, -1, 2, 0, 0, 1)}, - {ElementType::Th, ElementData("Th", 90, 232.04, 190, 245, -1, 2, 0, 0, 2)}, - {ElementType::Pa, ElementData("Pa", 91, 231.04, 184, 243, -1, 2, 0, 0, 3)}, - {ElementType::U, ElementData("U", 92, 238.03, 183, 241, -1, 2, 0, 0, 4)}, - {ElementType::Np, ElementData("Np", 93, 237.05, 180, 239, -1, 2, 0, 0, 5)}, - {ElementType::Pu, ElementData("Pu", 94, 244.10, 180, 243, -1, 2, 0, 0, 6)}, - {ElementType::Am, ElementData("Am", 95, 243.10, 173, 244, -1, 2, 0, 0, 7)}, - {ElementType::Cm, ElementData("Cm", 96, 247.10, 168, 245, -1, 2, 0, 0, 8)}, - {ElementType::Bk, ElementData("Bk", 97, 247.10, 168, 244, -1, 2, 0, 0, 9)}, - {ElementType::Cf, ElementData("Cf", 98, 251.10, 168, 245, -1, 2, 0, 0, 10)}, - {ElementType::Es, ElementData("Es", 99, 254.10, 165, 245, -1, 2, 0, 0, 11)}, - {ElementType::Fm, ElementData("Fm", 100, 257.10, 167, 245, -1, 2, 0, 0, 12)}, - {ElementType::Md, ElementData("Md", 101, 258, 173, 246, -1, 2, 0, 0, 13)}, - {ElementType::No, ElementData("No", 102, 259, 176, 246, -1, 2, 0, 0, 14)}, - {ElementType::Lr, ElementData("Lr", 103, 262, 161, 246, -1, 2, 0, 1, 14)}, - {ElementType::Rf, ElementData("Rf", 104, 261, 157, -1, -1, 2, 0, 2, 14)}, - {ElementType::Db, ElementData("Db", 105, 262, 149, -1, -1, 2, 0, 3, 14)}, - {ElementType::Sg, ElementData("Sg", 106, 266, 143, -1, -1, 2, 0, 4, 14)}, - {ElementType::Bh, ElementData("Bh", 107, 264, 141, -1, -1, 2, 0, 5, 14)}, - {ElementType::Hs, ElementData("Hs", 108, 277, 134, -1, -1, 2, 0, 6, 14)}, - {ElementType::Mt, ElementData("Mt", 109, 268, 129, -1, -1, 2, 0, 7, 14)}, - {ElementType::Ds, ElementData("Ds", 110, 281, 128, -1, -1, 2, 0, 8, 14)}, - {ElementType::Rg, ElementData("Rg", 111, 280, 121, -1, -1, 2, 0, 9, 14)}, - {ElementType::Cn, ElementData("Cn", 112, 285, 122, -1, -1, 2, 0, 10, 14)}, + {ElementType::none, ElementData("None", 0, 0.0, 0.0, 0.0, 0.0, -1, 0, 0, 0, 0)}, + {ElementType::H, ElementData("H", 1, 1.0079, 32, 109, 2.2, -1, 1, 0, 0, 0)}, + {ElementType::He, ElementData("He", 2, 4.0026, 37, 140, -1, -1, 2, 0, 0, 0)}, + {ElementType::Li, ElementData("Li", 3, 6.941, 130, 182, 0.98, -1, 1, 0, 0, 0)}, + {ElementType::Be, ElementData("Be", 4, 9.0122, 99, 153, 1.57, -1, 2, 0, 0, 0)}, + {ElementType::B, ElementData("B", 5, 10.811, 84, 192, 2.04, -1, 2, 1, 0, 0)}, + {ElementType::C, ElementData("C", 6, 12.011, 75, 170, 2.55, -1, 2, 2, 0, 0)}, + {ElementType::N, ElementData("N", 7, 14.007, 71, 155, 3.04, -1, 2, 3, 0, 0)}, + {ElementType::O, ElementData("O", 8, 15.999, 64, 152, 3.44, -1, 2, 4, 0, 0)}, + {ElementType::F, ElementData("F", 9, 18.988, 60, 147, 3.98, -1, 2, 5, 0, 0)}, + {ElementType::Ne, ElementData("Ne", 10, 20.180, 62, 154, -1, -1, 2, 6, 0, 0)}, + {ElementType::Na, ElementData("Na", 11, 22.990, 160, 227, 0.93, -1, 1, 0, 0, 0)}, + {ElementType::Mg, ElementData("Mg", 12, 24.305, 140, 173, 1.31, -1, 2, 0, 0, 0)}, + {ElementType::Al, ElementData("Al", 13, 26.982, 124, 184, 1.61, -1, 2, 1, 0, 0)}, + {ElementType::Si, ElementData("Si", 14, 28.086, 114, 210, 1.90, -1, 2, 2, 0, 0)}, + {ElementType::P, ElementData("P", 15, 30.974, 109, 180, 2.19, -1, 2, 3, 0, 0)}, + {ElementType::S, ElementData("S", 16, 32.065, 104, 180, 2.58, -1, 2, 4, 0, 0)}, + {ElementType::Cl, ElementData("Cl", 17, 35.453, 100, 175, 3.16, -1, 2, 5, 0, 0)}, + {ElementType::Ar, ElementData("Ar", 18, 39.948, 101, 188, -1, -1, 2, 6, 0, 0)}, + {ElementType::K, ElementData("K", 19, 39.098, 200, 275, 0.82, -1, 1, 0, 0, 0)}, + {ElementType::Ca, ElementData("Ca", 20, 40.078, 174, 231, 1.0, -1, 2, 0, 0, 0)}, + {ElementType::Sc, ElementData("Sc", 21, 44.956, 159, 215, 1.36, -1, 2, 0, 1, 0)}, + {ElementType::Ti, ElementData("Ti", 22, 47.867, 148, 211, 1.54, -1, 2, 0, 2, 0)}, + {ElementType::V, ElementData("V", 23, 50.942, 144, 207, 1.63, -1, 2, 0, 3, 0)}, + {ElementType::Cr, ElementData("Cr", 24, 51.996, 130, 206, 1.66, -1, 2, 0, 4, 0)}, + {ElementType::Mn, ElementData("Mn", 25, 54.938, 129, 205, 1.55, -1, 2, 0, 5, 0)}, + {ElementType::Fe, ElementData("Fe", 26, 55.938, 124, 204, 1.83, -1, 2, 0, 6, 0)}, + {ElementType::Co, ElementData("Co", 27, 58.933, 118, 200, 1.88, -1, 2, 0, 7, 0)}, + {ElementType::Ni, ElementData("Ni", 28, 58.693, 117, 197, 1.91, -1, 2, 0, 8, 0)}, + {ElementType::Cu, ElementData("Cu", 29, 63.546, 122, 196, 1.90, -1, 2, 0, 9, 0)}, + {ElementType::Zn, ElementData("Zn", 30, 65.38, 120, 201, 1.65, -1, 2, 0, 10, 0)}, + {ElementType::Ga, ElementData("Ga", 31, 69.723, 123, 187, 1.81, -1, 2, 1, 10, 0)}, + {ElementType::Ge, ElementData("Ge", 32, 72.64, 120, 211, 2.01, -1, 2, 2, 10, 0)}, + {ElementType::As, ElementData("As", 33, 74.922, 120, 185, 2.18, -1, 2, 3, 10, 0)}, + {ElementType::Se, ElementData("Se", 34, 78.96, 118, 190, 2.55, -1, 2, 4, 10, 0)}, + {ElementType::Br, ElementData("Br", 35, 79.904, 117, 185, 2.96, -1, 2, 5, 10, 0)}, + {ElementType::Kr, ElementData("Kr", 36, 83.798, 116, 202, -1.0, -1, 2, 6, 10, 0)}, + {ElementType::Rb, ElementData("Rb", 37, 83.468, 215, 303, 0.82, -1, 1, 0, 0, 0)}, + {ElementType::Sr, ElementData("Sr", 38, 87.62, 190, 249, 0.95, -1, 2, 0, 0, 0)}, + {ElementType::Y, ElementData("Y", 39, 88.906, 176, 232, 1.22, -1, 2, 0, 1, 0)}, + {ElementType::Zr, ElementData("Zr", 40, 91.224, 164, 223, 1.33, -1, 2, 0, 2, 0)}, + {ElementType::Nb, ElementData("Nb", 41, 92.906, 156, 218, 1.6, -1, 2, 0, 3, 0)}, + {ElementType::Mo, ElementData("Mo", 42, 95.96, 146, 217, 2.16, -1, 2, 0, 4, 0)}, + {ElementType::Tc, ElementData("Tc", 43, 98.91, 138, 216, 2.10, -1, 2, 0, 5, 0)}, + {ElementType::Ru, ElementData("Ru", 44, 101.07, 136, 213, 2.2, -1, 2, 0, 6, 0)}, + {ElementType::Rh, ElementData("Rh", 45, 102.91, 134, 210, 2.28, -1, 2, 0, 7, 0)}, + {ElementType::Pd, ElementData("Pd", 46, 106.42, 130, 210, 2.2, -1, 2, 0, 8, 0)}, + {ElementType::Ag, ElementData("Ag", 47, 107.87, 136, 211, 1.93, -1, 2, 0, 9, 0)}, + {ElementType::Cd, ElementData("Cd", 48, 112.41, 140, 218, 1.69, -1, 2, 0, 10, 0)}, + {ElementType::In, ElementData("In", 49, 114.82, 142, 193, 1.78, -1, 2, 1, 10, 0)}, + {ElementType::Sn, ElementData("Sn", 50, 118.71, 140, 217, 1.96, -1, 2, 2, 10, 0)}, + {ElementType::Sb, ElementData("Sb", 51, 121.76, 140, 206, 2.05, -1, 2, 3, 10, 0)}, + {ElementType::Te, ElementData("Te", 52, 127.60, 137, 206, 2.1, -1, 2, 4, 10, 0)}, + {ElementType::I, ElementData("I", 53, 126.90, 136, 198, 2.66, -1, 2, 5, 10, 0)}, + {ElementType::Xe, ElementData("Xe", 54, 131.29, 136, 216, 2.6, -1, 2, 6, 10, 0)}, + {ElementType::Cs, ElementData("Cs", 55, 132.91, 238, 343, 0.79, -1, 1, 0, 0, 0)}, + {ElementType::Ba, ElementData("Ba", 56, 137.33, 206, 268, 0.89, -1, 2, 0, 0, 0)}, + {ElementType::La, ElementData("La", 57, 138.91, 194, 243, 1.1, -1, 2, 0, 0, 1)}, + {ElementType::Ce, ElementData("Ce", 58, 140.12, 184, 242, 1.12, -1, 2, 0, 0, 2)}, + {ElementType::Pr, ElementData("Pr", 59, 140.91, 190, 240, 1.13, -1, 2, 0, 0, 3)}, + {ElementType::Nd, ElementData("Nd", 60, 144.24, 188, 239, 1.14, -1, 2, 0, 0, 4)}, + {ElementType::Pm, ElementData("Pm", 61, 146.90, 186, 238, -1.0, -1, 2, 0, 0, 5)}, + {ElementType::Sm, ElementData("Sm", 62, 150.36, 185, 236, 1.17, -1, 2, 0, 0, 6)}, + {ElementType::Eu, ElementData("Eu", 63, 151.96, 183, 235, -1.0, -1, 2, 0, 0, 7)}, + {ElementType::Gd, ElementData("Gd", 64, 157.25, 182, 234, 1.2, -1, 2, 0, 0, 8)}, + {ElementType::Tb, ElementData("Tb", 65, 158.93, 181, 233, -1.0, -1, 2, 0, 0, 9)}, + {ElementType::Dy, ElementData("Dy", 66, 162.50, 180, 231, 1.22, -1, 2, 0, 0, 10)}, + {ElementType::Ho, ElementData("Ho", 67, 164.93, 179, 230, 1.23, -1, 2, 0, 0, 11)}, + {ElementType::Er, ElementData("Er", 68, 167.26, 177, 229, 1.24, -1, 2, 0, 0, 12)}, + {ElementType::Tm, ElementData("Tm", 69, 168.93, 177, 227, 1.25, -1, 2, 0, 0, 13)}, + {ElementType::Yb, ElementData("Yb", 70, 173.05, 178, 226, -1.0, -1, 2, 0, 0, 14)}, + {ElementType::Lu, ElementData("Lu", 71, 174.97, 174, 224, 1.0, -1, 2, 0, 1, 14)}, + {ElementType::Hf, ElementData("Hf", 72, 178.49, 164, 223, 1.3, -1, 2, 0, 2, 14)}, + {ElementType::Ta, ElementData("Ta", 73, 180.95, 158, 222, 1.5, -1, 2, 0, 3, 14)}, + {ElementType::W, ElementData("W", 74, 183.84, 150, 218, 1.7, -1, 2, 0, 4, 14)}, + {ElementType::Re, ElementData("Re", 75, 186.21, 141, 216, 1.9, -1, 2, 0, 5, 14)}, + {ElementType::Os, ElementData("Os", 76, 190.23, 136, 216, 2.2, -1, 2, 0, 6, 14)}, + {ElementType::Ir, ElementData("Ir", 77, 192.22, 132, 213, 2.2, -1, 2, 0, 7, 14)}, + {ElementType::Pt, ElementData("Pt", 78, 195.08, 130, 213, 2.2, -1, 2, 0, 8, 14)}, + {ElementType::Au, ElementData("Au", 79, 196.97, 130, 214, 2.4, -1, 2, 0, 9, 14)}, + {ElementType::Hg, ElementData("Hg", 80, 200.59, 132, 223, 1.9, -1, 2, 0, 10, 14)}, + {ElementType::Tl, ElementData("Tl", 81, 204.38, 144, 196, 1.8, -1, 2, 1, 10, 14)}, + {ElementType::Pb, ElementData("Pb", 82, 207.2, 145, 202, 1.8, -1, 2, 2, 10, 14)}, + {ElementType::Bi, ElementData("Bi", 83, 208.98, 150, 207, 1.9, -1, 2, 3, 10, 14)}, + {ElementType::Po, ElementData("Po", 84, 209.98, 142, 197, 2.0, -1, 2, 4, 10, 14)}, + {ElementType::At, ElementData("At", 85, 210, 148, 202, 2.2, -1, 2, 5, 10, 14)}, + {ElementType::Rn, ElementData("Rn", 86, 222, 146, 220, -1.0, -1, 2, 6, 10, 14)}, + {ElementType::Fr, ElementData("Fr", 87, 223, 242, 348, 0.7, -1, 1, 0, 0, 0)}, + {ElementType::Ra, ElementData("Ra", 88, 226.03, 211, 283, 0.9, -1, 2, 0, 0, 0)}, + {ElementType::Ac, ElementData("Ac", 89, 227, 201, 247, 1.1, -1, 2, 0, 0, 1)}, + {ElementType::Th, ElementData("Th", 90, 232.04, 190, 245, 1.3, -1, 2, 0, 0, 2)}, + {ElementType::Pa, ElementData("Pa", 91, 231.04, 184, 243, 1.5, -1, 2, 0, 0, 3)}, + {ElementType::U, ElementData("U", 92, 238.03, 183, 241, 1.7, -1, 2, 0, 0, 4)}, + {ElementType::Np, ElementData("Np", 93, 237.05, 180, 239, 1.3, -1, 2, 0, 0, 5)}, + {ElementType::Pu, ElementData("Pu", 94, 244.10, 180, 243, 1.3, -1, 2, 0, 0, 6)}, + {ElementType::Am, ElementData("Am", 95, 243.10, 173, 244, -1.0, -1, 2, 0, 0, 7)}, + {ElementType::Cm, ElementData("Cm", 96, 247.10, 168, 245, -1.0, -1, 2, 0, 0, 8)}, + {ElementType::Bk, ElementData("Bk", 97, 247.10, 168, 244, -1.0, -1, 2, 0, 0, 9)}, + {ElementType::Cf, ElementData("Cf", 98, 251.10, 168, 245, -1.0, -1, 2, 0, 0, 10)}, + {ElementType::Es, ElementData("Es", 99, 254.10, 165, 245, -1.0, -1, 2, 0, 0, 11)}, + {ElementType::Fm, ElementData("Fm", 100, 257.10, 167, 245, -1.0, -1, 2, 0, 0, 12)}, + {ElementType::Md, ElementData("Md", 101, 258, 173, 246, -1.0, -1, 2, 0, 0, 13)}, + {ElementType::No, ElementData("No", 102, 259, 176, 246, -1.0, -1, 2, 0, 0, 14)}, + {ElementType::Lr, ElementData("Lr", 103, 262, 161, 246, -1.0, -1, 2, 0, 1, 14)}, + {ElementType::Rf, ElementData("Rf", 104, 261, 157, -1, -1.0, -1, 2, 0, 2, 14)}, + {ElementType::Db, ElementData("Db", 105, 262, 149, -1, -1.0, -1, 2, 0, 3, 14)}, + {ElementType::Sg, ElementData("Sg", 106, 266, 143, -1, -1.0, -1, 2, 0, 4, 14)}, + {ElementType::Bh, ElementData("Bh", 107, 264, 141, -1, -1.0, -1, 2, 0, 5, 14)}, + {ElementType::Hs, ElementData("Hs", 108, 277, 134, -1, -1.0, -1, 2, 0, 6, 14)}, + {ElementType::Mt, ElementData("Mt", 109, 268, 129, -1, -1.0, -1, 2, 0, 7, 14)}, + {ElementType::Ds, ElementData("Ds", 110, 281, 128, -1, -1.0, -1, 2, 0, 8, 14)}, + {ElementType::Rg, ElementData("Rg", 111, 280, 121, -1, -1.0, -1, 2, 0, 9, 14)}, + {ElementType::Cn, ElementData("Cn", 112, 285, 122, -1, -1.0, -1, 2, 0, 10, 14)}, }; // clang-format on return map; diff --git a/src/Utils/Utils/Geometry/ElementData.h b/src/Utils/Utils/Geometry/ElementData.h index 9680b6e..2b04892 100644 --- a/src/Utils/Utils/Geometry/ElementData.h +++ b/src/Utils/Utils/Geometry/ElementData.h @@ -59,6 +59,7 @@ class ElementDataSingleton { * @param mass The mass (isotope average, precision: 3 digits). * @param covalentRadiusInPicometers The covalent radius in pm. * @param vdWRadiusInPicometers The van der Waals radius in pm. + * @param paulingElectronegativity The Pauling electronegativity. * @param valElectrons The number of electrons in the valence shell. * @param sElectrons The number of s-electrons in the valence shell. * @param pElectrons The number of p-electrons in the valence shell. @@ -66,13 +67,14 @@ class ElementDataSingleton { * @param fElectrons The number of f-electrons in the valence shell. */ ElementData(std::string symbol, int Z, double mass, double covalentRadiusInPicometers = -1, - double vdWRadiusInPicometers = -1, int valElectrons = -1, int sElectrons = -1, int pElectrons = -1, - int dElectrons = -1, int fElectrons = -1) + double vdWRadiusInPicometers = -1, double paulingElectronegativity = -1, int valElectrons = -1, + int sElectrons = -1, int pElectrons = -1, int dElectrons = -1, int fElectrons = -1) : d_symbol(std::move(symbol)), d_Z(Z), d_mass(mass), d_covalentRadius(covalentRadiusInPicometers / 100 * bohr_per_angstrom), d_vdWRadius(vdWRadiusInPicometers / 100 * bohr_per_angstrom), + d_paulingElectronegativity(paulingElectronegativity), d_valElectrons(valElectrons), d_sElectrons(sElectrons), d_pElectrons(pElectrons), @@ -187,6 +189,21 @@ class ElementDataSingleton { throw DataNotAvailable(); } + /** + * @brief Getter for the Pauling electronegativity. + * + * CRC Handbook of Physics and Chemistry 85th Edition, + * CRC Press/Taylor and Francis, Boca Raton, FL, Molecular Structure and Spectroscopy, + * Sec. 9, page 74.\n + * Original references: + * Pauling, L., The Nature of the Chemical Bond, Third Edition, Cornell University Press, Ithaca, NY, 1960.\n + * Allen, L. C., J. Am. Chem. Soc., 111, 9003, 1989.\n + * Allred, A. L., J. Inorg. Nucl. Chem., 17, 215, 1961. + * + * @return double Returns the Pauling electronegativity. + */ + double paulingElectronegativity() const; + private: std::string d_symbol{}; @@ -194,6 +211,7 @@ class ElementDataSingleton { double d_mass{-1}; double d_covalentRadius{-1}; double d_vdWRadius{-1}; + double d_paulingElectronegativity{-1}; int d_valElectrons{-1}; int d_sElectrons{-1}; diff --git a/src/Utils/Utils/Geometry/ElementInfo.cpp b/src/Utils/Utils/Geometry/ElementInfo.cpp index 630b9a4..2fde6ed 100644 --- a/src/Utils/Utils/Geometry/ElementInfo.cpp +++ b/src/Utils/Utils/Geometry/ElementInfo.cpp @@ -209,6 +209,14 @@ ElementType ElementInfo::base(ElementType isotope) { } } +double ElementInfo::paulingElectronegativity(ElementType e) { + if (A(e) != 0) { + e = base(e); + } + + return Constants::ElementDataSingleton::lookup(e).paulingElectronegativity(); +} + double ElementInfo::vdwRadius(ElementType e) { if (A(e) != 0) { e = base(e); diff --git a/src/Utils/Utils/Geometry/ElementInfo.h b/src/Utils/Utils/Geometry/ElementInfo.h index e677c9c..51bbb2c 100644 --- a/src/Utils/Utils/Geometry/ElementInfo.h +++ b/src/Utils/Utils/Geometry/ElementInfo.h @@ -104,6 +104,22 @@ class ElementInfo { */ static double vdwRadius(ElementType e); + /** + * @brief Getter for the Pauling electronegativity. + * @param e The ElementType. + * + * CRC Handbook of Physics and Chemistry 85th Edition, + * CRC Press/Taylor and Francis, Boca Raton, FL, Molecular Structure and Spectroscopy, + * Sec. 9, page 74.\n + * Original references: + * Pauling, L., The Nature of the Chemical Bond, Third Edition, Cornell University Press, Ithaca, NY, 1960.\n + * Allen, L. C., J. Am. Chem. Soc., 111, 9003, 1989.\n + * Allred, A. L., J. Inorg. Nucl. Chem., 17, 215, 1961. + * + * @return double Returns the Pauling electronegativity. + */ + static double paulingElectronegativity(ElementType e); + /** * @brief Getter for the atomic number * @param e The ElementType. diff --git a/src/Utils/Utils/Geometry/PeriodicSystem.cpp b/src/Utils/Utils/Geometry/PeriodicSystem.cpp index 9785803..5a460c1 100644 --- a/src/Utils/Utils/Geometry/PeriodicSystem.cpp +++ b/src/Utils/Utils/Geometry/PeriodicSystem.cpp @@ -12,6 +12,10 @@ namespace Scine { namespace Utils { +PeriodicSystem::PeriodicSystem(const PeriodicSystem& other) + : PeriodicSystem(other.pbc, other.atoms, other.solidStateAtomIndices) { +} + PeriodicSystem::PeriodicSystem(const PeriodicBoundaries& pbc, int N, std::unordered_set solidStateAtomIndices) : PeriodicSystem(pbc, AtomCollection(N), std::move(solidStateAtomIndices)) { } @@ -28,6 +32,16 @@ PeriodicSystem::PeriodicSystem(const PeriodicBoundaries& pbc, AtomCollection ato [&](unsigned i) { return static_cast(i) < this->atoms.size(); })) { throw std::logic_error("At least one of the given solid state indices is not valid for the given AtomCollection."); } + canonicalize(); +} + +void PeriodicSystem::canonicalize() { + centerAndTranslateAtomsIntoCell(); + Eigen::Matrix3d rotation = pbc.getCanonicalizationRotationMatrix(); + if (rotation != Eigen::Matrix3d::Identity()) { + this->pbc.canonicalize(); + this->atoms.setPositions(this->atoms.getPositions() * rotation); + } centerAndTranslateAtomsIntoCell(); } @@ -214,8 +228,9 @@ bool PeriodicSystem::operator!=(const PeriodicSystem& other) const { } PeriodicSystem PeriodicSystem::operator*(const Eigen::Vector3i& scalingFactors) const { - PeriodicSystem super = PeriodicSystem(this->pbc, this->atoms, this->solidStateAtomIndices); - return super *= scalingFactors; + PeriodicSystem super = PeriodicSystem(*this); + super *= scalingFactors; + return super; } PeriodicSystem& PeriodicSystem::operator*=(const Eigen::Vector3i& scalingFactors) { @@ -223,6 +238,12 @@ PeriodicSystem& PeriodicSystem::operator*=(const Eigen::Vector3i& scalingFactors throw std::runtime_error( "At least one scaling factor of the periodic system is less than 1, which is not possible."); } + for (int i = 0; i < 3; ++i) { + if (scalingFactors[i] > 1 && !pbc.getPeriodicity()[i]) { + throw std::runtime_error( + "Scaling factor of the periodic system is larger than 1 for a dimension which is not periodic."); + } + } // do operations within cell -> save original state to translate back afterwards const PositionCollection originalPositions = atoms.getPositions(); centerAndTranslateAtomsIntoCell(); @@ -269,13 +290,23 @@ PeriodicSystem& PeriodicSystem::operator*=(const Eigen::Vector3i& scalingFactors } PeriodicSystem PeriodicSystem::operator*(int scalingFactor) const { - Eigen::Vector3i scaling = Eigen::Vector3i::Constant(scalingFactor); - return *this * scaling; + PeriodicSystem super = (*this); + super *= scalingFactor; + return super; } PeriodicSystem& PeriodicSystem::operator*=(int scalingFactor) { + if (scalingFactor < 1) { + throw std::runtime_error("Specified scaling factor of " + std::to_string(scalingFactor) + ", but it must be at least 1"); + } Eigen::Vector3i scaling = Eigen::Vector3i::Constant(scalingFactor); - return *this *= scaling; + for (int i = 0; i < 3; ++i) { + if (!pbc.getPeriodicity()[i]) { + scaling[i] = 1; + } + } + *this *= scaling; + return *this; } } /* namespace Utils */ diff --git a/src/Utils/Utils/Geometry/PeriodicSystem.h b/src/Utils/Utils/Geometry/PeriodicSystem.h index b2408dd..4dec0f0 100644 --- a/src/Utils/Utils/Geometry/PeriodicSystem.h +++ b/src/Utils/Utils/Geometry/PeriodicSystem.h @@ -32,6 +32,7 @@ class PeriodicSystem { std::tuple, std::unordered_map>; public: + PeriodicSystem(const PeriodicSystem& other); /** * @brief Construct a new PeriodicSystem object. * @@ -230,6 +231,8 @@ class PeriodicSystem { // member has been changed in the meantime of last image construction AtomCollection _lastImageConstructedAtoms; + void canonicalize(); + inline void clearImageAtoms() { _imageAtoms = nullptr; _imageBondOrderCollection = nullptr; diff --git a/src/Utils/Utils/Geometry/StructuralCompletion.h b/src/Utils/Utils/Geometry/StructuralCompletion.h index 0adfd36..1e16ddd 100644 --- a/src/Utils/Utils/Geometry/StructuralCompletion.h +++ b/src/Utils/Utils/Geometry/StructuralCompletion.h @@ -87,7 +87,7 @@ class StructuralCompletion { static void generate1TriangleCornerFrom2Others(const Eigen::Vector3d& v1, const Eigen::Vector3d& v2, Eigen::Ref v3); - /// @brief The standard tetrahedral angel in rad. + /// @brief The standard tetrahedral angle in rad. static constexpr double tetrahedronAngle = 109.4712 * Constants::rad_per_degree; }; diff --git a/src/Utils/Utils/GeometryOptimization/GeometryOptimization.h b/src/Utils/Utils/GeometryOptimization/GeometryOptimization.h index d5a95e1..2b2e060 100644 --- a/src/Utils/Utils/GeometryOptimization/GeometryOptimization.h +++ b/src/Utils/Utils/GeometryOptimization/GeometryOptimization.h @@ -44,12 +44,15 @@ template bool settingsMakeSense(const OptimizerType& optimizer) { const auto settings = optimizer.getCalculatorSettings(); if (!settings || !settings->valueExists(SettingsNames::selfConsistenceCriterion)) { - throw std::runtime_error("Necessary '" + std::string(SettingsNames::selfConsistenceCriterion) + - "' not implemented in the applied calculator."); + std::cout << "Warning: Setting '" + std::string(SettingsNames::selfConsistenceCriterion) + "' not implemented in the applied calculator." + << std::endl; + return true; + } + else { + const double energyThreshold = settings->getDouble(SettingsNames::selfConsistenceCriterion); + const GradientBasedCheck check = optimizer.getConvergenceCheck(); + return energyThresholdIsLowEnough(energyThreshold, check); } - const double energyThreshold = settings->getDouble(SettingsNames::selfConsistenceCriterion); - const GradientBasedCheck check = optimizer.getConvergenceCheck(); - return energyThresholdIsLowEnough(energyThreshold, check); } // @brief The Convergence of the NtOptimizer is independent of the energy convergence threshold diff --git a/src/Utils/Utils/GeometryOptimization/NtOptimizer.cpp b/src/Utils/Utils/GeometryOptimization/NtOptimizer.cpp index b93a1dd..685ddb2 100644 --- a/src/Utils/Utils/GeometryOptimization/NtOptimizer.cpp +++ b/src/Utils/Utils/GeometryOptimization/NtOptimizer.cpp @@ -184,6 +184,14 @@ void NtOptimizer::sanityCheck(const AtomCollection& atoms) const { } } } + if (std::find(possibleExtractionOptions.begin(), possibleExtractionOptions.end(), extractionCriterion) == + possibleExtractionOptions.end()) { + std::string options = "\nValid options are:\n"; + for (const auto& criterion : possibleExtractionOptions) { + options += "'" + criterion + "'\n"; + } + throw std::logic_error("Extract criterion " + extractionCriterion + " is not a valid option." + options); + } } void NtOptimizer::updateGradients(const AtomCollection& atoms, const double& /* energy */, @@ -337,7 +345,11 @@ PositionCollection NtOptimizer::extractTsGuess() const { if (maximaList.empty()) { throw std::runtime_error("No transition state guess was found in Newton Trajectory scan."); } - // Extract TS guess from point with highest energy + // Extract TS guess according to given criterion + if (extractionCriterion == ntExtractFirst) { + return _trajectory[maximaList.back()]; // back because maximalist starts from back + } + // Extract TS guess from point with the highest energy double maxValue = -std::numeric_limits::max(); int maxIndex = -1; for (auto& i : maximaList) { diff --git a/src/Utils/Utils/GeometryOptimization/NtOptimizer.h b/src/Utils/Utils/GeometryOptimization/NtOptimizer.h index 30c3f01..dbf373f 100644 --- a/src/Utils/Utils/GeometryOptimization/NtOptimizer.h +++ b/src/Utils/Utils/GeometryOptimization/NtOptimizer.h @@ -57,9 +57,12 @@ class NtOptimizer : public Optimizer { static constexpr const char* ntFixedNumberOfMicroCycles = "nt_fixed_number_of_micro_cycles"; static constexpr const char* ntNumberOfMicroCycles = "nt_number_of_micro_cycles"; static constexpr const char* ntFilterPasses = "nt_filter_passes"; + static constexpr const char* ntExtractionCriterion = "nt_extraction_criterion"; static constexpr const char* ntCoordinateSystemKey = "nt_coordinate_system"; static constexpr const char* ntFixedAtomsKey = "nt_constrained_atoms"; static constexpr const char* ntMovableSide = "nt_movable_side"; + static constexpr const char* ntExtractHighest = "highest_maximum"; + static constexpr const char* ntExtractFirst = "first_maximum"; /** * @brief Construct a new NtOptimizer object. * @param calculator The calculator to be used for the underlying single point/gradient calculations. @@ -137,6 +140,15 @@ class NtOptimizer : public Optimizer { * the reaction curve. */ int filterPasses = 10; + // @brief possible options for extraction + const std::vector possibleExtractionOptions = { + ntExtractHighest, + ntExtractFirst, + }; + /** + * @brief Criterion to extract a TS guess from the trajectory + */ + std::string extractionCriterion = possibleExtractionOptions.front(); /// @brief The special convergence settings for this optimizer. struct NtConvergenceStub { /// @brief The maximum number of iterations. diff --git a/src/Utils/Utils/GeometryOptimization/NtOptimizer2.cpp b/src/Utils/Utils/GeometryOptimization/NtOptimizer2.cpp index 56c41bd..4bc4e7e 100644 --- a/src/Utils/Utils/GeometryOptimization/NtOptimizer2.cpp +++ b/src/Utils/Utils/GeometryOptimization/NtOptimizer2.cpp @@ -53,7 +53,7 @@ int NtOptimizer2::optimize(AtomCollection& atoms, Core::Log& log) { // Apply Cartesian constraints auto bos = results.get(); auto gradientMatrix = results.get(); - this->updateGradients(atoms, value, gradientMatrix, bos); + this->updateGradients(atoms, value, gradientMatrix, bos, cycle); gradients = Eigen::Map(gradientMatrix.data(), nAtoms * 3); }; // Run micro iterations @@ -100,7 +100,7 @@ int NtOptimizer2::optimize(AtomCollection& atoms, Core::Log& log) { auto bos = results.get(); this->triggerObservers(cycle, value, Eigen::Map(coordinates.data(), nAtoms * 3)); // Evaluate additional force - this->updateGradients(atoms, value, gradients, bos, true); + this->updateGradients(atoms, value, gradients, bos, cycle, true); // Add energy to energy vector _values.push_back(value); // Add coordinates to trajectory @@ -176,15 +176,105 @@ void NtOptimizer2::sanityCheck(const AtomCollection& atoms) const { } } } + if (std::find(possibleExtractionOptions.begin(), possibleExtractionOptions.end(), extractionCriterion) == + possibleExtractionOptions.end()) { + std::string options = "\nValid options are:\n"; + for (const auto& criterion : possibleExtractionOptions) { + options += "'" + criterion + "'\n"; + } + throw std::logic_error("Extract criterion " + extractionCriterion + " is not a valid option." + options); + } } void NtOptimizer2::updateGradients(const AtomCollection& atoms, const double& /* energy */, GradientCollection& gradients, - const BondOrderCollection& bondOrders, bool addForce) const { - // define reaction coordinate based on LHS and RHS list + const BondOrderCollection& bondOrders, int cycle, bool addForce) { + const auto& positions = atoms.getPositions(); + eliminateReactiveAtomsGradients(positions, gradients); + auto reactions = inferReactions(bondOrders); + ReactionMapping associations = reactions.first; + ReactionMapping dissociations = reactions.second; + // define reaction coordinate GradientCollection reactionCoordinate(gradients); reactionCoordinate.setZero(); + double maxScale = 0.0; + // association + for (const auto& [left, right] : associations) { + double r12cov = NtUtils::smallestCovalentRadius(atoms, left) + NtUtils::smallestCovalentRadius(atoms, right); + Displacement c2c = NtUtils::centerToCenterVector(positions, left, right); + double dist = c2c.norm(); + double bo = 0.0; + for (const auto& l : left) { + for (const auto& r : right) { + bo += bondOrders.getOrder(l, r); + } + } + /* Check if one should keep pushing along the reaction coordinate */ + /* Keep pushing until 10% above attractive bond order stop OR distance 10% below distance stop */ + if (bo > 1.1 * check.attractiveBondOrderStop || dist < 0.9 * this->check.attractiveDistanceStop * r12cov) { + if (_firstCoordinateReachedIndex == -1) { + _firstCoordinateReachedIndex = cycle; + } + continue; + } + double scale = dist - r12cov; + if (scale < 0) { + scale = 0.1; + } + maxScale = std::max(maxScale, scale); + c2c *= scale / dist; + for (const auto& l : left) { + for (const auto& r : right) { + reactionCoordinate.row(l) += c2c; + reactionCoordinate.row(r) -= c2c; + } + } + } + // dissociation + for (const auto& [left, right] : dissociations) { + double bo = 0.0; + Displacement c2c = NtUtils::centerToCenterVector(positions, left, right); + double dist = c2c.norm(); + for (const auto& l : left) { + for (const auto& r : right) { + bo += bondOrders.getOrder(l, r); + } + } + /* Check if one should keep pulling along the reaction coordinate */ + /* Keep pulling until 30% bellow repulsive bond order stop */ + if (bo < 0.7 * this->check.repulsiveBondOrderStop) { + if (_firstCoordinateReachedIndex == -1) { + _firstCoordinateReachedIndex = cycle; + } + continue; + } + const double scale = 0.8; + maxScale = std::max(maxScale, scale); + c2c *= scale / dist; + for (const auto& l : left) { + for (const auto& r : right) { + reactionCoordinate.row(l) -= c2c; + reactionCoordinate.row(r) += c2c; + } + } + } + if (maxScale != 0.0) { + reactionCoordinate *= 0.5 * (totalForceNorm / maxScale); + } - const auto& positions = atoms.getPositions(); + /* Update gradients */ + // replace gradient along reaction coordinate with fixed force for both sides + if (addForce) { + gradients += reactionCoordinate; + } + // Apply Cartesian constraints + if (!this->fixedAtoms.empty()) { + for (const auto& a : fixedAtoms) { + gradients.row(a).setZero(); + } + } +} + +void NtOptimizer2::eliminateReactiveAtomsGradients(const PositionCollection& positions, GradientCollection& gradients) const { // remove gradient along reaction coordinates /* Loop over all reactive atoms */ for (auto const& reactiveAtomIndex : _reactiveAtomsList) { @@ -206,16 +296,14 @@ void NtOptimizer2::updateGradients(const AtomCollection& atoms, const double& /* PositionCollection constraintsMatrixOrtho = (constraintsMatrix.transpose() * saesEigenVec).transpose(); /* Determine rank of matrix (eigenvalues > 1e-6) */ unsigned int rank = 0; - int lastRowIndex = constraintsMatrixOrtho.rows() - 1; + long lastRowIndex = constraintsMatrixOrtho.rows() - 1; /* Determine rank of matrix */ - for (int i = constraintsMatrixOrtho.rows() - 1; i > -1; i--) { + for (long i = constraintsMatrixOrtho.rows() - 1; i > -1; i--) { /* Break loop over eigenvalues if it is < 1e-6 */ if (saesEigenVal.col(0)[i] < 1e-6) { break; } - else { - rank += 1; - } + rank += 1; } /* Manipulate gradients of reactive atom */ /* Can not be orthogonal to more than 2 axis, hence set the gradient to 0 */ @@ -230,94 +318,130 @@ void NtOptimizer2::updateGradients(const AtomCollection& atoms, const double& /* proj /= proj.norm(); gradients.row(reactiveAtomIndex) = (gradients.row(reactiveAtomIndex).dot(proj)) * proj; } - /* Substract projection along the axis from gradient */ + /* Subtract projection along the axis from gradient */ else if (rank == 1) { Eigen::Vector3d a2a = constraintsMatrixOrtho.row(lastRowIndex) / (constraintsMatrixOrtho.row(lastRowIndex).norm()); gradients.row(reactiveAtomIndex) -= (gradients.row(reactiveAtomIndex).dot(a2a)) * a2a; } } +} - double maxScale = 0.0; - for (unsigned int i = 0; i < associationList.size() / 2; i++) { - const auto& l = associationList[2 * i]; - const auto& r = associationList[2 * i + 1]; - const auto bo = bondOrders.getOrder(l, r); - double r12cov = ElementInfo::covalentRadius(atoms.getElement(l)) + ElementInfo::covalentRadius(atoms.getElement(r)); - double dist = (positions.row(l) - positions.row(r)).norm(); - /* Check if one should keep pushing along the reaction coordinate */ - /* Keep pushing until 10% above attractive bond order stop OR distance 10% below distance stop */ - if (bo > 1.1 * check.attractiveBondOrderStop || dist < 0.9 * this->check.attractiveDistanceStop * r12cov) { - continue; - } - Eigen::Vector3d a2a = positions.row(l) - positions.row(r); - const double norm = a2a.norm(); - r12cov = ElementInfo::covalentRadius(atoms.getElement(l)) + ElementInfo::covalentRadius(atoms.getElement(r)); - double scale = norm - r12cov; - if (scale < 0) { - scale = 0.1; +std::pair +NtOptimizer2::inferReactions(const BondOrderCollection& bondOrders) const { + std::vector maps; + // we do the exact same thing for associations and dissociations + std::vector> reactionLists = {associationList, dissociationList}; + for (const auto& reactionList : reactionLists) { + // we restructure the single list with included logic into list of pairs + std::vector> reactions; + for (unsigned int i = 0; i < reactionList.size() / 2; i++) { + const auto& left = reactionList[2 * i]; + const auto& right = reactionList[2 * i + 1]; + reactions.emplace_back(left, right); } - maxScale = std::max(maxScale, scale); - a2a *= scale / norm; - reactionCoordinate.row(l) += a2a; - reactionCoordinate.row(r) -= a2a; - } - for (unsigned int i = 0; i < dissociationList.size() / 2; i++) { - const auto& l = dissociationList[2 * i]; - const auto& r = dissociationList[2 * i + 1]; - const auto bo = bondOrders.getOrder(l, r); - /* Check if one should keep pulling along the reaction coordinate */ - /* Keep pulling until 30% bellow repulsive bond order stop */ - if (bo < 0.7 * this->check.repulsiveBondOrderStop) { - continue; + // now we need to combine duplicates + // we pretty much build now a map that allows duplicate keys + // reason is, if we have a duplicate index X on the left side, we want to distinguish if X takes part in genuinely + // different reaction coordinates, e.g. X - Y and X - Z, or if we have an eta bond reaction coordinate X - Y-Z + // we distinguish this in our data structure with: + // + // + // for different reaction coordinates and + // + // for an eta bond reaction coordinate. The value can have any length and is not limited to two + // we determine the case based on the fact if Y and Z are part of the same molecule (second case) or not (first case) + std::vector>> rightCorrectedReactions; + std::vector sanitizedLeftIndices; + for (unsigned long i = 0; i < reactions.size(); ++i) { + if (std::find(sanitizedLeftIndices.begin(), sanitizedLeftIndices.end(), reactions[i].first) != + sanitizedLeftIndices.end()) { + continue; // no need to check same index again + } + sanitizedLeftIndices.push_back(reactions[i].first); + std::vector values = {reactions[i].second}; + for (unsigned long j = i + 1; j < reactions.size(); ++j) { + if (reactions[i].first == reactions[j].first) { + values.push_back(reactions[j].second); + } + } + if (values.size() > 1) { + // we had a duplicate, let's check if the shared values are part of the same molecule + std::vector> molecules = NtUtils::connectedNuclei(values, bondOrders); + // add entry for each molecule, since each represents a genuine reaction coordinate + for (const auto& molecule : molecules) { + rightCorrectedReactions.emplace_back(reactions[i].first, molecule); + } + } + else { + // no duplicate simply copy into new data structure + rightCorrectedReactions.emplace_back(reactions[i].first, values); + } } - Eigen::Vector3d a2a = positions.row(l) - positions.row(r); - const double norm = a2a.norm(); - const double scale = 0.8; - maxScale = std::max(maxScale, scale); - a2a *= scale / norm; - reactionCoordinate.row(l) -= a2a; - reactionCoordinate.row(r) += a2a; - } - if (maxScale != 0.0) { - reactionCoordinate *= 0.5 * (totalForceNorm / maxScale); - } - - /* Update gradients */ - // replace gradient along reaction coordinate with fixed force for both sides - if (addForce) { - gradients += reactionCoordinate; - } - // Apply Cartesian constraints - if (!this->fixedAtoms.empty()) { - for (const auto& a : fixedAtoms) { - gradients.row(a).setZero(); + // now we need to do the exact same thing for the right side and check if there are indices on the right side + // or lists of indices that have different reaction partners and check if they are again different reaction + // coordinates or eta bond formations + ReactionMapping allCorrectedReactions; + std::vector> sanitizedRightIndices; + for (unsigned long i = 0; i < rightCorrectedReactions.size(); ++i) { + if (std::find(sanitizedRightIndices.begin(), sanitizedRightIndices.end(), rightCorrectedReactions[i].second) != + sanitizedRightIndices.end()) { + continue; + } + sanitizedRightIndices.push_back(rightCorrectedReactions[i].second); + std::vector values = {rightCorrectedReactions[i].first}; + for (unsigned long j = i + 1; j < rightCorrectedReactions.size(); ++j) { + if (rightCorrectedReactions[i].second == rightCorrectedReactions[j].second) { + values.push_back(rightCorrectedReactions[j].first); + } + } + if (values.size() > 1) { + std::vector> molecules = NtUtils::connectedNuclei(values, bondOrders); + for (const auto& molecule : molecules) { + allCorrectedReactions.emplace_back(molecule, rightCorrectedReactions[i].second); + } + } + else { + allCorrectedReactions.emplace_back(values, rightCorrectedReactions[i].second); + } } + maps.push_back(allCorrectedReactions); } + return std::make_pair(maps[0], maps[1]); } bool NtOptimizer2::convergedOptimization(const AtomCollection& atoms, const BondOrderCollection& bondOrders) const { + auto reactions = inferReactions(bondOrders); + ReactionMapping associations = reactions.first; + ReactionMapping dissociations = reactions.second; const auto& positions = atoms.getPositions(); - bool done = true; - for (unsigned int i = 0; i < this->associationList.size() / 2; i++) { - const auto& l = this->associationList[2 * i]; - const auto& r = this->associationList[2 * i + 1]; - const auto bo = bondOrders.getOrder(l, r); - const double r12cov = - ElementInfo::covalentRadius(atoms.getElement(l)) + ElementInfo::covalentRadius(atoms.getElement(r)); - double dist = (positions.row(l) - positions.row(r)).norm(); + // associations + for (const auto& [left, right] : associations) { + double r12cov = NtUtils::smallestCovalentRadius(atoms, left) + NtUtils::smallestCovalentRadius(atoms, right); + Displacement c2c = NtUtils::centerToCenterVector(positions, left, right); + double dist = c2c.norm(); + double bo = 0.0; + for (const auto& l : left) { + for (const auto& r : right) { + bo += bondOrders.getOrder(l, r); + } + } if (bo < this->check.attractiveBondOrderStop && dist > this->check.attractiveDistanceStop * r12cov) { return false; } } - for (unsigned int i = 0; i < this->dissociationList.size() / 2; i++) { - const auto& l = this->dissociationList[2 * i]; - const auto& r = this->dissociationList[2 * i + 1]; - const auto bo = bondOrders.getOrder(l, r); + // dissociations + for (const auto& [left, right] : dissociations) { + double bo = 0.0; + for (const auto& l : left) { + for (const auto& r : right) { + bo += bondOrders.getOrder(l, r); + } + } if (bo > this->check.repulsiveBondOrderStop) { return false; } } - return done; + return true; } PositionCollection NtOptimizer2::extractTsGuess() const { @@ -353,16 +477,33 @@ PositionCollection NtOptimizer2::extractTsGuess() const { if (maximaList.empty()) { throw std::runtime_error("No transition state guess was found in Newton Trajectory scan."); } - // Extract TS guess from point with highest energy - double maxValue = -std::numeric_limits::max(); - int maxIndex = -1; - for (auto& i : maximaList) { - if (_values[i] > maxValue) { - maxIndex = i; - maxValue = _values[i]; + // Extract TS guess according to given criterion + if (extractionCriterion == ntExtractFirst) { + return _trajectory[maximaList.back()]; // back because maximalist starts from back + } + // get the highest value if wished or coordinate was never stopped + if (extractionCriterion == ntExtractHighest || _firstCoordinateReachedIndex == -1) { + double maxValue = -std::numeric_limits::max(); + int maxIndex = -1; + for (auto& i : maximaList) { + if (_values[i] > maxValue) { + maxIndex = i; + maxValue = _values[i]; + } + } + return _trajectory[maxIndex]; + } + int firstMaximumAfterCoordinateReached = -1; + // maximaList lists maxima from back + for (auto& maximum : maximaList) { + if (maximum < _firstCoordinateReachedIndex) { + return _trajectory[maximum]; } + firstMaximumAfterCoordinateReached = maximum; } - return _trajectory[maxIndex]; + assert(firstMaximumAfterCoordinateReached > -1); + // got no maximum before coordinate was ended, so we pick the first after the coordinate was reached + return _trajectory[firstMaximumAfterCoordinateReached]; } void NtOptimizer2::updateCoordinates(PositionCollection& coordinates, const AtomCollection& atoms, @@ -406,7 +547,6 @@ void NtOptimizer2::addSettingsDescriptors(UniversalSettings::DescriptorCollectio } void NtOptimizer2::setSettings(const Settings& settings) { - // optimizer.applySettings(settings); if (!settings.valid()) { settings.throwIncorrectSettings(); } @@ -423,6 +563,7 @@ void NtOptimizer2::setSettings(const Settings& settings) { this->numberOfMicroCycles = settings.getInt(NtOptimizer2::ntNumberOfMicroCycles); this->filterPasses = settings.getInt(NtOptimizer2::ntFilterPasses); this->fixedAtoms = settings.getIntList(NtOptimizer2::ntFixedAtomsKey); + this->extractionCriterion = settings.getString(NtOptimizer2::ntExtractionCriterion); // Check whether constraints and coordinate transformations are both switched on: if (!this->fixedAtoms.empty() && this->coordinateSystem != CoordinateSystem::Cartesian) { @@ -490,5 +631,81 @@ const std::vector>& NtOptimizer2::getConstraintsMap() { const std::vector& NtOptimizer2::getReactiveAtomsList() { return _reactiveAtomsList; } + +namespace NtUtils { +std::vector> connectedNuclei(std::vector indices, const BondOrderCollection& bondOrders) { + std::sort(indices.begin(), indices.end()); // makes avoiding double counting easier + std::vector> result; + for (const auto& i : indices) { + // check if i is already in one of the entries + if (std::any_of(result.begin(), result.end(), [&](std::vector bondedIndices) { + return std::find(bondedIndices.begin(), bondedIndices.end(), i) != bondedIndices.end(); + })) { + continue; + } + if (i == indices.back()) { + // we reached last entry and this has never been included in the others, so it must be its own thing + result.push_back({i}); + continue; + } + std::vector connected; + std::vector nucleiToVisit = {i}; + while (!nucleiToVisit.empty()) { + // take last value in list to crawl + int current = nucleiToVisit.back(); + nucleiToVisit.pop_back(); + auto bondPartners = bondOrders.getBondPartners(current, 0.5); + for (const auto& partner : bondPartners) { + // add newly found indices to total list and list to further extend search + if (std::find(connected.begin(), connected.end(), partner) == connected.end()) { + connected.push_back(partner); + nucleiToVisit.push_back(partner); + } + } + } + // now check for all others if they are within the list of bond partners i.e. same molecule + std::vector connectedIndices = {i}; + for (const auto& j : indices) { + if (j > i && std::find(connected.begin(), connected.end(), j) != connected.end()) { + connectedIndices.push_back(j); + } + } + std::sort(connectedIndices.begin(), connectedIndices.end()); + result.push_back(connectedIndices); + } + return result; +} + +Displacement centerToCenterVector(const PositionCollection& positions, const std::vector& lhsList, + const std::vector& rhsList) { + Position lhsCenter = Position::Zero(); + Position rhsCenter = Position::Zero(); + for (const auto l : lhsList) { + lhsCenter += positions.row(l); + } + lhsCenter.array() /= static_cast(lhsList.size()); + for (const auto r : rhsList) { + rhsCenter += positions.row(r); + } + rhsCenter.array() /= static_cast(rhsList.size()); + Displacement c2c = lhsCenter - rhsCenter; + return c2c; +} + +double smallestCovalentRadius(const AtomCollection& atoms, const std::vector& indices) { + if (indices.empty()) { + throw std::runtime_error("Missing reactive atom indices"); + } + double result = std::numeric_limits::max(); + for (const auto& index : indices) { + double rad = ElementInfo::covalentRadius(atoms.getElement(index)); + if (rad < result) { + result = rad; + } + } + return result; +} + +} // namespace NtUtils } // namespace Utils } // namespace Scine diff --git a/src/Utils/Utils/GeometryOptimization/NtOptimizer2.h b/src/Utils/Utils/GeometryOptimization/NtOptimizer2.h index 17dab15..389dce5 100644 --- a/src/Utils/Utils/GeometryOptimization/NtOptimizer2.h +++ b/src/Utils/Utils/GeometryOptimization/NtOptimizer2.h @@ -37,6 +37,8 @@ class NtOptimizer2Settings; * */ class NtOptimizer2 : public Optimizer { + using ReactionMapping = std::vector, std::vector>>; + public: // Definition of Utils::Settings keys static constexpr const char* ntAssListKey = "nt_associations"; @@ -50,8 +52,12 @@ class NtOptimizer2 : public Optimizer { static constexpr const char* ntFixedNumberOfMicroCycles = "nt_fixed_number_of_micro_cycles"; static constexpr const char* ntNumberOfMicroCycles = "nt_number_of_micro_cycles"; static constexpr const char* ntFilterPasses = "nt_filter_passes"; + static constexpr const char* ntExtractionCriterion = "nt_extraction_criterion"; static constexpr const char* ntCoordinateSystemKey = "nt_coordinate_system"; static constexpr const char* ntFixedAtomsKey = "nt_constrained_atoms"; + static constexpr const char* ntExtractLastBeforeTarget = "last_maximum_before_first_target"; + static constexpr const char* ntExtractHighest = "highest_maximum"; + static constexpr const char* ntExtractFirst = "first_maximum"; /** * @brief Construct a new NtOptimizer object. * @param calculator The calculator to be used for the underlying single point/gradient calculations. @@ -129,6 +135,16 @@ class NtOptimizer2 : public Optimizer { * @brief Number of passes through a Savitzky-Golay filter before analyzing the reaction curve. */ int filterPasses = 10; + // @brief possible options for extraction + const std::vector possibleExtractionOptions = { + ntExtractLastBeforeTarget, + ntExtractHighest, + ntExtractFirst, + }; + /** + * @brief Criterion to extract a TS guess from the trajectory + */ + std::string extractionCriterion = possibleExtractionOptions.front(); /** * @brief The special convergence settings for this optimizer. */ @@ -181,25 +197,24 @@ class NtOptimizer2 : public Optimizer { * @throw std::logic_error if invalid index */ void sanityCheck(const AtomCollection& atoms) const; - /** - * @brief Calculates vector from geometric centers of lhs and rhs. - * - * @param positions The current positions. - * @return Displacement The vector from rhs to lhs - */ - Displacement centerToCenterVector(const PositionCollection& positions) const; /** * @brief The function evaluating all artificial forces for the given set of atoms. * * @param atoms The current atoms. * @param energy The current energy in hartree. * @param gradients The gradient to be updated (in a.u). - * @param bondOrders The bond orders used to determine contraints + * @param bondOrders The bond orders used to determine constraints * @param addForce Add external force. */ void updateGradients(const AtomCollection& atoms, const double& energy, GradientCollection& gradients, - const BondOrderCollection& bondOrders, bool addForce = false) const; - + const BondOrderCollection& bondOrders, int cycle, bool addForce = false); + /** + * @brief Avoids gradient manipulations of the reactive atoms + * + * @param positions The current positions. + * @param gradients The to be altered gradients. + */ + void eliminateReactiveAtomsGradients(const PositionCollection& positions, GradientCollection& gradients) const; /** * @brief The function determining whether the optimization is converged. * @@ -221,6 +236,14 @@ class NtOptimizer2 : public Optimizer { * @param gradients The current gradients for the SD step. */ void updateCoordinates(PositionCollection& coordinates, const AtomCollection& atoms, const GradientCollection& gradients) const; + /** + * @brief Extracts the map for each associations and dissociations including eta bond considerations. + * + * @param bondOrders The bonds defining separate molecules. + * + * @return reactionMaps The reaction mapping of both associative and dissociative reaction coordinates + */ + std::pair inferReactions(const BondOrderCollection& bondOrders) const; /** * @brief Set the reactive atoms list by combining the sorted association and dissociation list and * make the combination unique. @@ -241,8 +264,41 @@ class NtOptimizer2 : public Optimizer { // @brief A mapping of atom indices to the constraints they are involved in. std::vector> _constraintsMap; Core::Calculator& _calculator; + int _firstCoordinateReachedIndex = -1; }; +namespace NtUtils { +/** + * @brief Gives list of list, with each sublist containing the values of indices that are part of the same molecules + * based on the given bond orders. + * + * @param indices The nuclei indices that shall be sorted into molecules. + * @param bondOrders The bonds defining separate molecules. + * + * @return molecules Each sublist representing a molecule. + */ +std::vector> connectedNuclei(std::vector indices, const BondOrderCollection& bondOrders); +/** + * @brief Calculates vector from geometric centers of lhsList and rhsList. + * + * @note opposite direction to identical method in NtOptimizer due to legacy reasons. + * + * @param positions The current positions. + * @return Displacement The vector from lhs to rhs + */ +Displacement centerToCenterVector(const PositionCollection& positions, const std::vector& lhsList, + const std::vector& rhsList); +/** + * @brief Returns the smallest covalent radius from the given indices and atoms. + * + * @param atoms The atoms to pick from. + * @param indices The list of indices of possible atoms. + * + * @return covRad The smallest covalent radius in Bohr + */ +double smallestCovalentRadius(const AtomCollection& atoms, const std::vector& indices); + +} // namespace NtUtils } // namespace Utils } // namespace Scine diff --git a/src/Utils/Utils/GeometryOptimization/NtOptimizer2Settings.h b/src/Utils/Utils/GeometryOptimization/NtOptimizer2Settings.h index 66bdba0..4b6b3bd 100644 --- a/src/Utils/Utils/GeometryOptimization/NtOptimizer2Settings.h +++ b/src/Utils/Utils/GeometryOptimization/NtOptimizer2Settings.h @@ -84,6 +84,13 @@ class NtOptimizer2Settings : public Settings { nt_filter_passes.setMinimum(0); this->_fields.push_back(NtOptimizer2::ntFilterPasses, nt_filter_passes); + UniversalSettings::OptionListDescriptor nt_extraction("Sets the TS guess extraction criterion."); + for (const auto& criterion : nt.possibleExtractionOptions) { + nt_extraction.addOption(criterion); + } + nt_extraction.setDefaultOption(nt.possibleExtractionOptions.front()); + this->_fields.push_back(NtOptimizer2::ntExtractionCriterion, nt_extraction); + UniversalSettings::OptionListDescriptor nt_coordinate_system("Set the coordinate system."); nt_coordinate_system.addOption("internal"); nt_coordinate_system.addOption("cartesianWithoutRotTrans"); diff --git a/src/Utils/Utils/GeometryOptimization/NtOptimizerSettings.h b/src/Utils/Utils/GeometryOptimization/NtOptimizerSettings.h index 947acfb..c75e262 100644 --- a/src/Utils/Utils/GeometryOptimization/NtOptimizerSettings.h +++ b/src/Utils/Utils/GeometryOptimization/NtOptimizerSettings.h @@ -94,6 +94,13 @@ class NtOptimizerSettings : public Settings { nt_filter_passes.setMinimum(0); this->_fields.push_back(NtOptimizer::ntFilterPasses, nt_filter_passes); + UniversalSettings::OptionListDescriptor nt_extraction("Sets the TS guess extraction criterion."); + for (const auto& criterion : nt.possibleExtractionOptions) { + nt_extraction.addOption(criterion); + } + nt_extraction.setDefaultOption(nt.possibleExtractionOptions.front()); + this->_fields.push_back(NtOptimizer::ntExtractionCriterion, nt_extraction); + UniversalSettings::OptionListDescriptor nt_coordinate_system("Set the coordinate system."); nt_coordinate_system.addOption("internal"); nt_coordinate_system.addOption("cartesianWithoutRotTrans"); diff --git a/src/Utils/Utils/IO/MolecularTrajectoryIO.cpp b/src/Utils/Utils/IO/MolecularTrajectoryIO.cpp index 39328d4..f2f03ee 100644 --- a/src/Utils/Utils/IO/MolecularTrajectoryIO.cpp +++ b/src/Utils/Utils/IO/MolecularTrajectoryIO.cpp @@ -76,8 +76,13 @@ void MolecularTrajectoryIO::writeBinary(std::ostream& out, const MolecularTrajec void MolecularTrajectoryIO::writeXYZ(std::ostream& out, const MolecularTrajectory& m) { out.imbue(std::locale("C")); const auto& elements = m.getElementTypes(); + bool hasEnergies = !m.getEnergies().empty(); for (int i = 0; i < m.size(); ++i) { - out << m.molecularSize() << endl << endl; + out << m.molecularSize() << endl; + if (hasEnergies) { + out << m.getEnergies()[i]; + } + out << endl; for (int j = 0; j < m.molecularSize(); ++j) { writeXYZLine(out, elements[j], m[i].row(j)); } @@ -161,6 +166,7 @@ MolecularTrajectory MolecularTrajectoryIO::readXYZ(std::istream& in) { bool firstDone = false; ElementTypeCollection elements; MolecularTrajectory m; + std::vector energies; while (!in.eof()) { int nAtoms; in >> nAtoms; @@ -170,7 +176,18 @@ MolecularTrajectory MolecularTrajectoryIO::readXYZ(std::istream& in) { std::string unnecessaryLine; getline(in, unnecessaryLine); - getline(in, unnecessaryLine); + std::string commentLine; + getline(in, commentLine); + bool gotEnergy = !commentLine.empty() && commentLine.find_first_not_of("-.0123456789") == string::npos; + if (gotEnergy) { + // in case we got funky comment, that really isn't a double + try { + energies.push_back(std::stod(commentLine)); + } + catch (...) { + ; + } + } std::string element; PositionCollection s(nAtoms, 3); for (int i = 0; i < nAtoms; ++i) { @@ -195,6 +212,9 @@ MolecularTrajectory MolecularTrajectoryIO::readXYZ(std::istream& in) { m.push_back(s * Constants::bohr_per_angstrom); } m.setElementTypes(elements); + if (energies.size() == static_cast(m.size())) { + m.setEnergies(energies); + } return m; } diff --git a/src/Utils/Utils/MolecularTrajectory.cpp b/src/Utils/Utils/MolecularTrajectory.cpp index 95e8994..63aa3a6 100644 --- a/src/Utils/Utils/MolecularTrajectory.cpp +++ b/src/Utils/Utils/MolecularTrajectory.cpp @@ -10,6 +10,20 @@ namespace Scine { namespace Utils { +MolecularTrajectory::MolecularTrajectory(const ElementTypeCollection& elements) { + elements_ = elements; +} + +MolecularTrajectory::MolecularTrajectory(double minimumRmsdForAddition) + : MolecularTrajectory(elements_, minimumRmsdForAddition) { +} + +MolecularTrajectory::MolecularTrajectory(const ElementTypeCollection& elements, double minimumRmsdForAddition) { + elements_ = elements; + minMeanSquareDeviation_ = minimumRmsdForAddition * minimumRmsdForAddition; + respectMinRmsd_ = true; +} + void MolecularTrajectory::setElementType(int index, ElementType e) { elements_[index] = e; } @@ -85,7 +99,9 @@ void MolecularTrajectory::push_back(PositionCollection p) { if (!pbcs_.empty()) { throw std::runtime_error("Pbc container is not empty. Clear pbcs first, before a structure only is added."); } - structureVector_.push_back(std::move(p)); + if (additionOfPositionCollectionIsAllowedBasedOnRmsd(p)) { + structureVector_.push_back(std::move(p)); + } } void MolecularTrajectory::push_back(PositionCollection p, double e) { @@ -94,8 +110,10 @@ void MolecularTrajectory::push_back(PositionCollection p, double e) { throw std::runtime_error( "The number of energies does not match the number of structures in this molecular trajectory."); } - structureVector_.push_back(std::move(p)); - energies_.push_back(e); + if (additionOfPositionCollectionIsAllowedBasedOnRmsd(p)) { + structureVector_.push_back(std::move(p)); + energies_.push_back(e); + } } void MolecularTrajectory::push_back(PositionCollection p, PeriodicBoundaries pbc) { @@ -104,8 +122,10 @@ void MolecularTrajectory::push_back(PositionCollection p, PeriodicBoundaries pbc throw std::runtime_error( "The number of pbcs does not match the number of structures in this molecular trajectory."); } - structureVector_.push_back(std::move(p)); - pbcs_.push_back(pbc.getCellMatrix()); + if (additionOfPositionCollectionIsAllowedBasedOnRmsd(p)) { + structureVector_.push_back(std::move(p)); + pbcs_.push_back(pbc.getCellMatrix()); + } } void MolecularTrajectory::push_back(PositionCollection p, double e, PeriodicBoundaries pbc) { @@ -119,9 +139,11 @@ void MolecularTrajectory::push_back(PositionCollection p, double e, PeriodicBoun throw std::runtime_error( "The number of pbcs does not match the number of structures in this molecular trajectory."); } - structureVector_.push_back(std::move(p)); - energies_.push_back(e); - pbcs_.push_back(pbc.getCellMatrix()); + if (additionOfPositionCollectionIsAllowedBasedOnRmsd(p)) { + structureVector_.push_back(std::move(p)); + energies_.push_back(e); + pbcs_.push_back(pbc.getCellMatrix()); + } } bool MolecularTrajectory::empty() const { @@ -237,5 +259,14 @@ bool MolecularTrajectory::additionOfPositionCollectionIsAllowed(const PositionCo return validGivenThatElementCollectionIsAlreadySet || validGivenThatElementCollectionIsNotSet; } +bool MolecularTrajectory::additionOfPositionCollectionIsAllowedBasedOnRmsd(const PositionCollection& p) const { + if (!respectMinRmsd_ || structureVector_.empty()) { + return true; + } + PositionCollection lastPos = structureVector_.back(); + double meanSquareDeviation = (lastPos - p).rowwise().squaredNorm().sum() / lastPos.rows(); + return meanSquareDeviation > minMeanSquareDeviation_; +} + } /* namespace Utils */ } /* namespace Scine */ diff --git a/src/Utils/Utils/MolecularTrajectory.h b/src/Utils/Utils/MolecularTrajectory.h index fd8857e..7e31f22 100644 --- a/src/Utils/Utils/MolecularTrajectory.h +++ b/src/Utils/Utils/MolecularTrajectory.h @@ -33,6 +33,9 @@ class MolecularTrajectory { /// @brief Default constructor. MolecularTrajectory() = default; + explicit MolecularTrajectory(const ElementTypeCollection& elements); + explicit MolecularTrajectory(double minimumRmsdForAddition); + explicit MolecularTrajectory(const ElementTypeCollection& elements, double minimumRmsdForAddition); /** * @brief Set a single element type. * @param index The index of the atom to be altered @@ -282,10 +285,13 @@ class MolecularTrajectory { * @return false If the addition is not allowed. */ bool additionOfPositionCollectionIsAllowed(const PositionCollection& p) const; + bool additionOfPositionCollectionIsAllowedBasedOnRmsd(const PositionCollection& p) const; Container structureVector_; ElementTypeCollection elements_; EnergyContainer energies_; PbcContainer pbcs_; + double minMeanSquareDeviation_; + bool respectMinRmsd_ = false; }; } /* namespace Utils */ diff --git a/src/Utils/Utils/Scf/OrbitalPerturbation/UniqueRandomNumbersGenerator.h b/src/Utils/Utils/Scf/OrbitalPerturbation/UniqueRandomNumbersGenerator.h index 04cddc3..f0d8347 100644 --- a/src/Utils/Utils/Scf/OrbitalPerturbation/UniqueRandomNumbersGenerator.h +++ b/src/Utils/Utils/Scf/OrbitalPerturbation/UniqueRandomNumbersGenerator.h @@ -11,6 +11,7 @@ #include #include #include +#include #include namespace Scine { @@ -62,7 +63,9 @@ std::vector UniqueRandomNumbersGenerator::generate(uns std::vector shuffler(max_ - min_ + 1); std::iota(shuffler.begin(), shuffler.end(), min_); - std::random_shuffle(shuffler.begin(), shuffler.end()); + std::random_device rd; + std::mt19937 g(rd()); + std::shuffle(shuffler.begin(), shuffler.end(), g); return {shuffler.begin(), shuffler.begin() + N}; } diff --git a/src/Utils/Utils/Solvation/RandomIndexGenerator.cpp b/src/Utils/Utils/Solvation/RandomIndexGenerator.cpp index 93af42d..92ad76c 100644 --- a/src/Utils/Utils/Solvation/RandomIndexGenerator.cpp +++ b/src/Utils/Utils/Solvation/RandomIndexGenerator.cpp @@ -15,6 +15,14 @@ SoluteSolventComplex::RandomIndexGenerator::RandomIndexGenerator(int size, int s distribution_ = std::uniform_int_distribution<>(0, size - 1); } +void SoluteSolventComplex::RandomIndexGenerator::setSeed(int seed) { + generator_ = std::mt19937(seed); +} + +void SoluteSolventComplex::RandomIndexGenerator::setSize(int size) { + distribution_ = std::uniform_int_distribution<>(0, size - 1); +} + int SoluteSolventComplex::RandomIndexGenerator::next() { return distribution_(generator_); } diff --git a/src/Utils/Utils/Solvation/RandomIndexGenerator.h b/src/Utils/Utils/Solvation/RandomIndexGenerator.h index a3e4e88..a7138e4 100644 --- a/src/Utils/Utils/Solvation/RandomIndexGenerator.h +++ b/src/Utils/Utils/Solvation/RandomIndexGenerator.h @@ -21,12 +21,27 @@ namespace SoluteSolventComplex { */ class RandomIndexGenerator { public: + /** + * @brief Default empty constructor + */ + RandomIndexGenerator() = default; /** * @brief Initialising the RandomIndexGenerator, assigning size and seed to its generator and distribution * @param size Size of the wanted distribution, meaning the size of the related vector. * @param seed Seed for the generator. */ RandomIndexGenerator(int size, int seed); + + /** + * @brief Set seed for generator + * @param seed Seed for generator. + */ + void setSeed(int seed); + /** + * @brief Set size of wanted distribution. + * @param size Size of the wanted distribution. + */ + void setSize(int size); /** * @brief Generate next random index from initialised RandomIndexGenerator. * @return Integer between 0 and size - 1 with uniform distribution. diff --git a/src/Utils/Utils/Solvation/SoluteSolventComplex.cpp b/src/Utils/Utils/Solvation/SoluteSolventComplex.cpp index 3d7eaaf..2b465b6 100644 --- a/src/Utils/Utils/Solvation/SoluteSolventComplex.cpp +++ b/src/Utils/Utils/Solvation/SoluteSolventComplex.cpp @@ -14,44 +14,56 @@ namespace Utils { namespace SoluteSolventComplex { std::vector> solvate(const AtomCollection& soluteComplex, int soluteSize, - const AtomCollection& solvent, int numSolvents, int seed, int resolution, - double solventOffset, double maxDistance, double stepSize, - int numRotamers, bool strategicSolv, double coverageThreshold) { - return SoluteSolventComplex::solvate(soluteComplex, soluteSize, solvent, numSolvents, std::numeric_limits::max(), - seed, resolution, solventOffset, maxDistance, stepSize, numRotamers, - strategicSolv, coverageThreshold); + const AtomCollection& solvent, int numSolvents, int seed, + SolventPlacementSettings placementSettings) { + auto resultTuple = + SoluteSolventComplex::solvate(soluteComplex, soluteSize, std::vector{solvent}, std::vector{1}, + numSolvents, std::numeric_limits::max(), seed, placementSettings); + + return std::get<0>(resultTuple); } -std::vector> solvateShells(const AtomCollection& soluteComplex, int soluteSize, - const AtomCollection& solvent, int numShells, int seed, int resolution, - double solventOffset, double maxDistance, double stepSize, - int numRotamers, bool strategicSolv, double coverageThreshold) { - return SoluteSolventComplex::solvate(soluteComplex, soluteSize, solvent, std::numeric_limits::max(), numShells, - seed, resolution, solventOffset, maxDistance, stepSize, numRotamers, - strategicSolv, coverageThreshold); +std::tuple>, std::vector>> +solvateMix(const AtomCollection& soluteComplex, int soluteSize, const std::vector& solvents, + const std::vector& solventRatios, int numSolvents, int seed, SolventPlacementSettings placementSettings) { + return SoluteSolventComplex::solvate(soluteComplex, soluteSize, solvents, solventRatios, numSolvents, + std::numeric_limits::max(), seed, placementSettings); } -std::vector> -solvate(const AtomCollection& soluteComplex, const int soluteSize, const AtomCollection& solvent, const int numSolvents, - const int numShells, const int seed, const int resolution, const double solventOffset, const double maxDistance, - const double stepSize, const int numRotamers, const bool strategicSolv, const double coverageThreshold) { - auto solventSites = MolecularSurface::getPrunedMolecularSurface(solvent, resolution); +std::vector> solvateShells(const AtomCollection& soluteComplex, int soluteSize, + const AtomCollection& solvent, int numShells, int seed, + SolventPlacementSettings placementSettings) { + auto resultTuple = + SoluteSolventComplex::solvate(soluteComplex, soluteSize, std::vector{solvent}, std::vector{1}, + std::numeric_limits::max(), numShells, seed, placementSettings); - double dMin = solventOffset; - double dMax = maxDistance; + return std::get<0>(resultTuple); +} + +std::tuple>, std::vector>> +solvateShellsMix(const AtomCollection& soluteComplex, int soluteSize, const std::vector& solvents, + const std::vector& solventRatios, int numShells, int seed, SolventPlacementSettings placementSettings) { + return SoluteSolventComplex::solvate(soluteComplex, soluteSize, solvents, solventRatios, + std::numeric_limits::max(), numShells, seed, placementSettings); +} +std::tuple>, std::vector>> +solvate(const AtomCollection& soluteComplex, int soluteSize, const std::vector& solvents, + const std::vector& solventRatios, int numSolvents, int numShells, int seed, + SolventPlacementSettings placementSettings) { + double dMin = placementSettings.solventOffset; + double dMax = placementSettings.maxDistance; auto complex = soluteComplex; std::vector solventVector; + std::vector solventIndexVector; std::vector> shellVector; + std::vector> solventIndexShellVector; int numAddedSolvents = 0; int numAddedShells = 0; int soluteStart = 0; int soluteEnd = soluteSize; // random number generator with given seed std::mt19937 shuffleGen(seed); - // set up index generator with seed and solvent - SoluteSolventComplex::RandomIndexGenerator randSolventSite(solventSites.size(), seed); - // bool to remember if previous addition was successful or not bool prevAddition = false; int reCalcCounter = 0; @@ -59,63 +71,98 @@ solvate(const AtomCollection& soluteComplex, const int soluteSize, const AtomCol // trigger to force recalculation of the visible surface int reCalcTrigger = 1; // calculate first surface before adding anything - auto visibleComplexSites = MolecularSurface::getVisibleMolecularSurface(complex, soluteStart, soluteEnd, resolution); - // save original number of sites and initialize coverag - double sizeVisibleComplexSites = visibleComplexSites.size(); + auto visibleComplexSites = + MolecularSurface::getVisibleMolecularSurface(complex, soluteStart, soluteEnd, placementSettings.resolution); + // save original number of sites and initialize coverage + auto sizeVisibleComplexSites = double(visibleComplexSites.size()); bool resetSize = false; double coverage = 0.0; + // Create solvent index list and shuffle it; hard limit of 1e6 solvent molecules + if (numSolvents == std::numeric_limits::max()) { + numSolvents = 1000000; + } + bool multipleSolvents = true; + // Declaration of solvent related variables + AtomCollection solvent; + std::vector flatSolventIndexVector; + std::vector solventSites; + SoluteSolventComplex::RandomIndexGenerator randSolventSite; + // Check, if only one solvent in solvents + // Set solvent and solvent sites, if only one solvent is given + if (solvents.size() == 1) { + solvent = solvents.at(0); + solventSites = MolecularSurface::getPrunedMolecularSurface(solvent, placementSettings.resolution); + // set index generator with seed and solvent + randSolventSite.setSize(solventSites.size()); + randSolventSite.setSeed(seed); + multipleSolvents = false; + } + else { + flatSolventIndexVector = getSolventIndices(numSolvents, solventRatios, solvents.size()); + std::shuffle(std::begin(flatSolventIndexVector), std::end(flatSolventIndexVector), shuffleGen); + } // loop till required number of solvent molecules has been added while (numAddedSolvents < numSolvents && numAddedShells < numShells) { bool addition = false; // check, if previous addition round was successful; if yes, recalculate surface if (prevAddition && (numAddedSolvents % reCalcTrigger == 0)) { - visibleComplexSites = MolecularSurface::getVisibleMolecularSurface(complex, soluteStart, soluteEnd, resolution); + visibleComplexSites = + MolecularSurface::getVisibleMolecularSurface(complex, soluteStart, soluteEnd, placementSettings.resolution); // recalculate number of visible complex sites to calculate coverage if (resetSize) { - sizeVisibleComplexSites = visibleComplexSites.size(); + sizeVisibleComplexSites = double(visibleComplexSites.size()); resetSize = false; } // calculate coverage of current solute coverage = 1.0 - visibleComplexSites.size() / sizeVisibleComplexSites; - if (strategicSolv) { + if (placementSettings.strategicSolv) { reCalcTrigger = SoluteSolventComplex::solvationStrategy(visibleComplexSites.size()); } // shuffle visible solute sites std::shuffle(std::begin(visibleComplexSites), std::end(visibleComplexSites), shuffleGen); reCalcCounter += 1; } + // choose solvent here according to index + if (multipleSolvents) { + solvent = solvents.at(flatSolventIndexVector.at(numAddedSolvents)); + solventSites = MolecularSurface::getPrunedMolecularSurface(solvent, placementSettings.resolution); + // set index generator with seed and solvent + randSolventSite.setSize(solventSites.size()); + randSolventSite.setSeed(seed); + } // choose random solvent surface index int randSolventSiteIndex = randSolventSite.next(); // loop over visible complex sites and try to add solvent molecule for (const auto& visibleComplexSite : visibleComplexSites) { addition = SoluteSolventComplex::add(complex, solvent, visibleComplexSite, solventSites.at(randSolventSiteIndex), - dMin, dMax, stepSize, numRotamers); + dMin, dMax, placementSettings.stepSize, placementSettings.numRotamers); // if addition was successful, break loop over complex sites and try to add next solvent molecule if (addition) { - numAddedSolvents += 1; - prevAddition = true; AtomCollection tmpSolvent; - int tmpIndex = 0; // extract lastly added solvent and push to solvent list for (int i = complex.size() - solvent.size(); i < complex.size(); i++) { tmpSolvent.push_back(complex.at(i)); - tmpIndex += 1; } solventVector.push_back(tmpSolvent); + if (multipleSolvents) { + solventIndexVector.push_back(flatSolventIndexVector.at(numAddedSolvents)); + } + numAddedSolvents += 1; + prevAddition = true; break; } } // if addition over all sites was not successful, try next round by increasing the dMin and dMax values if (!addition) { dMin = dMax; - dMax += maxDistance; + dMax += placementSettings.maxDistance; prevAddition = false; failedAdd += 1; } // if no more complex sites are visible, return solventVector - if (visibleComplexSites.empty() || coverage >= coverageThreshold) { - dMin = solventOffset; - dMax = maxDistance; + if (visibleComplexSites.empty() || (placementSettings.coverageThreshold - coverage) <= 1e-12) { + dMin = placementSettings.solventOffset; + dMax = placementSettings.maxDistance; soluteStart = soluteEnd; soluteEnd = complex.size(); resetSize = true; @@ -124,14 +171,18 @@ solvate(const AtomCollection& soluteComplex, const int soluteSize, const AtomCol coverage = 0.0; numAddedShells += 1; shellVector.push_back(solventVector); + solventIndexShellVector.push_back(solventIndexVector); solventVector.clear(); + solventIndexVector.clear(); } } // only add last solventVector, if it has not been added due to a full shell if (!visibleComplexSites.empty()) { shellVector.push_back(solventVector); + solventIndexShellVector.push_back(solventIndexVector); } - return shellVector; + + return {shellVector, solventIndexShellVector}; } AtomCollection mergeAtomCollectionVector(const std::vector& atomCollList) { @@ -150,6 +201,14 @@ AtomCollection mergeSolventShellVector(const std::vector mergeSolventShellIndices(const std::vector>& solventShellIndices) { + std::vector flat; + for (const std::vector& shell : solventShellIndices) { + flat.insert(flat.end(), shell.begin(), shell.end()); + } + return flat; +} + std::vector transferSolventShellVector(const std::vector>& shellVector) { std::vector solventSizeVector; for (const auto& shell : shellVector) { @@ -257,6 +316,32 @@ int solvationStrategy(int numberSurfPoints) { return 50; } +std::vector getSolventIndices(int numSolvents, const std::vector& ratios, unsigned long int differentSolvents) { + // Check if ratios and solvents are of the same size. + if (ratios.size() != differentSolvents) { + throw std::runtime_error("The same number of ratios and solvents has to be given."); + } + std::vector solventIndices(numSolvents); + int indexedSolvents = 0; + while (indexedSolvents < numSolvents) { + // Index of solvent in solvents + int solventIndex = 0; + for (const int& ratio : ratios) { + // Break loop over ratios if number of indexed solvents exceeds number of solvents + if (indexedSolvents >= numSolvents) { + break; + } + int lastIndex = (indexedSolvents + ratio >= numSolvents) ? numSolvents : indexedSolvents + ratio; + std::fill(solventIndices.begin() + indexedSolvents, solventIndices.begin() + lastIndex, solventIndex); + // Increase index of solvent in solvents + solventIndex += 1; + // Add number of indexed solvents to count + indexedSolvents += lastIndex - indexedSolvents; + } + } + return solventIndices; +} + bool checkDistances(const AtomCollection& atoms1, const AtomCollection& atoms2) { // loop over atom collection, assuming atoms2 is smaller than atoms1 for (const auto& atom2 : atoms2) { diff --git a/src/Utils/Utils/Solvation/SoluteSolventComplex.h b/src/Utils/Utils/Solvation/SoluteSolventComplex.h index 2ee459e..7f85bdb 100644 --- a/src/Utils/Utils/Solvation/SoluteSolventComplex.h +++ b/src/Utils/Utils/Solvation/SoluteSolventComplex.h @@ -25,55 +25,104 @@ namespace Utils { */ namespace SoluteSolventComplex { /** - * @brief Add systematically a number of solvents to solute + * @brief Structure of settings for solvent placement * - * @param soluteComplex AtomCollection of the solute (solute complex, if already partially solvated) - * @param soluteSize Size of the solute, assuming the solute starts with index 0. - * @param solvent AtomCollection of solvent. - * @param numSolvents Number of solvent molecules which should be added. - * @param seed Seed for shuffeling the surface sites and choosing the surface site of the solvent * @param resolution Number of surface sites per atom * @param solventOffset Initial distance to arrange solvent molecule. * @param maxDistance Maximum distance below it should be attempted to add a solvent molecule before the next surface * point is tested out. - * @param stepSize Increament of solvent offset if no solvent molcule could be added. + * @param stepSize Increment of solvent offset if no solvent molecule could be added. * @param numRotamers Number of rotations to be tried for adding solvent. * @param strategicSolv Bool to turn on or off the solvationStrategy function. Reduces the number of surface * recalculations. * @param coverageThreshold Ratio of visible solvent sites to be covered by solvents per total number of visible solvent * sites for the given solute. Once the threshold is reached the shell is considered complete. - * @return A solvent shell vector. Each entry is one vector of atom collections, containg the solvent molecules of this - * shell. */ -std::vector> -solvate(const AtomCollection& soluteComplex, int soluteSize, const AtomCollection& solvent, int numSolvents, int seed, - int resolution = 32, double solventOffset = 0.0, double maxDistance = 10.0, double stepSize = 0.25, - int numRotamers = 3, bool strategicSolv = false, double coverageThreshold = 1.0); +struct SolventPlacementSettings { + int resolution = 18; + double solventOffset = 0.0; + double maxDistance = 5.0; + double stepSize = 1.0; + int numRotamers = 3; + bool strategicSolv = true; + double coverageThreshold = 0.85; +}; +/** + * @brief Add systematically a number of one type of solvent to solute + * + * @param soluteComplex AtomCollection of the solute (solute complex, if already partially solvated) + * @param soluteSize Size of the solute, assuming the solute starts with index 0. + * @param solvent AtomCollection of the solvent. + * @param numSolvents Number of solvent molecules which should be added. + * @param seed Seed for shuffling the surface sites and choosing the surface site of the solvent. + * @param placementSettings Settings for the placement of one solvent. + * + * @return A solvent shell vector. Each entry is one vector of atom collections, containing the solvent molecules of + * this shell. + */ +std::vector> solvate(const AtomCollection& soluteComplex, int soluteSize, + const AtomCollection& solvent, int numSolvents, int seed, + SolventPlacementSettings placementSettings); /** - * @brief Add number of solvent shells to solute. + * @brief Add systematically a number of different types of solvent with the given ratio to the solute. + * + * @param soluteComplex AtomCollection of the solute (solute complex, if already partially solvated) + * @param soluteSize Size of the solute, assuming the solute starts with index 0. + * @param solvents A vector of AtomCollections of different solvents. + * @param solventRatios A vector of ratios of the different solvents. + * @param numSolvents Number of solvent molecules which should be added. + * @param seed Seed for shuffling the surface sites, shuffling the occurrence of the different solvents and choosing + * the surface site of the solvent. + * @param placementSettings Settings for the placement of one solvent. + * + * @return A tuple of a solvent shell vector and a solvent shell indices vector (to determine where which solvent is). + * Each entry of the solvent shell vector is one vector of atom collections, containing the solvent molecules of + * this shell. + * Each entry of the solvent shell indices vector is one vector of indices, containing the solvent indices of the + * solvents given in the original input. + */ +std::tuple>, std::vector>> +solvateMix(const AtomCollection& soluteComplex, int soluteSize, const std::vector& solvents, + const std::vector& solventRatios, int numSolvents, int seed, SolventPlacementSettings placementSettings); +/** + * @brief Add number of shells of one type of solvent to solute. * * @param soluteComplex AtomCollection of the solute (solute complex, if already partially solvated) * @param soluteSize Size of the solute, assuming the solute starts with index 0. * @param solvent AtomCollection of solvent. - * @param numShells Number of shells to be completed. One shell means that no surface point is 'visible' anymore. - * @param seed Seed for shuffeling the surface sites and choosing the surface site of the solvent - * @param resolution Number of surface sites per atom - * @param solventOffset Initial distance to arrange solvent molecule. - * @param maxDistance Maximum distance below it should be attempted to add a solvent molecule before the next surface - * point is tested out. - * @param stepSize Increament of solvent offset if no solvent molcule could be added. - * @param numRotamers Number of rotations to be tried for adding solvent. - * @param strategicSolv Bool to turn on or off the solvationStrategy function. Reduces the number of surface - * recalculations. - * @param coverageThreshold Ratio of visible solvent sites to be covered by solvents per total number of visible solvent - * sites for the given solute. Once the threshold is reached the shell is considered complete. - * @return A solvent shell vector. Each entry is one vector of atom collections, containg the solvent molecules of this - * shell. + * @param numShells Number of shells to be completed by one type of solvent. One shell means that no surface point + * is 'visible' anymore. + * @param seed Seed for shuffling the surface sites and choosing the surface site of the solvent + * @param placementSettings Settings for the placement of one solvent. + * + * @return A solvent shell vector. Each entry is one vector of atom collections, containing the solvent molecules of + * this shell. */ -std::vector> -solvateShells(const AtomCollection& soluteComplex, int soluteSize, const AtomCollection& solvent, int numShells, - int seed, int resolution = 32, double solventOffset = 0.0, double maxDistance = 10.0, - double stepSize = 0.25, int numRotamers = 3, bool strategicSolv = false, double coverageThreshold = 1.0); +std::vector> solvateShells(const AtomCollection& soluteComplex, int soluteSize, + const AtomCollection& solvent, int numShells, int seed, + SolventPlacementSettings placementSettings); +/** + * @brief Add number of shells of different solvents to solute. + * + * @param soluteComplex AtomCollection of the solute (solute complex, if already partially solvated) + * @param soluteSize Size of the solute, assuming the solute starts with index 0. + * @param solvents A vector of AtomCollections of different solvents. + * @param solventRatios A vector of ratios of the different solvents. + * @param numShells Number of shells to be completed by solvents with given ratio. One shell means that no surface point + * is 'visible' anymore. + * @param seed Seed for shuffling the surface sites, shuffling the occurrence of the different solvents and choosing + * the surface site of the solvent. + * @param placementSettings Settings for the placement of one solvent. + * + * @return A tuple of a solvent shell vector and a solvent shell indices vector (to determine where which solvent is). + * Each entry of the solvent shell vector is one vector of atom collections, containing the solvent molecules of + * this shell. + * Each entry of the solvent shell indices vector is one vector of indices, containing the solvent indices of the + * solvents given in the original input. + */ +std::tuple>, std::vector>> +solvateShellsMix(const AtomCollection& soluteComplex, int soluteSize, const std::vector& solvents, + const std::vector& solventRatios, int numShells, int seed, SolventPlacementSettings placementSettings); /** * @brief Add systematically a number of solvents and a number of solvent shells to solute * @@ -82,24 +131,30 @@ solvateShells(const AtomCollection& soluteComplex, int soluteSize, const AtomCol * @param solvent AtomCollection of solvent. * @param numShells Number of shells to be completed. One shell means that no surface point is 'visible' anymore. * @param numSolvents Number of solvent molecules which should be added. - * @param seed Seed for shuffeling the surface sites and choosing the surface site of the solvent - * @param resolution Number of surface sites per atom - * @param solventOffset Initial distance to arrange solvent molecule. - * @param maxDistance Maximum distance below it should be attempted to add a solvent molecule before the next surface - * point is tested out. - * @param stepSize Increament of solvent offset if no solvent molcule could be added. - * @param numRotamers Number of rotations to be tried for adding solvent. - * @param strategicSolv Bool to turn on or off the solvationStrategy function. Reduces the number of surface - * recalculations. - * @param coverageThreshold Ratio of visible solvent sites to be covered by solvents per total number of visible solvent - * sites for the given solute. Once the threshold is reached the shell is considered complete. - * @return A solvent shell vector. Each entry is one vector of atom collections, containg the solvent molecules of this - * shell. + * @param seed Seed for shuffling the surface sites and choosing the surface site of the solvent + * @param placementSettings Settings for the placement of one solvent. + * + * @return A tuple of a solvent shell vector and a solvent shell indices vector (to determine where which solvent is). + * Each entry of the solvent shell vector is one vector of atom collections, containing the solvent molecules of + * this shell. + * Each entry of the solvent shell indices vector is one vector of indices, containing the solvent indices of the + * solvents given in the original input. */ -std::vector> solvate(const AtomCollection& soluteComplex, int soluteSize, - const AtomCollection& solvent, int numSolvents, int numShells, int seed, - int resolution, double solventOffset, double maxDistance, double stepSize, - int numRotamers, bool strategicSolv, double coverageThreshold); +std::tuple>, std::vector>> +solvate(const AtomCollection& soluteComplex, int soluteSize, const std::vector& solvents, + const std::vector& solventRatios, int numSolvents, int numShells, int seed, + SolventPlacementSettings placementSettings); +/** + * @brief Create a vector of indices of the solvent molecules. + * + * @param numSolvents Number of solvent molecules which should be added. + * @param ratios A vector of ratios of the different solvents. + * @param differentSolvents Number of different solvents. + * + * @return A vector with indices of solvents. + */ +std::vector getSolventIndices(int numSolvents, const std::vector& ratios, unsigned long int differentSolvents); + /** * @brief Analyze given solute-solvent complex to return solvent shell vector. Basically reverse solvate function, * giving the solvent shell vector of a given cluster. @@ -112,12 +167,13 @@ std::vector> solvate(const AtomCollection& soluteCom * recalculations when a solvent has been added. * @param coverageThreshold Ratio of visible solvent sites to be covered by solvents per total number of visible solvent * sites for the given solute. Once the threshold is reached the shell is considered complete. - * @return A sovlent shell vector in a solute solvent complex. Each entry is one vector of atom collections, containg + + * @return A solvent shell vector in a solute solvent complex. Each entry is one vector of atom collections, containing * the solvent molecules of this shell. Threshold allows of more loose definition of solvent shell vector. */ std::vector> giveSolventShellVector(const AtomCollection& complex, int soluteSize, const std::vector& solventSizeVector, - int resolution, Core::Log& log, bool strategicSolv = true, double coverageThreshold = 1); + int resolution, Core::Log& log, bool strategicSolv = true, double coverageThreshold = 1.0); /** * @brief Merge a vector of atom collections to one atom collection. * @@ -132,6 +188,15 @@ AtomCollection mergeAtomCollectionVector(const std::vector& atom * @return One AtomCollection, ordered as in the given vector. */ AtomCollection mergeSolventShellVector(const std::vector>& shellVector); +/** + * @brief Flatten the solvent shell indices vector to one dimension + * + * @param solventShellIndices A vector of vectors of indices. + * + * @return A vector of indices to track which input solvent is where. + */ +std::vector mergeSolventShellIndices(const std::vector>& solventShellIndices); + /** * @brief Transform solvent shell vector to solvent size vector. Needed for giveSolventShellVector calculation. * @@ -177,7 +242,7 @@ PositionCollection arrange(const Position& surfPoint1, const Position& surfNorma * @param complexSurfSite SurfaceSite of complex along which additive should be arranged. * @param additiveSurfSite SurfaceSite of additive along which additive should be arranged. * @param minDistance Minimal distance between two surface sites attempted to add. - * @param maxDistance Maxium distance between two surface sites attempted to add. + * @param maxDistance Maximum distance between two surface sites attempted to add. * @param incrementDistance Distance increment, if addition was unsuccessful and distance is still smaller than * maxDistance. * @param numRotationAttempts Number of rotations of additive to be tried to add additive. diff --git a/src/Utils/Utils/UniversalSettings/ValueCollection.cpp b/src/Utils/Utils/UniversalSettings/ValueCollection.cpp index 74ea40a..914ecc6 100644 --- a/src/Utils/Utils/UniversalSettings/ValueCollection.cpp +++ b/src/Utils/Utils/UniversalSettings/ValueCollection.cpp @@ -315,30 +315,44 @@ inline bool equalToDefault(const std::string& key, const GenericValue& value) { } bool operator==(const ValueCollection& lhs, const ValueCollection& rhs) { - if (lhs.empty() && rhs.empty()) { - return true; + return lhs.getDivergingKeys(rhs, true).empty(); +} + +bool operator!=(const ValueCollection& lhs, const ValueCollection& rhs) { + return !(lhs == rhs); +} + +std::vector ValueCollection::getDivergingKeys(const ValueCollection& otherSettings, + bool onlyFirstDivergingKey) const { + std::vector result; + if (this->empty() && otherSettings.empty()) { + return result; } - auto lhsKeys = lhs.getKeys(); - auto rhsKeys = rhs.getKeys(); + auto lhsKeys = this->getKeys(); + auto rhsKeys = otherSettings.getKeys(); for (const auto& key : lhsKeys) { // key not present or key present in both, but value is different - if (!rhs.valueExists(key) || lhs.getValue(key) != rhs.getValue(key)) { - return false; + if (!otherSettings.valueExists(key) || this->getValue(key) != otherSettings.getValue(key)) { + result.push_back(key); + if (onlyFirstDivergingKey) { + return result; + } } } // identical the other way around for (const auto& key : rhsKeys) { - if (!lhs.valueExists(key) || lhs.getValue(key) != rhs.getValue(key)) { - return false; + if (!this->valueExists(key) || this->getValue(key) != otherSettings.getValue(key)) { + result.push_back(key); + if (onlyFirstDivergingKey) { + return result; + } } } - // all checks survived -> settings are equal - return true; -} - -bool operator!=(const ValueCollection& lhs, const ValueCollection& rhs) { - return !(lhs == rhs); + // erase duplicates + std::sort(result.begin(), result.end()); + result.erase(std::unique(result.begin(), result.end()), result.end()); + return result; } } /* namespace UniversalSettings */ diff --git a/src/Utils/Utils/UniversalSettings/ValueCollection.h b/src/Utils/Utils/UniversalSettings/ValueCollection.h index 435c9e0..fe3c5d9 100644 --- a/src/Utils/Utils/UniversalSettings/ValueCollection.h +++ b/src/Utils/Utils/UniversalSettings/ValueCollection.h @@ -122,6 +122,8 @@ class ValueCollection { return std::end(values_); } + std::vector getDivergingKeys(const ValueCollection& otherSettings, bool onlyFirstDivergingKey = false) const; + private: const GenericValue& getGenericValue(const std::string& name) const; Container::iterator findName(const std::string& name);