Release 4.0.0

.gitignore

@@ -14,6 +14,8 @@ tmp_clang_*
 compile_commands.json
 .lvimrc
 *.vim
+.DS_Store
+*~
 
 # manual output
 manual/readuct_manual.aux
@@ -30,3 +32,8 @@ manual/readuct_manual.ilg
 
 # conan ci files
 conan/*
+
+# python
+*/__pycache__/
+*.pyc
+.pytest_cache/

CHANGELOG.rst

@@ -1,6 +1,15 @@
 Changelog
 =========
 
+Release 4.0.0
+-------------
+
+- Add 2nd Newton trajectory scan algorithm (NT2)
+- Update automatic TS mode picking to respect frequencies to some degree
+- Deprecate BondOrderTask and allow bond order calculation in SinglePointTask
+- Add option for spin propensity check
+- Add option to optimize periodic boundaries in geometry optimization
+
 Release 3.0.0
 -------------
 

CMakeLists.txt

@@ -4,7 +4,7 @@ cmake_minimum_required(VERSION 3.9)
 # tree must then provide a properly namespaced target with the same name as
 # your project.
 project(Readuct
-  VERSION 3.0.0
+  VERSION 4.0.0
   DESCRIPTION "This is the SCINE module Readuct."
 )
 
@@ -21,12 +21,18 @@ if(SCINE_BUILD_TESTS)
 endif()
 
 # Enable testing
-option(BUILD_SPARROW "Will download and build Sparrow (the Scine semi-empirical module)." OFF)
+option(BUILD_SPARROW "Will download and build Sparrow (the SCINE semi-empirical module)." OFF)
 if(BUILD_SPARROW)
   include(ImportSparrow)
   import_sparrow()
 endif()
 
+option(BUILD_XTB "Will download and build SCINE Xtb (the SCINE Wrapper around the xtb API by the Grimme group)." OFF)
+if(BUILD_XTB)
+  include(ImportXtb)
+  import_xtb()
+endif()
+
 option(SCINE_USE_MKL "Use the optimized MKL library for linear algebra operations of Eigen" OFF)
 
 # Subdirectories

README.rst

@@ -18,7 +18,7 @@ SCINE ReaDuct is a command-line tool that allows you to carry out
 
 For these calculations, it relies on a backend program to provide the necessary
 quantum chemical properties (such as nuclear gradients). Currently, SCINE Sparrow,
-SCINE XTB, CP2K, Gaussian, ORCA, and Turbomole are supported as backend programs.
+XTB, CP2K, Gaussian, ORCA, and Turbomole are supported as backend programs.
 
 License and Copyright Information
 ---------------------------------
@@ -33,8 +33,8 @@ Installation and Usage
 For instructions on how to install and use ReaDuct as well as for a detailed
 documentation of the entire functionality of ReaDuct, please consult the user
 manual found in the ``manual`` directory in in the repository.
-Alternatively the manual can also be found on the official GitHub website,
-on the SCINE website and in the hosted documentation.
+Alternatively the manual can also be found on the official GitHub website
+and on the SCINE website.
 
 How to Cite
 -----------

conanfile.py

@@ -3,7 +3,7 @@
 
 class ScineReaductConan(ScineConan):
     name = "scine_readuct"
-    version = "3.0.0"
+    version = "4.0.0"
     url = "https://github.com/qcscine/readuct"
     description = """
 SCINE ReaDuct is a command-line tool that allows to carry out:
@@ -18,7 +18,7 @@ class ScineReaductConan(ScineConan):
 - Newton trajectory scans searching for transition state guesses.
 For these calculations, it relies on a backend program to provide the necessary
 quantum chemical properties (such as nuclear gradients). Currently, SCINE Sparrow
-SCINE XTB, CP2K, Gaussian, ORCA, and Turbomole are supported as backend programs."""
+XTB, CP2K, Gaussian, ORCA, and Turbomole are supported as backend programs."""
     options = {
         "shared": [True, False],
         "python": [True, False],
@@ -38,7 +38,7 @@ class ScineReaductConan(ScineConan):
         "dev/cmake/*", "src/*", "CMakeLists.txt", "README.rst",
         "LICENSE.txt", "dev/conan/hook.cmake", "dev/conan/glue/*"
     ]
-    requires = ["scine_utilities/[~=4.0.0]",
+    requires = ["scine_utilities/[=5.0.0]",
                 "boost/[>1.65.0]",
                 "yaml-cpp/0.6.3"]
     cmake_name = "Readuct"
@@ -55,6 +55,6 @@ def configure(self):
 
     def build_requirements(self):
         if self.options.tests:
-            self.build_requires("scine_sparrow/[~=3.0.0]")
+            self.build_requires("scine_sparrow/[=3.0.1]")
 
         super().build_requirements()

src/Readuct/App/Tasks/AfirOptimizationTask.h

@@ -48,10 +48,12 @@ class AfirOptimizationTask : public Task {
     warningIfMultipleInputsGiven();
     warningIfMultipleOutputsGiven();
 
+    bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true);
     std::shared_ptr<Core::Calculator> calc;
     if (!testRunOnly) { // leave out in case of task chaining --> attention calc is NULL
       // Note: _input is guaranteed not to be empty by Task constructor
       calc = copyCalculator(systems, _input.front(), name());
+      Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator);
 
       // Check system size
       if (calc->getStructure()->size() == 1) {
@@ -66,11 +68,13 @@ class AfirOptimizationTask : public Task {
     if (optimizertype == "LBFGS") {
       auto tmp = std::make_shared<Utils::AfirOptimizer<Utils::Lbfgs>>(*calc);
       tmp->optimizer.useTrustRadius = true;
+      tmp->optimizer.trustRadius = 0.1;
       optimizer = std::move(tmp);
     }
     else if (optimizertype == "BFGS") {
       auto tmp = std::make_shared<Utils::AfirOptimizer<Utils::Bfgs>>(*calc);
       tmp->optimizer.useTrustRadius = true;
+      tmp->optimizer.trustRadius = 0.1;
       optimizer = std::move(tmp);
     }
     else if (optimizertype == "SD" || optimizertype == "STEEPESTDESCENT") {

src/Readuct/App/Tasks/BSplineInterpolationTask.h

@@ -63,6 +63,7 @@ class BSplineInterpolationTask : public Task {
     tangentFileName_ = taskSettings.extract("tangent_file", tangentFileName_);
     coordinateThresholdForMaximumExtraction_ =
         taskSettings.extract("extract_threshold", coordinateThresholdForMaximumExtraction_);
+    bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true);
 
     // If no errors encountered until here, the basic settings should be alright
     if (testRunOnly) {
@@ -72,6 +73,8 @@ class BSplineInterpolationTask : public Task {
     // Note: _input is guaranteed not to be empty by Task constructor
     auto calc = copyCalculator(systems, _input.front(), name());
     auto secondCalculator = systems.at(_input.back());
+    Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator);
+    Utils::CalculationRoutines::setLog(*secondCalculator, true, true, !silentCalculator);
     if (calc->settings() != secondCalculator->settings()) {
       _logger->warning
           << "  Warning: The given systems have different settings. Only taking first and ignoring second.\n";
@@ -296,7 +299,7 @@ class BSplineInterpolationTask : public Task {
                              "this optimization Task");
     }
 
-    const int cycles = optimizer->optimize();
+    const int cycles = optimizer->optimize(*_logger);
     const int maxiter = settings.getInt("convergence_max_iterations");
 
     bool converged = cycles < maxiter;

src/Readuct/App/Tasks/BondOrderTask.h

@@ -11,7 +11,7 @@
 #include "Tasks/Task.h"
 /* Scine */
 #include <Core/Interfaces/Calculator.h>
-#include <Utils/CalculatorBasics/Results.h>
+#include <Utils/CalculatorBasics.h>
 /* External */
 #include "boost/exception/diagnostic_information.hpp"
 #include <boost/filesystem.hpp>
@@ -27,6 +27,9 @@ class BondOrderTask : public Task {
  public:
   /**
    * @brief Construct a new BondOrderTask.
+   *
+   * @deprecated Please use the SinglePointTask with the `require_bond_orders` setting.
+   *
    * @param input  The input system names for the task.
    * @param output The output system names for the task.
    * @param logger The logger to/through which all text output will be handled.
@@ -45,6 +48,7 @@ class BondOrderTask : public Task {
 
     // Read and delete special settings
     bool stopOnError = stopOnErrorExtraction(taskSettings);
+    bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true);
     if (!taskSettings.empty()) {
       throw std::logic_error(falseTaskSettingsErrorMessage(name()));
     }
@@ -55,6 +59,7 @@ class BondOrderTask : public Task {
 
     // Note: _input is guaranteed not to be empty by Task constructor
     auto calc = copyCalculator(systems, _input.front(), name());
+    Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator);
 
     // Calculate bond orders and energy if not present in the results yet
     if (!calc->results().has<Utils::Property::BondOrderMatrix>()) {

src/Readuct/App/Tasks/GeometryOptimizationTask.h

@@ -13,6 +13,7 @@
 #include <Utils/GeometryOptimization/CoordinateSystem.h>
 #include <Utils/GeometryOptimization/GeometryOptimization.h>
 #include <Utils/GeometryOptimization/GeometryOptimizer.h>
+#include <Utils/GeometryOptimization/UnitCellGeometryOptimizer.h>
 #include <Utils/IO/ChemicalFileFormats/XyzStreamHandler.h>
 #include <Utils/Optimizer/GradientBased/Bfgs.h>
 #include <Utils/Optimizer/GradientBased/Lbfgs.h>
@@ -48,41 +49,19 @@ class GeometryOptimizationTask : public Task {
   bool run(SystemsMap& systems, Utils::UniversalSettings::ValueCollection taskSettings, bool testRunOnly = false) const final {
     warningIfMultipleInputsGiven();
     warningIfMultipleOutputsGiven();
+    bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true);
     // Get/Copy Calculator
     std::shared_ptr<Core::Calculator> calc;
     if (!testRunOnly) { // leave out in case of task chaining --> attention calc is NULL
       // Note: _input is guaranteed not to be empty by Task constructor
       calc = copyCalculator(systems, _input.front(), name());
+      Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator);
     }
 
     // Generate optimizer
     auto optimizertype = taskSettings.extract("optimizer", std::string{"BFGS"});
-    std::transform(optimizertype.begin(), optimizertype.end(), optimizertype.begin(), ::toupper);
-    std::shared_ptr<Utils::GeometryOptimizerBase> optimizer;
-    if (optimizertype == "LBFGS") {
-      auto tmp = std::make_shared<Utils::GeometryOptimizer<Utils::Lbfgs>>(*calc);
-      // Default convergence options
-      optimizer = std::move(tmp);
-    }
-    else if (optimizertype == "BFGS") {
-      auto tmp = std::make_shared<Utils::GeometryOptimizer<Utils::Bfgs>>(*calc);
-      // Default convergence options
-      optimizer = std::move(tmp);
-    }
-    else if (optimizertype == "SD" || optimizertype == "STEEPESTDESCENT") {
-      auto tmp = std::make_shared<Utils::GeometryOptimizer<Utils::SteepestDescent>>(*calc);
-      // Default convergence options
-      optimizer = std::move(tmp);
-    }
-    else if (optimizertype == "NR" || optimizertype == "NEWTONRAPHSON") {
-      auto tmp = std::make_shared<Utils::GeometryOptimizer<Utils::NewtonRaphson>>(*calc);
-      // Default convergence options
-      optimizer = std::move(tmp);
-    }
-    else {
-      throw std::runtime_error(
-          "Unknown Optimizer requested for a geometry optimization, available are: SD, NR, BFGS and LBFGS!");
-    }
+    auto unitcelloptimizertype = taskSettings.extract("unitcelloptimizer", std::string{""});
+    auto optimizer = constructOptimizer(*calc, optimizertype, unitcelloptimizertype);
 
     // Read and delete special settings
     bool stopOnError = stopOnErrorExtraction(taskSettings);
@@ -186,6 +165,86 @@ class GeometryOptimizationTask : public Task {
 
     return cycles < maxiter;
   }
+
+  inline std::shared_ptr<Utils::GeometryOptimizerBase> constructOptimizer(Core::Calculator& calc, std::string type,
+                                                                          std::string cellType) const {
+    std::transform(type.begin(), type.end(), type.begin(), ::toupper);
+    std::transform(cellType.begin(), cellType.end(), cellType.begin(), ::toupper);
+    if (type == "LBFGS") {
+      if (cellType.empty()) {
+        return std::make_shared<Utils::GeometryOptimizer<Utils::Lbfgs>>(calc);
+      }
+      else {
+        if (cellType == "LBFGS") {
+          return std::make_shared<Utils::UnitCellGeometryOptimizer<Utils::Lbfgs, Utils::Lbfgs>>(calc);
+        }
+        else if (cellType == "BFGS") {
+          return std::make_shared<Utils::UnitCellGeometryOptimizer<Utils::Lbfgs, Utils::Bfgs>>(calc);
+        }
+        else if (cellType == "SD" || cellType == "STEEPESTDESCENT") {
+          return std::make_shared<Utils::UnitCellGeometryOptimizer<Utils::Lbfgs, Utils::SteepestDescent>>(calc);
+        }
+        throw std::runtime_error(
+            "Unknown CellOptimizer requested for a geometry optimization, available are: SD, BFGS and LBFGS!");
+      }
+    }
+    else if (type == "BFGS") {
+      if (cellType.empty()) {
+        return std::make_shared<Utils::GeometryOptimizer<Utils::Bfgs>>(calc);
+      }
+      else {
+        if (cellType == "LBFGS") {
+          return std::make_shared<Utils::UnitCellGeometryOptimizer<Utils::Bfgs, Utils::Lbfgs>>(calc);
+        }
+        else if (cellType == "BFGS") {
+          return std::make_shared<Utils::UnitCellGeometryOptimizer<Utils::Bfgs, Utils::Bfgs>>(calc);
+        }
+        else if (cellType == "SD" || cellType == "STEEPESTDESCENT") {
+          return std::make_shared<Utils::UnitCellGeometryOptimizer<Utils::Bfgs, Utils::SteepestDescent>>(calc);
+        }
+        throw std::runtime_error(
+            "Unknown CellOptimizer requested for a geometry optimization, available are: SD, BFGS and LBFGS!");
+      }
+    }
+    else if (type == "SD" || type == "STEEPESTDESCENT") {
+      if (cellType.empty()) {
+        return std::make_shared<Utils::GeometryOptimizer<Utils::SteepestDescent>>(calc);
+      }
+      else {
+        if (cellType == "LBFGS") {
+          return std::make_shared<Utils::UnitCellGeometryOptimizer<Utils::SteepestDescent, Utils::Lbfgs>>(calc);
+        }
+        else if (cellType == "BFGS") {
+          return std::make_shared<Utils::UnitCellGeometryOptimizer<Utils::SteepestDescent, Utils::Bfgs>>(calc);
+        }
+        else if (cellType == "SD" || cellType == "STEEPESTDESCENT") {
+          return std::make_shared<Utils::UnitCellGeometryOptimizer<Utils::SteepestDescent, Utils::SteepestDescent>>(calc);
+        }
+        throw std::runtime_error(
+            "Unknown CellOptimizer requested for a geometry optimization, available are: SD, BFGS and LBFGS!");
+      }
+    }
+    else if (type == "NR" || type == "NEWTONRAPHSON") {
+      if (cellType.empty()) {
+        return std::make_shared<Utils::GeometryOptimizer<Utils::NewtonRaphson>>(calc);
+      }
+      else {
+        if (cellType == "LBFGS") {
+          return std::make_shared<Utils::UnitCellGeometryOptimizer<Utils::NewtonRaphson, Utils::Lbfgs>>(calc);
+        }
+        else if (cellType == "BFGS") {
+          return std::make_shared<Utils::UnitCellGeometryOptimizer<Utils::NewtonRaphson, Utils::Bfgs>>(calc);
+        }
+        else if (cellType == "SD" || cellType == "STEEPESTDESCENT") {
+          return std::make_shared<Utils::UnitCellGeometryOptimizer<Utils::NewtonRaphson, Utils::SteepestDescent>>(calc);
+        }
+        throw std::runtime_error(
+            "Unknown CellOptimizer requested for a geometry optimization, available are: SD, BFGS and LBFGS!");
+      }
+    }
+    throw std::runtime_error(
+        "Unknown Optimizer requested for a geometry optimization, available are: SD, NR, BFGS and LBFGS!");
+  };
 };
 
 } // namespace Readuct

src/Readuct/App/Tasks/HessianTask.h

@@ -50,6 +50,7 @@ class HessianTask : public Task {
 
     // Read and delete special settings
     bool stopOnError = stopOnErrorExtraction(taskSettings);
+    bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true);
     if (!taskSettings.empty()) {
       throw std::logic_error(falseTaskSettingsErrorMessage(name()));
     }
@@ -61,6 +62,7 @@ class HessianTask : public Task {
 
     // Note: _input is guaranteed not to be empty by Task constructor
     auto calc = copyCalculator(systems, _input.front(), name());
+    Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator);
 
     // Check system size
     if (calc->getStructure()->size() == 1) {

src/Readuct/App/Tasks/IrcTask.h

@@ -53,10 +53,12 @@ class IrcTask : public Task {
           << "  Warning: To store IRC results in a new location two output systems have to be specified.\n";
     }
 
+    bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true);
     std::shared_ptr<Core::Calculator> calc;
     if (!testRunOnly) { // leave out in case of task chaining --> attention calc is NULL
       // Note: _input is guaranteed not to be empty by Task constructor
       calc = copyCalculator(systems, _input.front(), name());
+      Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator);
 
       // Check system size
       if (calc->getStructure()->size() == 1) {
@@ -216,6 +218,9 @@ class IrcTask : public Task {
     // Add observer
     oldEnergy = 0.0;
     oldParams.resize(0);
+    // Reset optimizer
+    optimizer->setSettings(settings);
+    optimizer->reset();
     // Trajectory stream
     boost::filesystem::path dirB(((_output.size() > 1) ? _output[1] : _input[0]));
     boost::filesystem::create_directory(dirB);
@@ -241,6 +246,7 @@ class IrcTask : public Task {
 
     // Run optimization
     structure = systems.at(_input[0])->getStructure();
+    cycles = 0;
     try {
       cycles = optimizer->optimize(*structure, *_logger, ircModeVector, false);
     }