Skip to content

Commit

Permalink
Release 6.0.0
Browse files Browse the repository at this point in the history
  • Loading branch information
reiher-research-group committed Aug 26, 2024
1 parent 0bfa853 commit 69dbed5
Show file tree
Hide file tree
Showing 33 changed files with 874 additions and 234 deletions.
12 changes: 12 additions & 0 deletions CHANGELOG.rst
Original file line number Diff line number Diff line change
@@ -1,6 +1,18 @@
Changelog
=========

Release 6.0.0
-------------

- Add test for QM/MM transition state optimization
- Improve automatic mode selection by introducing a mode score considering a weighted sum of
the contributions of the relevant atoms and the wavenumber of the mode.
- Add the option to export the selected mode of a transition state optimization.
- Enable thermochemistry calculations for single atoms
- The one-electron integrals may be written to file through a dedicated integral evaluation task
- Separated energy and gradient contributions in QM/MM can be requested with the keyword "require_partial_energies",
and "require_partial_gradients" in the single-point task.

Release 5.1.0
-------------

Expand Down
23 changes: 16 additions & 7 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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 5.1.0
VERSION 6.0.0
DESCRIPTION "This is the SCINE module ReaDuct."
)

Expand All @@ -15,12 +15,9 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_CURRENT_SOURCE_DIR}/dev/cmake
include(ComponentSetup)
scine_setup_component()

# Enable testing
if(SCINE_BUILD_TESTS)
enable_testing()
endif()

# Enable testing
# Don't build tests of dependencies
set(_build_tests ${SCINE_BUILD_TESTS})
set(SCINE_BUILD_TESTS OFF)
option(BUILD_SPARROW "Will download and build Sparrow (the SCINE semi-empirical module)." OFF)
if(BUILD_SPARROW)
include(ImportSparrow)
Expand All @@ -33,6 +30,18 @@ if(BUILD_XTB)
import_xtb()
endif()

option(BUILD_SWOOSE "Will download and build Swoose (the SCINE QM/MM module)." OFF)
if(BUILD_SWOOSE)
include(ImportSwoose)
import_swoose()
endif()
set(SCINE_BUILD_TESTS ${_build_tests})

# Enable testing
if(SCINE_BUILD_TESTS)
enable_testing()
endif()

option(SCINE_USE_MKL "Use the optimized MKL library for linear algebra operations of Eigen" OFF)

# Subdirectories
Expand Down
7 changes: 7 additions & 0 deletions README.rst
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,13 @@ In addition, we kindly request you to cite the following article when using ReaD
A. C. Vaucher, M. Reiher, "Minimum Energy Paths and Transition States by Curve
Optimization", *J. Chem. Theory Comput.*, **2018**, *16*, 3091.

Furthermore, when publishing results obtained with any SCINE module, please cite the following paper:

T. Weymuth, J. P. Unsleber, P. L. Türtscher, M. Steiner, J.-G. Sobez, C. H. Müller, M. Mörchen,
V. Klasovita, S. A. Grimmel, M. Eckhoff, K.-S. Csizi, F. Bosia, M. Bensberg, M. Reiher,
"SCINE—Software for chemical interaction networks", *J. Chem. Phys.*, **2024**, *160*, 222501
(DOI `10.1063/5.0206974 <https://doi.org/10.1063/5.0206974>`_).

Support and Contact
-------------------

Expand Down
6 changes: 3 additions & 3 deletions conanfile.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@

class ScineReaductConan(ScineConan):
name = "scine_readuct"
version = "5.1.0"
version = "6.0.0"
url = "https://github.com/qcscine/readuct"
description = """
SCINE ReaDuct is a command-line tool that allows to carry out:
Expand Down Expand Up @@ -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/9.0.0",
requires = ["scine_utilities/[=10.0.0]",
"boost/[>1.65.0]",
"yaml-cpp/0.6.3"]
cmake_name = "Readuct"
Expand All @@ -55,6 +55,6 @@ def configure(self):

def build_requirements(self):
if self.options.tests:
self.build_requires("scine_sparrow/5.0.0")
self.build_requires("scine_sparrow/[=5.1.0]")

super().build_requirements()
279 changes: 148 additions & 131 deletions manual/readuct_manual.tex

Large diffs are not rendered by default.

24 changes: 19 additions & 5 deletions src/Readuct/App/Tasks/HessianTask.h
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
/* Scine */
#include <Core/Interfaces/Calculator.h>
#include <Utils/CalculatorBasics/Results.h>
#include <Utils/DataStructures/PartialHessian.h>
#include <Utils/GeometricDerivatives/NormalModeAnalysis.h>
#include <Utils/GeometricDerivatives/NormalModesContainer.h>
#include <Utils/IO/FormattedString.h>
Expand Down Expand Up @@ -67,21 +68,26 @@ class HessianTask : public Task {

// Note: _input is guaranteed not to be empty by Task constructor
auto calc = copyCalculator(systems, _input.front(), name());
const auto previousResults = calc->results();
auto previousResults = calc->results();
Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator);

// Check system size
// ToDo: Remove as soon as all calculators are able to handle single atom systems.
const bool singleAtom = calc->getStructure()->size() == 1;
if (calc->getStructure()->size() == 1) {
throw std::runtime_error("Cannot perform Hessian task for monoatomic systems.");
const auto hessian = Eigen::MatrixXd::Zero(3, 3);
calc->results().set<Utils::Property::Hessian>(hessian);
calc->results().set<Utils::Property::PartialHessian>(Utils::PartialHessian(hessian, {0}));
// Some calculators may remove the Hessian if only the energy is calculated.
previousResults.set<Utils::Property::Hessian>(hessian);
previousResults.set<Utils::Property::PartialHessian>(Utils::PartialHessian(hessian, {0}));
}
// Get/calculate Hessian
if (!(calc->results().has<Utils::Property::Hessian>() || calc->results().has<Utils::Property::PartialHessian>()) ||
!calc->results().has<Utils::Property::Energy>()) {
// has neither type of Hessian or has no energy
Utils::PropertyList requiredProperties = Utils::Property::Energy;
if (!calc->results().has<Utils::Property::Thermochemistry>() &&
calc->possibleProperties().containsSubSet(Utils::Property::Thermochemistry)) {
calc->possibleProperties().containsSubSet(Utils::Property::Thermochemistry) && !singleAtom) {
requiredProperties.addProperty(Utils::Property::Thermochemistry);
// many calculators cannot handle that only Thermochemistry is requested
// because they delete their results and then assume that Hessian was also calculated
Expand Down Expand Up @@ -156,7 +162,7 @@ class HessianTask : public Task {

// Setup output directory if needed
boost::filesystem::path dir(outputSystem);
if (wavenumbers[0] < 0.0) {
if (!singleAtom && wavenumbers[0] < 0.0) {
boost::filesystem::create_directory(dir);
}
// Print Imaginary modes
Expand All @@ -171,6 +177,14 @@ class HessianTask : public Task {
// Store result
systems[outputSystem] = calc;

if (singleAtom) {
const auto hessian = calc->results().get<Utils::Property::Hessian>();
auto thermoCalc = Scine::Utils::ThermochemistryCalculator(
hessian, *calc->getStructure(), calc->settings().getInt(Utils::SettingsNames::spinMultiplicity),
calc->results().get<Utils::Property::Energy>());
calc->results().set<Utils::Property::Thermochemistry>(thermoCalc.calculate());
}

// Print thermochemistry results
if (calc->results().has<Utils::Property::Thermochemistry>()) {
auto thermo = calc->results().get<Utils::Property::Thermochemistry>();
Expand Down
144 changes: 144 additions & 0 deletions src/Readuct/App/Tasks/IntegralTask.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,144 @@
/**
* @file
* @copyright This code is licensed under the 3-clause BSD license.\n
* Copyright ETH Zurich, Department of Chemistry and Applied Biosciences, Reiher Group.\n
* See LICENSE.txt for details.
*/

#ifndef READUCT_INTEGRALSTASK_H
#define READUCT_INTEGRALSTASK_H

/* Readuct */
#include "Tasks/Task.h"
/* Scine */
#include <Core/Interfaces/Calculator.h>
#include <Utils/CalculatorBasics.h>
#include <Utils/Settings.h>
#include <Utils/UniversalSettings/SettingsNames.h>
/* External */
#include "boost/exception/diagnostic_information.hpp"
#include <boost/filesystem.hpp>
/* std c++ */
#include <iomanip>
#include <string>
#include <vector>

namespace Scine::Readuct {

/**
* @class IntegralTask IntegralTask.h
* @brief This tasks writes integrals, such as the one-electron integrals, to file.
*/
class IntegralTask : public Task {
public:
/**
* @brief Construct a new IntegralTask.
* @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.
*/
IntegralTask(std::vector<std::string> input, std::vector<std::string> output, std::shared_ptr<Core::Log> logger = nullptr)
: Task(std::move(input), std::move(output), std::move(logger)) {
}
/**
* @brief Getter for the task's name.
* @return The task's name.
*/
std::string name() const override {
return "Integral Calculation";
}
/**
* @brief Run the task.
* @param systems The input systems.
* @param taskSettings The task settings.
* @param testRunOnly If true, the main part of the task is not executed. Only used for testing.
* @param observers Optional observers. This task does not support observers at the moment.
* @return Returns true if the task terminated successfully. False, otherwise.
*/
bool run(SystemsMap& systems, Utils::UniversalSettings::ValueCollection taskSettings, bool testRunOnly = false,
std::vector<std::function<void(const int&, const Utils::AtomCollection&, const Utils::Results&, const std::string&)>>
observers = {}) const final {
// Read and delete special settings
const bool stopOnError = stopOnErrorExtraction(taskSettings);
const bool requireHCore = taskSettings.extract("require_one_electron_integrals", true);
const bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true);
if (!taskSettings.empty()) {
std::string keyListing = "\n";
auto keys = taskSettings.getKeys();
for (const auto& key : keys) {
keyListing += "'" + key + "'\n";
}
throw std::logic_error("Specified one or more task settings that are not available for this task:" + keyListing);
}
if (observers.size() > 0) {
throw std::logic_error("IntegralTask does not feature algorithm accepting observers, yet one was given");
}
if (testRunOnly) {
return true; // leave out rest in case of task chaining
}

// Note: _input is guaranteed not to be empty by Task constructor
auto calc = copyCalculator(systems, _input.front(), name());
const auto previousResults = calc->results();
Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator);
// Check for available properties
const bool hCoreAvailable = calc->possibleProperties().containsSubSet(Utils::Property::OneElectronMatrix);

if (requireHCore && !hCoreAvailable) {
throw std::logic_error(
"One-electron integrals required, but chosen calculator does not provide them.\n"
"If you do not need one-electron integrals, set 'require_one_electron_integrals' to 'false'"
" in the task settings");
}
Utils::PropertyList requiredProperties;
if (requireHCore) {
requiredProperties.addProperty(Utils::Property::OneElectronMatrix);
}

try {
calc->setRequiredProperties(requiredProperties);
calc->calculate(name());
if (!calc->results().get<Utils::Property::SuccessfulCalculation>()) {
throw std::runtime_error(name() + " was not successful");
}
}
catch (...) {
if (stopOnError) {
throw;
}
_logger->error
<< " " + name() + " was not successful with error:\n " + boost::current_exception_diagnostic_information()
<< Core::Log::endl;
calc->results() = previousResults + calc->results();
return false;
}

if (requireHCore) {
const std::string& outputSystem = (!_output.empty() ? _output.front() : _input.front());
const boost::filesystem::path dir(outputSystem);
boost::filesystem::create_directory(dir);
const boost::filesystem::path asciiFile(outputSystem + ".hcore.dat");
const std::string outPath = (dir / asciiFile).string();
_logger->output << "Writing core Hamiltonian integrals to file " << outPath << Core::Log::endl;
std::ofstream outFile(outPath, std::ofstream::out);
const auto hCoreMatrix = calc->results().template get<Utils::Property::OneElectronMatrix>();
outFile << "AO core Hamiltonian integrals.\n";
outFile << std::scientific << std::setprecision(12) << hCoreMatrix << "\n";
outFile.close();

// Store result
if (!_output.empty()) {
systems[_output[0]] = calc;
}
else {
systems[_input[0]] = calc;
}
}

return true;
}
};

} // namespace Scine::Readuct

#endif // READUCT_INTEGRALSTASK_H
51 changes: 46 additions & 5 deletions src/Readuct/App/Tasks/SinglePointTask.h
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,8 @@ class SinglePointTask : public Task {
bool requireGradients = taskSettings.extract("require_gradients", false);
bool requireStressTensor = taskSettings.extract("require_stress_tensor", false);
bool requireBondOrders = taskSettings.extract("require_bond_orders", false);
bool requirePartialEnergies = taskSettings.extract("require_partial_energies", false);
bool requirePartialGradients = taskSettings.extract("require_partial_gradients", false);
bool requireOrbitalEnergies = taskSettings.extract("orbital_energies", false);
bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true);
int spinPropensityCheck = taskSettings.extract("spin_propensity_check", 0);
Expand All @@ -84,11 +86,13 @@ class SinglePointTask : public Task {
const auto previousResults = calc->results();
Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator);
// Check for available properties
bool chargesAvailable = calc->possibleProperties().containsSubSet(Utils::Property::AtomicCharges);
bool gradientsAvailable = calc->possibleProperties().containsSubSet(Utils::Property::Gradients);
bool stressTensorAvailable = calc->possibleProperties().containsSubSet(Utils::Property::StressTensor);
bool orbitalEnergiesAvailable = calc->possibleProperties().containsSubSet(Utils::Property::OrbitalEnergies);
bool bondOrdersAvailable = calc->possibleProperties().containsSubSet(Utils::Property::BondOrderMatrix);
const bool chargesAvailable = calc->possibleProperties().containsSubSet(Utils::Property::AtomicCharges);
const bool gradientsAvailable = calc->possibleProperties().containsSubSet(Utils::Property::Gradients);
const bool stressTensorAvailable = calc->possibleProperties().containsSubSet(Utils::Property::StressTensor);
const bool orbitalEnergiesAvailable = calc->possibleProperties().containsSubSet(Utils::Property::OrbitalEnergies);
const bool bondOrdersAvailable = calc->possibleProperties().containsSubSet(Utils::Property::BondOrderMatrix);
const bool partialEnergiesAvailable = calc->possibleProperties().containsSubSet(Utils::Property::PartialEnergies);
const bool partialGradientsAvailable = calc->possibleProperties().containsSubSet(Utils::Property::PartialGradients);
if (requireCharges && !chargesAvailable) {
throw std::logic_error("Charges required, but chosen calculator does not provide them.\n"
"If you do not need charges, set 'require_charges' to 'false' in the task settings");
Expand All @@ -105,6 +109,12 @@ class SinglePointTask : public Task {
if (requireBondOrders && !bondOrdersAvailable) {
throw std::logic_error("Bond orders required, but chosen calculator does not provide them.");
}
if (requirePartialEnergies && !partialEnergiesAvailable) {
throw std::logic_error("Partial energies were required, but the chosen calculator does not provide them.");
}
if (requirePartialGradients && !partialGradientsAvailable) {
throw std::logic_error("Partial gradients were required, but the chosen calculator does not provide them.");
}
// Calculate energy, and possibly atomic charges and/or gradients
Utils::PropertyList requiredProperties = Utils::Property::Energy;
if (requireCharges) {
Expand All @@ -122,6 +132,12 @@ class SinglePointTask : public Task {
if (requireBondOrders) {
requiredProperties.addProperty(Utils::Property::BondOrderMatrix);
}
if (requirePartialEnergies) {
requiredProperties.addProperty(Utils::Property::PartialEnergies);
}
if (requirePartialGradients) {
requiredProperties.addProperty(Utils::Property::PartialGradients);
}

try {
calc->setRequiredProperties(requiredProperties);
Expand Down Expand Up @@ -263,6 +279,31 @@ class SinglePointTask : public Task {
}
}
}
if (requirePartialEnergies) {
auto partiaEnergies = calc->results().get<Utils::Property::PartialEnergies>();
cout << " Partial Energies:\n\n";
if (partiaEnergies.find("qm_energy") != partiaEnergies.end()) {
cout.printf(" The QM energy is: %+16.9f hartree\n", partiaEnergies["qm_energy"]);
}
if (partiaEnergies.find("mm_energy") != partiaEnergies.end()) {
cout.printf(" The MM energy is: %+16.9f hartree\n", partiaEnergies["mm_energy"]);
}
}
if (requirePartialGradients) {
auto partialGradients = calc->results().get<Utils::Property::PartialGradients>();
if (partialGradients.find("qm_gradients") != partialGradients.end()) {
const auto& qmGradients = partialGradients["qm_gradients"];
cout << " QM gradients (hartree / bohr):\n\n";
cout << [&qmGradients](std::ostream& os) { Utils::matrixPrettyPrint(os, qmGradients); };
cout << Core::Log::nl;
}
if (partialGradients.find("mm_gradients") != partialGradients.end()) {
const auto& mmGradients = partialGradients["mm_gradients"];
cout << " MM gradients (hartree / bohr):\n\n";
cout << [&mmGradients](std::ostream& os) { Utils::matrixPrettyPrint(os, mmGradients); };
cout << Core::Log::nl;
}
}

cout << Core::Log::nl << Core::Log::endl;
return true;
Expand Down
Loading

0 comments on commit 69dbed5

Please sign in to comment.