diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 026676c..f529556 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,6 +1,10 @@ Changelog ========= +Release 2.1.0 +------------- +- rework non-bonded interaction calculation to increase memory efficiency. + Release 2.0.0 ------------- - enable electrostatic embedding QM/MM with Turbomole and xtb as QM calculators diff --git a/CMakeLists.txt b/CMakeLists.txt index edd4867..e2a4785 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(Swoose - VERSION 2.0.0 + VERSION 2.1.0 DESCRIPTION "This is the SCINE module Swoose." ) diff --git a/README.rst b/README.rst index 998343c..bd40e18 100644 --- a/README.rst +++ b/README.rst @@ -55,8 +55,15 @@ K.-S. Csizi, M. Reiher, "Automated preparation of nanoscopic structures: Graph-b sequence analysis, mismatch detection, and pH-consistent protonation with uncertainty estimates" *J. Comput. Chem*, **2024**, *45*, 761-776. -K.-S. Csizi, M. Steiner, M. Reiher, "Quantum Magnifying Glass for Chemistry at the Nanoscale", -*chemRxiv", **2023**, DOI 10.26434/chemrxiv-2023-t10sc. +K.-S. Csizi, M. Steiner, M. Reiher, "Nanoscale chemical reaction exploration with a quantum magnifying glass", +*Nat. Commun.*, **2024**, 15, 5320. + +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 `_). Support and Contact ------------------- diff --git a/conanfile.py b/conanfile.py index b5daef1..eec437e 100644 --- a/conanfile.py +++ b/conanfile.py @@ -3,7 +3,7 @@ class ScineSwooseConan(ScineConan): name = "scine_swoose" - version = "2.0.0" + version = "2.1.0" url = "https://github.com/qcscine/swoose" description = """ Treat large molecular systems with self-parametrizing system-focused atomistic @@ -22,8 +22,8 @@ class ScineSwooseConan(ScineConan): exports_sources = ["dev/cmake/*", "src/*", "CMakeLists.txt", "README.rst", "LICENSE.txt", "dev/conan/hook.cmake", "dev/conan/glue/*"] requires = ["yaml-cpp/0.6.3", - "scine_utilities/[==9.0.0]", - "scine_molassembler/[==2.0.1]"] + "scine_utilities/10.0.0", + "scine_molassembler/3.0.0"] cmake_name = "Swoose" def requirements(self): @@ -32,7 +32,7 @@ def requirements(self): } if self.options.get_safe("database"): - self.requires("scine_database/[==1.3.0]") + self.requires("scine_database/1.4.0") if hasattr(super(), "requirements"): super().requirements() diff --git a/dev b/dev index 518ab3c..2712fd3 160000 --- a/dev +++ b/dev @@ -1 +1 @@ -Subproject commit 518ab3c7f8a0a724081fcd7ed518c669724bcd37 +Subproject commit 2712fd3204d703c777244b077066ad75c6bc0db1 diff --git a/manual/swoose_manual.tex b/manual/swoose_manual.tex index efe7b79..d84263a 100644 --- a/manual/swoose_manual.tex +++ b/manual/swoose_manual.tex @@ -18,8 +18,8 @@ %% % Book metadata -\title[SCINE Swoose manual]{User Manual \vskip 0.5em {\setlength{\parindent}{0pt} \Huge SCINE Swoose 2.0.0}} -\author[The SCINE Swoose Developers]{\noindent Christoph Brunken, Katja-Sophia Csizi, Miguel Steiner, Thomas Weymuth, and Markus Reiher} +\title[SCINE Swoose manual]{User Manual \vskip 0.5em {\setlength{\parindent}{0pt} \Huge SCINE Swoose 2.1.0}} +\author[The SCINE Swoose Developers]{\noindent Moritz Bensberg, Christoph Brunken, Katja-Sophia Csizi, Miguel Steiner, Thomas Weymuth, and Markus Reiher} \publisher{ETH Z\"urich} %% @@ -276,7 +276,7 @@ \section{Conan package} \begin{mdframed}[backgroundcolor=LightSteelBlue!25, linewidth=0pt] \begin{verbatim} -conan install scine_swoose/2.0.0@ +conan install scine_swoose/2.1.0@ \end{verbatim} \end{mdframed} @@ -286,7 +286,7 @@ \section{Conan package} \begin{mdframed}[backgroundcolor=LightSteelBlue!25, userdefinedwidth=12cm, linewidth=0pt] \begin{verbatim} -conan install -o scine_swoose:python=True scine_swoose/2.0.0@ +conan install -o scine_swoose:python=True scine_swoose/2.1.0@ \end{verbatim} \end{mdframed} @@ -302,7 +302,7 @@ \section{Conan package} \begin{mdframed}[backgroundcolor=LightSteelBlue!25, linewidth=0pt] \begin{verbatim} -conan install scine_swoose/2.0.0@ -g virtualrunenv +conan install scine_swoose/2.1.0@ -g virtualrunenv \end{verbatim} \end{mdframed} @@ -321,7 +321,7 @@ \section{Conan package} \begin{mdframed}[backgroundcolor=LightSteelBlue!25, linewidth=0pt] \begin{verbatim} [requires] -scine_swoose/[>=2.0.0]@ +scine_swoose/[>=2.1.0]@ scine_sparrow/[>=5.0.0]@ scine_xtb_wrapper/[>=3.0.0]@ diff --git a/src/Swoose/CMakeLists.txt b/src/Swoose/CMakeLists.txt index b3092df..a28d29e 100644 --- a/src/Swoose/CMakeLists.txt +++ b/src/Swoose/CMakeLists.txt @@ -162,6 +162,57 @@ if(SCINE_BUILD_PYTHON_BINDINGS) import_pybind11() set(PYBIND11_PYTHON_VERSION ${PYTHONVERSION}) + # Figure out which targets the python bindings are going to need copied + include(TargetLibName) + set(_py_targets_to_copy Swoose) + target_lib_type(Scine::UtilsOS _utils_libtype) + if(_utils_libtype STREQUAL "SHARED_LIBRARY") + list(APPEND _py_targets_to_copy Scine::UtilsOS) + endif() + unset(_utils_libtype) + + # Generate generator expressions for each targets and figure out filenames + # for the python setup file + set(swoose_PY_DEPS "") + foreach(target ${_py_targets_to_copy}) + list(APPEND _py_target_gen_exprs "\$") + target_lib_filename(${target} _target_filename) + string(APPEND swoose_PY_DEPS ", \"${_target_filename}\"") + endforeach() + message(STATUS "Targets to copy for Python bindings: ${_py_targets_to_copy}") + unset(_py_targets_to_copy) + + add_custom_command(TARGET Swoose POST_BUILD + COMMAND ${CMAKE_COMMAND} -E copy ${_py_target_gen_exprs} ${CMAKE_CURRENT_BINARY_DIR}/scine_swoose + COMMENT "Copying required shared libraries into Python package directory" + ) + unset(_py_target_gen_exprs) + + # Utils python dependency + include(FindPythonModule) + if(scine-utils-os_BINARY_DIR) + set(UTILS_PYTHONPATH ${scine-utils-os_BINARY_DIR}/src/Utils) + else() + find_python_module(scine_utilities) + if(PY_SCINE_UTILITIES) + set(UTILS_PYTHONPATH ${PY_SCINE_UTILITIES}) + else() + message(WARNING "Utilities Python module not found. Cannot test Python module or generate stubs") + endif() + endif() + + # Molassembler python dependency + if(molassembler_BINARY_DIR) + set(MASM_PYTHONPATH ${scine-utils-os_BINARY_DIR}/src/Utils) + else() + find_python_module(scine_molassembler) + if(PY_SCINE_MOLASSEMBLER) + set(MASM_PYTHONPATH ${PY_SCINE_MOLASSEMBLER}) + else() + message(WARNING "Molassembler Python module not found. Cannot test Python module or generate stubs") + endif() + endif() + # Python module pybind11_add_module(scine_swoose ${SWOOSE_PYTHON_FILES} @@ -238,15 +289,19 @@ if(SCINE_BUILD_PYTHON_BINDINGS) endif() # Typing stubs - include(FindPythonModule) find_python_module(pybind11_stubgen) - if(PY_PYBIND11_STUBGEN) + if(PY_PYBIND11_STUBGEN AND UTILS_PYTHONPATH AND MASM_PYTHONPATH) add_custom_command(TARGET scine_swoose POST_BUILD - COMMAND ${PYTHON_EXECUTABLE} -m pybind11_stubgen -o . --root-module-suffix \"\" --no-setup-py --bare-numpy-ndarray scine_swoose + COMMAND ${CMAKE_COMMAND} -E env PYTHONPATH=${UTILS_PYTHONPATH}:$ENV{PYTHONPATH} ${PYTHON_EXECUTABLE} -m pybind11_stubgen -o . --root-module-suffix \"\" --no-setup-py --bare-numpy-ndarray scine_swoose COMMENT "Generating python package typechecking stubs with pybind11-stubgen" - ) - else() + BYPRODUCTS ${CMAKE_CURRENT_BINARY_DIR}/scine_swoose/__init__.pyi + ) + elseif(UTILS_PYTHONPATH AND MASM_PYTHONPATH) message(STATUS "Not generating typechecking stubs for python package as pybind11-stubgen was not found") + elseif(UTILS_PYTHONPATH) + message(STATUS "Not generating typechecking stubs for Python package as SCINE Molassembler dependency was not found") + else() + message(STATUS "Not generating typechecking stubs for Python package as SCINE Utilities dependency was not found") endif() # Add setuptools file diff --git a/src/Swoose/Files.cmake b/src/Swoose/Files.cmake index 5a7cd01..e052091 100644 --- a/src/Swoose/Files.cmake +++ b/src/Swoose/Files.cmake @@ -9,6 +9,10 @@ set(SWOOSE_MODULE_FILES ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/GAFF/GaffCalculatorSettings.h ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/GAFF/GaffMolecularMechanicsCalculator.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/GAFF/GaffMolecularMechanicsCalculator.h + ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/InteractionExclusion.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/InteractionExclusion.h + ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/ScaledInteractions.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/ScaledInteractions.h ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/MMExceptions.h ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/MMParameters.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/MMParameters.h @@ -24,8 +28,6 @@ set(SWOOSE_MODULE_FILES ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/SFAM/SfamAtomTypeIdentifier.h ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/Electrostatic.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/Electrostatic.h - ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/ElectrostaticTerm.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/ElectrostaticTerm.h ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/ElectrostaticEvaluator.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/ElectrostaticEvaluator.h ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/Repulsion.cpp @@ -75,8 +77,6 @@ set(SWOOSE_MODULE_FILES ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/LennardJones.h ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/LennardJonesEvaluator.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/LennardJonesEvaluator.h - ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/LennardJonesTerm.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/LennardJonesTerm.h ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/HydrogenBond.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/HydrogenBond.h ${CMAKE_CURRENT_SOURCE_DIR}/Swoose/MolecularMechanics/Interactions/HydrogenBondTerm.cpp diff --git a/src/Swoose/Python/setup.py b/src/Swoose/Python/setup.py index f88d935..1ba4d4b 100644 --- a/src/Swoose/Python/setup.py +++ b/src/Swoose/Python/setup.py @@ -9,7 +9,7 @@ import os # Read README.rst for the long description -with open("README.rst", "r") as fh: +with open("README.rst", "r", encoding="utf-8") as fh: long_description = fh.read() diff --git a/src/Swoose/Swoose/MMParametrization/ParametrizationUtils/ReparametrizationHelper.h b/src/Swoose/Swoose/MMParametrization/ParametrizationUtils/ReparametrizationHelper.h index 634f719..ea85f1f 100644 --- a/src/Swoose/Swoose/MMParametrization/ParametrizationUtils/ReparametrizationHelper.h +++ b/src/Swoose/Swoose/MMParametrization/ParametrizationUtils/ReparametrizationHelper.h @@ -14,7 +14,7 @@ namespace Scine { namespace Core { -class Log; +struct Log; } // namespace Core namespace MMParametrization { diff --git a/src/Swoose/Swoose/MMParametrization/ParametrizationUtils/SuperfluousFragmentIdentifier.h b/src/Swoose/Swoose/MMParametrization/ParametrizationUtils/SuperfluousFragmentIdentifier.h index 5b34795..ab5489c 100644 --- a/src/Swoose/Swoose/MMParametrization/ParametrizationUtils/SuperfluousFragmentIdentifier.h +++ b/src/Swoose/Swoose/MMParametrization/ParametrizationUtils/SuperfluousFragmentIdentifier.h @@ -13,7 +13,7 @@ namespace Scine { namespace Core { -class Log; +struct Log; } // namespace Core namespace MMParametrization { diff --git a/src/Swoose/Swoose/MMParametrization/ParametrizationUtils/TitrationHelper.cpp b/src/Swoose/Swoose/MMParametrization/ParametrizationUtils/TitrationHelper.cpp index 1c1cd4a..ce923a3 100644 --- a/src/Swoose/Swoose/MMParametrization/ParametrizationUtils/TitrationHelper.cpp +++ b/src/Swoose/Swoose/MMParametrization/ParametrizationUtils/TitrationHelper.cpp @@ -92,6 +92,7 @@ double TitrationHelper::getPkaOfTrainingMolecule(std::string dataDirectory) { } double TitrationHelper::getEnergyOfDeprotonationForTrainingMolecule(std::string dataDirectory) { + (void)dataDirectory; return 0.0; } diff --git a/src/Swoose/Swoose/MMParametrization/ReferenceCalculationHelpers/DirectCalculationsHelper.h b/src/Swoose/Swoose/MMParametrization/ReferenceCalculationHelpers/DirectCalculationsHelper.h index ec51cbd..be8e707 100644 --- a/src/Swoose/Swoose/MMParametrization/ReferenceCalculationHelpers/DirectCalculationsHelper.h +++ b/src/Swoose/Swoose/MMParametrization/ReferenceCalculationHelpers/DirectCalculationsHelper.h @@ -13,7 +13,7 @@ namespace Scine { namespace Core { -class Log; +struct Log; } // namespace Core namespace Utils { diff --git a/src/Swoose/Swoose/MachineLearning/MolecularMachineLearningModel.cpp b/src/Swoose/Swoose/MachineLearning/MolecularMachineLearningModel.cpp index 9cb9b26..0669bde 100644 --- a/src/Swoose/Swoose/MachineLearning/MolecularMachineLearningModel.cpp +++ b/src/Swoose/Swoose/MachineLearning/MolecularMachineLearningModel.cpp @@ -191,9 +191,9 @@ Eigen::MatrixXd MolecularMachineLearningModel::getSingleForceFeatures(int atomIn } Eigen::MatrixXd MolecularMachineLearningModel::getSingleForceTargets(int atomIndex) { - int nDataPoints = forcesFeatures_.size(); + const size_t nDataPoints = forcesFeatures_.size(); Eigen::MatrixXd targets(nDataPoints, 3); - for (long unsigned int j = 0; j < nDataPoints; ++j) { + for (size_t j = 0; j < nDataPoints; ++j) { targets.row(j) = refForces_[j].row(atomIndex); } return targets; diff --git a/src/Swoose/Swoose/MolecularMechanics/AtomTypesHolder.h b/src/Swoose/Swoose/MolecularMechanics/AtomTypesHolder.h index 0b24fbd..e3902c3 100644 --- a/src/Swoose/Swoose/MolecularMechanics/AtomTypesHolder.h +++ b/src/Swoose/Swoose/MolecularMechanics/AtomTypesHolder.h @@ -28,7 +28,7 @@ class AtomTypesHolder { /** * @brief Getter for the atom type for an atom with a certain index. */ - std::string getAtomType(int index) const; + const std::string& getAtomType(unsigned int index) const; /** * @brief Returns the number of atom types stored in this object. @@ -42,7 +42,7 @@ class AtomTypesHolder { inline AtomTypesHolder::AtomTypesHolder(std::vector atomTypes) : atomTypes_(std::move(atomTypes)) { } -inline std::string AtomTypesHolder::getAtomType(int index) const { +inline const std::string& AtomTypesHolder::getAtomType(unsigned int index) const { return atomTypes_.at(index); } diff --git a/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffMolecularMechanicsCalculator.cpp b/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffMolecularMechanicsCalculator.cpp index 16d7c87..6ef2059 100644 --- a/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffMolecularMechanicsCalculator.cpp +++ b/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffMolecularMechanicsCalculator.cpp @@ -110,7 +110,6 @@ const Utils::Results& GaffMolecularMechanicsCalculator::calculate(std::string de const Utils::Results& GaffMolecularMechanicsCalculator::calculateImpl(std::string description) { int nAtoms = structure_.size(); double energy = 0.0; - // Initialize the atomic derivatives container Utils::AtomicSecondDerivativeCollection derivativesForBondedInteractions(nAtoms); derivativesForBondedInteractions.setZero(); @@ -119,9 +118,10 @@ const Utils::Results& GaffMolecularMechanicsCalculator::calculateImpl(std::strin int nAtomsInitialization = 0; if (!onlyCalculateBondedContribution_) nAtomsInitialization = nAtoms; - Utils::FullSecondDerivativeCollection fullDerivatives(nAtomsInitialization); - fullDerivatives.setZero(); + const unsigned int derivativeOrder = (requiredProperties_.containsSubSet(Utils::Property::Hessian)) ? 2 : 1; + Utils::DerivativeCollection fullDerivatives(nAtomsInitialization, derivativeOrder); + fullDerivatives.setZero(); double energyBonds = bondsEvaluator_->evaluate(derivativesForBondedInteractions); energy += energyBonds; double energyAngles = anglesEvaluator_->evaluate(derivativesForBondedInteractions); @@ -130,7 +130,6 @@ const Utils::Results& GaffMolecularMechanicsCalculator::calculateImpl(std::strin energy += energyDihedrals; double energyImproperDihedrals = improperDihedralsEvaluator_->evaluate(derivativesForBondedInteractions); energy += energyImproperDihedrals; - double energyLennardJones = 0.0; double energyElectro = 0.0; if (!onlyCalculateBondedContribution_) { @@ -218,7 +217,6 @@ void GaffMolecularMechanicsCalculator::generatePotentialTerms(const std::string& GaffAtomTypeIdentifier atomTypeGenerator(structure_.getElements().size(), structure_.getElements(), listsOfNeighbors_, atomTypesFile_); auto atomTypes = atomTypeGenerator.getAtomTypes(); - if (!parametersHaveBeenSetInternally_ || parameterFilePathHasBeenChanged_) { if (parameterPath.empty()) { GaffParameterDefaultsProvider parameterProvider; @@ -234,7 +232,8 @@ void GaffMolecularMechanicsCalculator::generatePotentialTerms(const std::string& parametersHaveBeenSetInternally_ = true; parameterFilePathHasBeenChanged_ = false; } - + lennardJonesEvaluator_->setAtomTypesHolder(std::make_shared(atomTypes)); + lennardJonesEvaluator_->setParameters(std::make_shared(parameters_)); generatePotentialTerms(parameters_, topology, atomTypes); } @@ -256,15 +255,24 @@ void GaffMolecularMechanicsCalculator::generatePotentialTerms(const GaffParamete // Generate the potential terms GaffPotentialTermsGenerator generator(structure_.size(), atomTypes, topology, parameters, structure_.getPositions(), - nonCovalentCutoffRadius_, this->getLog()); + this->getLog()); bondsEvaluator_->setBondTerms(generator.getBondedTerms()); anglesEvaluator_->setAngleTerms(generator.getAngleTerms()); dihedralsEvaluator_->setDihedralTerms(generator.getDihedralTerms()); improperDihedralsEvaluator_->setDihedralTerms(generator.getImproperDihedralTerms()); if (!onlyCalculateBondedContribution_) { - lennardJonesEvaluator_->setLennardJonesTerms(generator.getLennardJonesTerms(applyCutoffDuringInitialization_)); - electrostaticEvaluator_->setElectrostaticTerms(generator.getElectrostaticTerms(applyCutoffDuringInitialization_)); + lennardJonesEvaluator_->setCutOffRadius(std::make_shared(nonCovalentCutoffRadius_)); + lennardJonesEvaluator_->resetExclusions(structure_.size()); + lennardJonesEvaluator_->addExclusions(topology); + lennardJonesEvaluator_->resetScaledInteractions(structure_.size()); + lennardJonesEvaluator_->addScaledInteractionPairs(topology); + + electrostaticEvaluator_->setCutOffRadius(std::make_shared(nonCovalentCutoffRadius_)); + electrostaticEvaluator_->resetExclusions(structure_.size()); + electrostaticEvaluator_->addExclusions(topology); + electrostaticEvaluator_->resetScaledInteractions(structure_.size()); + electrostaticEvaluator_->addScaledInteractionPairs(topology); } } @@ -289,7 +297,6 @@ void GaffMolecularMechanicsCalculator::initialize() { SwooseUtilities::TopologyUtils::generateListsOfNeighborsFromBondOrderMatrix(structure_.size(), bondOrders, 0.5); } } - generatePotentialTerms(parameterFilePath_); } diff --git a/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffMolecularMechanicsCalculator.h b/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffMolecularMechanicsCalculator.h index c451769..9122439 100644 --- a/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffMolecularMechanicsCalculator.h +++ b/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffMolecularMechanicsCalculator.h @@ -27,10 +27,10 @@ class InteractionTermEliminator; namespace MolecularMechanics { class GaffParameters; class AtomTypesHolder; -class BondType; -class AngleType; -class DihedralType; -class ImproperDihedralType; +struct BondType; +struct AngleType; +struct DihedralType; +struct ImproperDihedralType; /** * @class GaffMolecularMechanicsCalculator GaffMolecularMechanicsCalculator.h * @brief Calculator for the GAFF Molecular Mechanics method. diff --git a/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffParameters.cpp b/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffParameters.cpp index 0846b34..3ca2ed8 100644 --- a/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffParameters.cpp +++ b/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffParameters.cpp @@ -70,7 +70,7 @@ std::vector GaffParameters::getMMImproperDihedrals(std::string central return improperDihedralsForGivenAtomTypes; } -LennardJones GaffParameters::getMMLennardJones(std::string t1, std::string t2, double scalingFactor) const { +LennardJones GaffParameters::getMMLennardJones(const std::string& t1, const std::string& t2, double scalingFactor) const { auto vdwAtom1 = lennardJonesPairs_.find(t1); auto vdwAtom2 = lennardJonesPairs_.find(t2); diff --git a/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffParameters.h b/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffParameters.h index 2f535f0..53899f7 100644 --- a/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffParameters.h +++ b/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffParameters.h @@ -47,7 +47,7 @@ class GaffParameters : public MMParameters { /** * @brief Get Lennard-Jones for two atom types t1 and t2. */ - LennardJones getMMLennardJones(std::string t1, std::string t2, double scalingFactor) const; + LennardJones getMMLennardJones(const std::string& t1, const std::string& t2, double scalingFactor) const; // These functions add certain parameters of the MM model void addLennardJones(std::string atomType, LennardJonesParameters lennardJonesParameters); diff --git a/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffPotentialTermsGenerator.cpp b/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffPotentialTermsGenerator.cpp index 90d8a1d..bf5ebd1 100644 --- a/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffPotentialTermsGenerator.cpp +++ b/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffPotentialTermsGenerator.cpp @@ -19,15 +19,8 @@ namespace MolecularMechanics { GaffPotentialTermsGenerator::GaffPotentialTermsGenerator(int nAtoms, const AtomTypesHolder& atomTypes, const IndexedStructuralTopology& topology, const GaffParameters& parameters, - const Utils::PositionCollection& positions, - const double& nonCovalentCutoffRadius, Core::Log& log) - : nAtoms_(nAtoms), - atomTypesHolder_(atomTypes), - topology_(topology), - parameters_(parameters), - positions_(positions), - cutoff_(std::make_shared(nonCovalentCutoffRadius)), - log_(log) { + const Utils::PositionCollection& positions, Core::Log& log) + : nAtoms_(nAtoms), atomTypesHolder_(atomTypes), topology_(topology), parameters_(parameters), positions_(positions), log_(log) { } std::vector GaffPotentialTermsGenerator::getBondedTerms() { @@ -97,45 +90,5 @@ std::vector GaffPotentialTermsGenerator::getImproperDihedralTerms( return improperDihedralList; } -std::vector GaffPotentialTermsGenerator::getLennardJonesTerms(bool applyCutoff) { - std::vector ljList; - - // 0 = excluded, 1 = included, -1 = scaled - Eigen::MatrixXi exclusionType = PotentialTermsHelper::getExclusionTypeMatrix(topology_, nAtoms_); - - for (int a = 0; a < nAtoms_; ++a) { - for (int b = 0; b < a; ++b) { - if (applyCutoff) { - const auto& R = (positions_.row(b) - positions_.row(a)).norm(); - if (R > *cutoff_) - continue; - } - - auto type = exclusionType(a, b); - if (type == 0) - continue; - - auto atomType1 = atomTypesHolder_.getAtomType(a); - auto atomType2 = atomTypesHolder_.getAtomType(b); - - double factor = 1.0; - if (type == -1) - factor = scalingFactorForLennardJonesOneFourTerms_; - LennardJones lj = parameters_.getMMLennardJones(atomType1, atomType2, factor); - LennardJonesTerm term(a, b, lj, cutoff_); - ljList.push_back(term); - } - } - - return ljList; -} - -std::vector GaffPotentialTermsGenerator::getElectrostaticTerms(bool applyCutoff) { - // 0 = excluded, 1 = included, -1 = scaled - Eigen::MatrixXi exclusionType = PotentialTermsHelper::getExclusionTypeMatrix(topology_, nAtoms_); - return PotentialTermsHelper::getElectrostaticTerms(applyCutoff, cutoff_, scalingFactorForElectrostaticOneFourTerms_, - exclusionType, positions_); -} - } // namespace MolecularMechanics } // namespace Scine diff --git a/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffPotentialTermsGenerator.h b/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffPotentialTermsGenerator.h index 613ec28..cc04eca 100644 --- a/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffPotentialTermsGenerator.h +++ b/src/Swoose/Swoose/MolecularMechanics/GAFF/GaffPotentialTermsGenerator.h @@ -11,9 +11,7 @@ #include "../Interactions/AngleTerm.h" #include "../Interactions/BondedTerm.h" #include "../Interactions/DihedralTerm.h" -#include "../Interactions/ElectrostaticTerm.h" #include "../Interactions/ImproperDihedralTerm.h" -#include "../Interactions/LennardJonesTerm.h" #include namespace Scine { @@ -42,12 +40,10 @@ class GaffPotentialTermsGenerator { * @param atomTypes The atom types. * @param topology The topology. * @param parameters The parameters of GAFF. - * @param nonCovalentCutoffRadius The cutoff radius for non-covalent interactions. * @param log The logger. */ GaffPotentialTermsGenerator(int nAtoms, const AtomTypesHolder& atomTypes, const IndexedStructuralTopology& topology, - const GaffParameters& parameters, const Utils::PositionCollection& positions, - const double& nonCovalentCutoffRadius, Core::Log& log); + const GaffParameters& parameters, const Utils::PositionCollection& positions, Core::Log& log); /** @brief Getter for bonded terms */ std::vector getBondedTerms(); @@ -57,10 +53,6 @@ class GaffPotentialTermsGenerator { std::vector getDihedralTerms(); /** @brief Getter for improper dihedral terms. In GAFF, an improper dihedral is treated as a normal dihedral. */ std::vector getImproperDihedralTerms(); - /** @brief Getter for Lennard-Jones terms with decision whether to use a cutoff radius. */ - std::vector getLennardJonesTerms(bool applyCutoff); - /** @brief Getter for electrostatic terms with decision whether to use a cutoff radius. */ - std::vector getElectrostaticTerms(bool applyCutoff); private: int nAtoms_; @@ -68,19 +60,8 @@ class GaffPotentialTermsGenerator { const IndexedStructuralTopology& topology_; const GaffParameters& parameters_; const Utils::PositionCollection& positions_; - /* - * @brief The cutoff radius for the non-covalent interactions. - * - * It is a pointer so that it is updated with the settings. - * This is important since it is passed on from here to the actual DispersionTerms, etc. - * A reference does not work because then the assignment operator of DispersionTerms, etc. vanishes. - */ - std::shared_ptr cutoff_; // The logger Core::Log& log_; - - static constexpr double scalingFactorForElectrostaticOneFourTerms_ = 0.5; - static constexpr double scalingFactorForLennardJonesOneFourTerms_ = 0.5; }; } // namespace MolecularMechanics diff --git a/src/Swoose/Swoose/MolecularMechanics/InteractionExclusion.cpp b/src/Swoose/Swoose/MolecularMechanics/InteractionExclusion.cpp new file mode 100644 index 0000000..ffda731 --- /dev/null +++ b/src/Swoose/Swoose/MolecularMechanics/InteractionExclusion.cpp @@ -0,0 +1,55 @@ +/** + * @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. + */ + +#include "Swoose/MolecularMechanics/InteractionExclusion.h" +#include "Topology/IndexedStructuralTopology.h" + +namespace Scine { +namespace MolecularMechanics { + +InteractionExclusion::InteractionExclusion(unsigned int nAtoms) + : exclusions_(Eigen::SparseMatrix(nAtoms, nAtoms)) { +} + +void InteractionExclusion::setExclusions(const Eigen::SparseMatrix& excludedPairs) { + if (excludedPairs.cols() != exclusions_.cols() || excludedPairs.rows() != exclusions_.rows()) { + throw std::runtime_error("Matrix dimensions are incorrect. The number of atoms in the system is" + " not allowed to change!"); + } + exclusions_ = excludedPairs; +} + +void InteractionExclusion::addExclusions(const Eigen::SparseMatrix& excludedPairs) { + if (excludedPairs.cols() != exclusions_.cols() || excludedPairs.rows() != exclusions_.rows()) { + throw std::runtime_error("Matrix dimensions are incorrect. The added exclusions must refer to the same" + " atom indices as the existing ones!"); + } + exclusions_ = exclusions_ || excludedPairs; +} + +const Eigen::SparseMatrix& InteractionExclusion::getExclusions() { + return exclusions_; +} +void InteractionExclusion::addExclusions(const IndexedStructuralTopology& topology) { + std::vector> triplets; + for (const auto& excluded : topology.getExcludedNonBondedContainer()) { + if (excluded.atom1 > this->exclusions_.cols() || excluded.atom2 > this->exclusions_.cols()) { + throw std::runtime_error("Exclusion lists does not fit to the atom indices!"); + } + triplets.emplace_back(excluded.atom1, excluded.atom2, true); + triplets.emplace_back(excluded.atom2, excluded.atom1, true); + } + Eigen::SparseMatrix exclusions(this->exclusions_.cols(), this->exclusions_.cols()); + exclusions.setFromTriplets(triplets.begin(), triplets.end()); + this->addExclusions(exclusions); +} +void InteractionExclusion::resetExclusions(unsigned int nAtoms) { + exclusions_ = Eigen::SparseMatrix(nAtoms, nAtoms); +} + +} // namespace MolecularMechanics +} // namespace Scine diff --git a/src/Swoose/Swoose/MolecularMechanics/InteractionExclusion.h b/src/Swoose/Swoose/MolecularMechanics/InteractionExclusion.h new file mode 100644 index 0000000..f6ed1e5 --- /dev/null +++ b/src/Swoose/Swoose/MolecularMechanics/InteractionExclusion.h @@ -0,0 +1,68 @@ +/** + * @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 SWOOSE_INTERACTIONEXCLUSION_H +#define SWOOSE_INTERACTIONEXCLUSION_H + +#include + +namespace Scine { +namespace MolecularMechanics { +class IndexedStructuralTopology; + +/** + * @class InteractionExclusion InteractionExclusion.h + * @brief This class provides a framework to store exclusions of interactions between atoms. + * Such exclusions are required for the non-covalent interactions between bonded atoms or + * to eliminate contributions to the energy for the interaction between QM atoms in QM/MM. + * + * The exclusions are stored as a sparse matrix atom x atom of booleans. If the entry is true + * for two atoms, their interaction is excluded. + */ +class InteractionExclusion { + public: + /** + * @brief Constructor. + * @param nAtoms The number of atoms. + * + * Note that the number of atoms may be changed later by calling resetExclusions(...). + */ + InteractionExclusion(unsigned int nAtoms); + /** + * @brief Setter for the exclusion matrix. + * @param excludedPairs The exclusion matrix. + */ + void setExclusions(const Eigen::SparseMatrix& excludedPairs); + /** + * @brief Add new exclusions. + * @param excludedPairs The exclusions. + */ + void addExclusions(const Eigen::SparseMatrix& excludedPairs); + /** + * @brief Add the exclusions encoded in the topology, e.g., 1-3 and 1-2 interactions. + * @param topology The topology. + */ + void addExclusions(const IndexedStructuralTopology& topology); + /** + * @brief Getter for the exclusions. + * @return The exclusions. + */ + const Eigen::SparseMatrix& getExclusions(); + /** + * @brief Resize the exclusion matrix to the given number of atoms and remove all existing exlcuions. + * @param nAtoms The new number of atoms. + */ + void resetExclusions(unsigned int nAtoms); + + private: + Eigen::SparseMatrix exclusions_; +}; + +} // namespace MolecularMechanics +} // namespace Scine + +#endif // SWOOSE_INTERACTIONEXCLUSION_H diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/AnglesEvaluator.h b/src/Swoose/Swoose/MolecularMechanics/Interactions/AnglesEvaluator.h index 1eced95..8aed026 100644 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/AnglesEvaluator.h +++ b/src/Swoose/Swoose/MolecularMechanics/Interactions/AnglesEvaluator.h @@ -19,7 +19,7 @@ class InteractionTermEliminator; namespace MolecularMechanics { -class AngleType; +struct AngleType; class AnglesEvaluator { public: diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/DihedralsEvaluator.h b/src/Swoose/Swoose/MolecularMechanics/Interactions/DihedralsEvaluator.h index 0c80854..507ba23 100644 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/DihedralsEvaluator.h +++ b/src/Swoose/Swoose/MolecularMechanics/Interactions/DihedralsEvaluator.h @@ -20,7 +20,7 @@ class InteractionTermEliminator; namespace MolecularMechanics { -class DihedralType; +struct DihedralType; /** * @class DihedralsEvaluator DihedralsEvaluator.h diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/Electrostatic.h b/src/Swoose/Swoose/MolecularMechanics/Interactions/Electrostatic.h index b5145ab..d36b938 100644 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/Electrostatic.h +++ b/src/Swoose/Swoose/MolecularMechanics/Interactions/Electrostatic.h @@ -9,6 +9,7 @@ #define MOLECULARMECHANICS_ELECTROSTATIC_H #include +#include namespace Scine { namespace MolecularMechanics { diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticEvaluator.cpp b/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticEvaluator.cpp index d0df347..ba37576 100644 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticEvaluator.cpp +++ b/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticEvaluator.cpp @@ -5,26 +5,84 @@ * See LICENSE.txt for details. */ -#include "ElectrostaticEvaluator.h" +#include "Swoose/MolecularMechanics/Interactions/ElectrostaticEvaluator.h" +#include +#include +#include namespace Scine { namespace MolecularMechanics { + ElectrostaticEvaluator::ElectrostaticEvaluator(const Utils::PositionCollection& positions, const std::vector& atomicCharges) - : positions_(positions), atomicCharges_(atomicCharges) { + : InteractionExclusion(positions.rows()), + ScaledInteractions(positions.rows()), + positions_(positions), + atomicCharges_(atomicCharges) { } -double ElectrostaticEvaluator::evaluate(Utils::FullSecondDerivativeCollection& derivatives) { - double energy = 0.0; +double ElectrostaticEvaluator::evaluate(Utils::DerivativeCollection& derivatives) { + const unsigned int nAtoms = this->positions_.rows(); + const Eigen::Matrix allAtoms = Eigen::Matrix::Constant(nAtoms, 1, true); + const Eigen::SparseMatrix& exclusions = this->getExclusions(); + const Eigen::SparseMatrix& scaledTerms = this->getScaledInteractionPairs(); + const double& scalingFactor = this->getInteractionScalingFactor(); - for (const auto& electrostaticTerm : electrostaticTerms_) { - energy += electrostaticTerm.evaluateElectrostaticTerm(positions_, derivatives, atomicCharges_, scalingFactor_); + // Avoid nested parallelization. + omp_set_max_active_levels(1); + const unsigned int nThreads = omp_get_max_threads(); + Eigen::VectorXd energySet = Eigen::VectorXd::Zero(nThreads); + std::vector derivativeSet(nThreads, + Utils::DerivativeCollection(int(nAtoms), derivatives.getOrder())); + for (auto& deriv : derivativeSet) { + deriv.setZero(); } - - return energy; +#pragma omp parallel for schedule(dynamic) + for (unsigned int iAtom = 0; iAtom < nAtoms; ++iAtom) { + const unsigned int threadID = omp_get_thread_num(); + // All atoms for which the interaction is neither scaled nor excluded. + // We must call the method pruned() here to ensure that any entry with false is removed from the sparse vector. + Eigen::SparseVector notScaledOrExcluded = + (allAtoms - (exclusions.col(iAtom) || scaledTerms.col(iAtom))).eval().pruned(); + // Evaluate non-scaled terms. + energySet(threadID) += evaluateTermsForAtom(iAtom, derivativeSet[threadID], notScaledOrExcluded, 1.0); + // Evaluate scaled terms. + const Eigen::SparseVector scaledButExcluded = scaledTerms.col(iAtom) && exclusions.col(iAtom); + energySet(threadID) += evaluateTermsForAtom(iAtom, derivativeSet[threadID], + (scaledTerms.col(iAtom) - scaledButExcluded).pruned(), scalingFactor); + } // for iAtom + for (const auto& deriv : derivativeSet) { + derivatives += deriv; + } + return energySet.sum(); } -void ElectrostaticEvaluator::setElectrostaticTerms(std::vector&& electrostaticTerms) { - electrostaticTerms_ = electrostaticTerms; +double ElectrostaticEvaluator::evaluateTermsForAtom(unsigned int atomIndex, Utils::DerivativeCollection& derivatives, + const Eigen::SparseVector& otherAtoms, double scaling) { + double energyIncrement = 0.0; + const auto positionsI = this->positions_.row(atomIndex); + const double& chargeI = this->atomicCharges_[atomIndex]; + const double totalScalingFactor = scalingFactor_ * scalingFactor_ * scaling; + for (Eigen::SparseVector::InnerIterator it(otherAtoms); it; ++it) { + if (it.row() >= atomIndex) { + break; + } + const unsigned int atomIndexII = it.row(); + const auto positionDifference = this->positions_.row(atomIndexII) - positionsI; + const double distance = positionDifference.norm(); + if (distance > *this->cutOffRadius_) { + continue; + } + const double& chargeII = this->atomicCharges_[atomIndexII]; + const auto contribution = + totalScalingFactor / + Utils::AutomaticDifferentiation::variableWithUnitDerivative(distance) * + (chargeI * chargeII); + const auto value = + Utils::AutomaticDifferentiation::get3Dfrom1D(contribution, positionDifference); + derivatives.addDerivative(int(atomIndex), int(atomIndexII), value); + energyIncrement += value.value(); + } + return energyIncrement; } const std::vector& ElectrostaticEvaluator::getAtomicCharges() { @@ -34,6 +92,9 @@ const std::vector& ElectrostaticEvaluator::getAtomicCharges() { void ElectrostaticEvaluator::setScalingFactor(const double& scalingFactor) { scalingFactor_ = scalingFactor; } +void ElectrostaticEvaluator::setCutOffRadius(std::shared_ptr cutOffRadius) { + cutOffRadius_ = cutOffRadius; +} } // namespace MolecularMechanics } // namespace Scine diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticEvaluator.h b/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticEvaluator.h index 1567b67..f1803cd 100644 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticEvaluator.h +++ b/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticEvaluator.h @@ -8,7 +8,11 @@ #ifndef MOLECULARMECHANICS_ELECTROSTATICEVALUATOR_H #define MOLECULARMECHANICS_ELECTROSTATICEVALUATOR_H -#include "ElectrostaticTerm.h" +#include "Swoose/MolecularMechanics/InteractionExclusion.h" +#include "Swoose/MolecularMechanics/ScaledInteractions.h" +#include +#include +#include #include namespace Scine { @@ -27,7 +31,7 @@ static constexpr double defaultAtomicChargesScalingFactor = 1.0; * @class ElectrostaticEvaluator ElectrostaticEvaluator.h * @brief Class evaluating the total energy and derivatives from the electrostatic interactions. */ -class ElectrostaticEvaluator { +class ElectrostaticEvaluator : public InteractionExclusion, public ScaledInteractions { public: /** * @brief Constructor from positions and partial atomic charges. @@ -36,29 +40,40 @@ class ElectrostaticEvaluator { /** * @brief Function to evaluate and return the total electrostatic energy and update the derivatives. */ - double evaluate(Utils::FullSecondDerivativeCollection& derivatives); - /** - * @brief Sets a vector of instances of the ElectrostaticTerm class. - */ - void setElectrostaticTerms(std::vector&& electrostaticTerms); + double evaluate(Utils::DerivativeCollection& derivatives); /** * @brief Constant accessor for the partial atomic charges. */ const std::vector& getAtomicCharges(); - /** * @brief Setter for the scaling factor for each atomic charge. */ void setScalingFactor(const double& scalingFactor); + /** + * @brief Set the cut off radius for the electrostatic interactions. + * @param cutOffRadius The cut of radius in atomic units. + */ + void setCutOffRadius(std::shared_ptr cutOffRadius); private: // friend class declaration friend class Qmmm::InteractionTermEliminator; const Utils::PositionCollection& positions_; const std::vector& atomicCharges_; - std::vector electrostaticTerms_; // Scaling factor applied to each atomic charge. double scalingFactor_ = Defaults::defaultAtomicChargesScalingFactor; + /** + * @brief Calculate the energy and force contribution for a single atom. + * @param atomIndex The atom index. + * @param derivatives The derivative object to which the gradient contributions are added. + * @param otherAtoms Sparse list of atoms for which the interaction is calculated. + * @param scaling Interaction scaling. + * @return The energy contribution. + */ + double evaluateTermsForAtom(unsigned int atomIndex, Utils::DerivativeCollection& derivatives, + const Eigen::SparseVector& otherAtoms, double scaling); + ///@brief The cut off radius. + std::shared_ptr cutOffRadius_ = std::make_shared(std::numeric_limits::infinity()); }; } // namespace MolecularMechanics diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticTerm.cpp b/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticTerm.cpp deleted file mode 100644 index 31b0f98..0000000 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticTerm.cpp +++ /dev/null @@ -1,50 +0,0 @@ -/** - * @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. - */ - -#include "ElectrostaticTerm.h" -#include -#include - -namespace Scine { -namespace MolecularMechanics { -ElectrostaticTerm::ElectrostaticTerm(AtomIndex firstAtom, AtomIndex secondAtom, const Electrostatic& electrostatic, - std::shared_ptr cutoffRadius) - : firstAtom_(firstAtom), secondAtom_(secondAtom), electrostatic_(electrostatic), cutoffRadius_(cutoffRadius) { -} - -double ElectrostaticTerm::evaluateElectrostaticTerm(const Utils::PositionCollection& positions, - Utils::FullSecondDerivativeCollection& derivatives, - const std::vector& atomicCharges, - const double& scalingFactorForEachCharge) const { - if (this->disabled_) - return 0.0; - - auto R = positions.row(secondAtom_) - positions.row(firstAtom_); - - if (R.norm() > *cutoffRadius_) - return 0.0; - - auto v = Utils::AutomaticDifferentiation::get3Dfrom1D( - electrostatic_.getInteraction(R.norm(), scalingFactorForEachCharge * atomicCharges[firstAtom_], - scalingFactorForEachCharge * atomicCharges[secondAtom_]), - R); - - derivatives.addDerivative(firstAtom_, secondAtom_, v); - - return v.value(); -} - -int ElectrostaticTerm::getFirstAtom() const { - return firstAtom_; -} - -int ElectrostaticTerm::getSecondAtom() const { - return secondAtom_; -} - -} // namespace MolecularMechanics -} // namespace Scine diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticTerm.h b/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticTerm.h deleted file mode 100644 index 779efa5..0000000 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/ElectrostaticTerm.h +++ /dev/null @@ -1,60 +0,0 @@ -/** - * @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 MOLECULARMECHANICS_ELECTROSTATICTERM_H -#define MOLECULARMECHANICS_ELECTROSTATICTERM_H - -#include "Electrostatic.h" -#include "InteractionTermBase.h" -#include -#include -#include - -namespace Scine { - -namespace Utils { -class FullSecondDerivativeCollection; -} // namespace Utils - -namespace MolecularMechanics { -/** - * @class ElectrostaticTerm ElectrostaticTerm.h - * @brief Class evaluating electrostatic interaction between two atoms. - */ -class ElectrostaticTerm : public InteractionTermBase { - public: - using AtomIndex = int; - - /** - * @brief Constructor from two atom indices and an instance of the Electrostatic class. - */ - ElectrostaticTerm(AtomIndex firstAtom, AtomIndex secondAtom, const Electrostatic& electrostatic, - std::shared_ptr cutoffRadius); - - /** - * @brief Evaluates energy contribution and adds the derivatives. - */ - double evaluateElectrostaticTerm(const Utils::PositionCollection& positions, Utils::FullSecondDerivativeCollection& derivatives, - const std::vector& atomicCharges, const double& scalingFactorForEachCharge) const; - - /** @brief Getter for index of first atom. */ - int getFirstAtom() const; - /** @brief Getter for index of second atom. */ - int getSecondAtom() const; - - private: - AtomIndex firstAtom_, secondAtom_; - Electrostatic electrostatic_; - // Pointer so that it is updated whenever the settings are updated. - // Reference does not work because then the class loses its assignment operator. - std::shared_ptr cutoffRadius_; -}; - -} // namespace MolecularMechanics -} // namespace Scine - -#endif // MOLECULARMECHANICS_ELECTROSTATICTERM_H diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/ImproperDihedralTerm.cpp b/src/Swoose/Swoose/MolecularMechanics/Interactions/ImproperDihedralTerm.cpp index cc6cc4f..768f963 100644 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/ImproperDihedralTerm.cpp +++ b/src/Swoose/Swoose/MolecularMechanics/Interactions/ImproperDihedralTerm.cpp @@ -100,7 +100,6 @@ double ImproperDihedralTerm::evaluateImproperDihedralTerm(const Utils::PositionC h4.setFirst3D(firstDer4); double theta = getTheta(A, B, G); - // std::cout << "Theta: " << theta * 57.2958 << std::endl; auto result = improperDihedral_.getInteraction(theta); // Apply chain rule to get derivatives with respect to the energy diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/ImproperDihedralsEvaluator.h b/src/Swoose/Swoose/MolecularMechanics/Interactions/ImproperDihedralsEvaluator.h index 4c2b432..637862a 100644 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/ImproperDihedralsEvaluator.h +++ b/src/Swoose/Swoose/MolecularMechanics/Interactions/ImproperDihedralsEvaluator.h @@ -20,7 +20,7 @@ class InteractionTermEliminator; namespace MolecularMechanics { -class ImproperDihedralType; +struct ImproperDihedralType; /** * @class ImproperDihedralsEvaluator ImproperDihedralsEvaluator.h diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesEvaluator.cpp b/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesEvaluator.cpp index 75f61be..8c3e247 100644 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesEvaluator.cpp +++ b/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesEvaluator.cpp @@ -6,25 +6,90 @@ */ #include "LennardJonesEvaluator.h" +#include "Swoose/MolecularMechanics/AtomTypesHolder.h" +#include "Swoose/MolecularMechanics/GAFF/GaffParameters.h" +#include +#include namespace Scine { namespace MolecularMechanics { -LennardJonesEvaluator::LennardJonesEvaluator(const Utils::PositionCollection& positions) : positions_(positions) { +LennardJonesEvaluator::LennardJonesEvaluator(const Utils::PositionCollection& positions) + : InteractionExclusion(positions.rows()), ScaledInteractions(positions.rows()), positions_(positions) { } -double LennardJonesEvaluator::evaluate(Utils::FullSecondDerivativeCollection& derivatives) { - double energy = 0.0; +double LennardJonesEvaluator::evaluate(Utils::DerivativeCollection& derivatives) { + if (!parameters_ || !atomTypesHolder_) { + throw std::runtime_error("Parameters or atom types not available for Lennard Jones energy evaluation."); + } + const unsigned int nAtoms = this->positions_.rows(); + const Eigen::Matrix allAtoms = Eigen::Matrix::Constant(nAtoms, 1, true); + const Eigen::SparseMatrix& exclusions = this->getExclusions(); + const Eigen::SparseMatrix& scaledTerms = this->getScaledInteractionPairs(); + const double& scalingFactor = this->getInteractionScalingFactor(); - for (auto& lj : ljTerms_) { - energy += lj.evaluateLennardJonesTerm(positions_, derivatives); + // Avoid nested parallelization. + omp_set_max_active_levels(1); + const unsigned int nThreads = omp_get_max_threads(); + Eigen::VectorXd energySet = Eigen::VectorXd::Zero(nThreads); + std::vector derivativeSet(nThreads, + Utils::DerivativeCollection(int(nAtoms), derivatives.getOrder())); + for (auto& deriv : derivativeSet) { + deriv.setZero(); + } +#pragma omp parallel for schedule(dynamic) + for (unsigned int iAtom = 0; iAtom < nAtoms; ++iAtom) { + const unsigned int threadID = omp_get_thread_num(); + // All atoms for which the interaction is neither scaled nor excluded. + // We must call the method pruned() here to ensure that any entry with false is removed from the sparse vector. + Eigen::SparseVector notScaledOrExcluded = + (allAtoms - (exclusions.col(iAtom) || scaledTerms.col(iAtom))).eval().pruned(); + // Evaluate non-scaled terms. + energySet(threadID) += evaluateTermsForAtom(iAtom, derivativeSet[threadID], notScaledOrExcluded, 1.0); + // Evaluate scaled terms. + const Eigen::SparseVector scaledButExcluded = scaledTerms.col(iAtom) && exclusions.col(iAtom); + energySet(threadID) += evaluateTermsForAtom(iAtom, derivativeSet[threadID], + (scaledTerms.col(iAtom) - scaledButExcluded).pruned(), scalingFactor); + } // for iAtom + for (const auto& deriv : derivativeSet) { + derivatives += deriv; } + return energySet.sum(); +} - return energy; +double LennardJonesEvaluator::evaluateTermsForAtom(unsigned int atomIndex, Utils::DerivativeCollection& derivatives, + const Eigen::SparseVector& otherAtoms, double scaling) { + double energyIncrement = 0.0; + const auto positionsI = this->positions_.row(atomIndex); + const auto& atomTypeI = this->atomTypesHolder_->getAtomType(atomIndex); + for (Eigen::SparseVector::InnerIterator it(otherAtoms); it; ++it) { + if (it.row() >= atomIndex) { + break; + } + const unsigned int atomIndexII = it.row(); + const auto positionDifference = this->positions_.row(atomIndexII) - positionsI; + const double distance = positionDifference.norm(); + if (distance > *this->cutOffRadius_) { + continue; + } + const auto& atomTypeII = this->atomTypesHolder_->getAtomType(atomIndexII); + const auto combinedParameters = parameters_->getMMLennardJones(atomTypeI, atomTypeII, scaling); + const auto value = Utils::AutomaticDifferentiation::get3Dfrom1D( + combinedParameters.getInteraction(distance), positionDifference); + derivatives.addDerivative(int(atomIndex), int(atomIndexII), value); + energyIncrement += value.value(); + } + return energyIncrement; } -void LennardJonesEvaluator::setLennardJonesTerms(std::vector&& ljTerms) { - ljTerms_ = ljTerms; +void LennardJonesEvaluator::setCutOffRadius(std::shared_ptr cutOffRadius) { + cutOffRadius_ = cutOffRadius; +} +void LennardJonesEvaluator::setParameters(std::shared_ptr parameters) { + parameters_ = parameters; +} +void LennardJonesEvaluator::setAtomTypesHolder(std::shared_ptr atomTypesHolder) { + atomTypesHolder_ = atomTypesHolder; } } // namespace MolecularMechanics diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesEvaluator.h b/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesEvaluator.h index 381aa3d..d4a39d6 100644 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesEvaluator.h +++ b/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesEvaluator.h @@ -8,8 +8,10 @@ #ifndef MOLECULARMECHANICS_LENNARDJONESEVALUATOR_H #define MOLECULARMECHANICS_LENNARDJONESEVALUATOR_H -#include "LennardJonesTerm.h" +#include "Swoose/MolecularMechanics/InteractionExclusion.h" +#include "Swoose/MolecularMechanics/ScaledInteractions.h" #include +#include #include namespace Scine { @@ -19,15 +21,17 @@ class InteractionTermEliminator; } // namespace Qmmm namespace Utils { -class FullSecondDerivativeCollection; +class DerivativeCollection; } // namespace Utils namespace MolecularMechanics { +class GaffParameters; +class AtomTypesHolder; /** * @brief LennardJonesEvaluator LennardJonesEvaluator.h * @brief This class evaluates the overall energy and derivatives of Lennard-Jones interactions. */ -class LennardJonesEvaluator { +class LennardJonesEvaluator : public InteractionExclusion, public ScaledInteractions { public: /** * @brief Constructor from positions. @@ -36,17 +40,44 @@ class LennardJonesEvaluator { /** * @brief This function evaluates and returns the energy for all LJ interactions and updates the derivatives. */ - double evaluate(Utils::FullSecondDerivativeCollection& derivatives); + double evaluate(Utils::DerivativeCollection& derivatives); /** - * @brief Sets a vector of instances of the LennardJonesTerm class. + * @brief Set the cut off radius for the electrostatic interactions. + * @param cutOffRadius The cut of radius in atomic units. */ - void setLennardJonesTerms(std::vector&& ljTerms); + void setCutOffRadius(std::shared_ptr cutOffRadius); + /** + * @brief Setter for the parameter object. + * @param parameters The parameters. + */ + void setParameters(std::shared_ptr parameters); + /** + * @brief Setter for the atom types. + * @param atomTypesHolder The atom tyoes. + */ + void setAtomTypesHolder(std::shared_ptr atomTypesHolder); private: // friend class declaration friend class Qmmm::InteractionTermEliminator; const Utils::PositionCollection& positions_; - std::vector ljTerms_; + + ///@brief The Parameters. + std::shared_ptr parameters_ = nullptr; + ///@brief The atom types. + std::shared_ptr atomTypesHolder_ = nullptr; + ///@brief The cut off radius. + std::shared_ptr cutOffRadius_ = std::make_shared(std::numeric_limits::infinity()); + /** + * @brief Calculate the energy and force contribution for a single atom. + * @param atomIndex The atom index. + * @param derivatives The derivative object to which the gradient contributions are added. + * @param otherAtoms Sparse list of atoms for which the interaction is calculated. + * @param scaling Interaction scaling. + * @return The energy contribution. + */ + double evaluateTermsForAtom(unsigned int atomIndex, Utils::DerivativeCollection& derivatives, + const Eigen::SparseVector& otherAtoms, double scaling); }; } // namespace MolecularMechanics diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesTerm.cpp b/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesTerm.cpp deleted file mode 100644 index 0a4d9b8..0000000 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesTerm.cpp +++ /dev/null @@ -1,46 +0,0 @@ -/** - * @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. - */ - -#include "LennardJonesTerm.h" -#include -#include - -namespace Scine { -namespace MolecularMechanics { - -LennardJonesTerm::LennardJonesTerm(AtomIndex firstAtom, AtomIndex secondAtom, const LennardJones& lj, - std::shared_ptr cutoffRadius) - : firstAtom_(firstAtom), secondAtom_(secondAtom), lj_(lj), cutoffRadius_(cutoffRadius) { -} - -double LennardJonesTerm::evaluateLennardJonesTerm(const Utils::PositionCollection& positions, - Utils::FullSecondDerivativeCollection& derivatives) const { - if (this->disabled_) - return 0.0; - - auto R = positions.row(secondAtom_) - positions.row(firstAtom_); - - if (R.norm() > *cutoffRadius_) - return 0.0; - - auto v = Utils::AutomaticDifferentiation::get3Dfrom1D(lj_.getInteraction(R.norm()), R); - - derivatives.addDerivative(firstAtom_, secondAtom_, v); - - return v.value(); -} - -int LennardJonesTerm::getFirstAtom() const { - return firstAtom_; -} - -int LennardJonesTerm::getSecondAtom() const { - return secondAtom_; -} - -} // namespace MolecularMechanics -} // namespace Scine diff --git a/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesTerm.h b/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesTerm.h deleted file mode 100644 index a356f4d..0000000 --- a/src/Swoose/Swoose/MolecularMechanics/Interactions/LennardJonesTerm.h +++ /dev/null @@ -1,59 +0,0 @@ -/** - * @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 MOLECULARMECHANICS_LENNARDJONESTERM_H -#define MOLECULARMECHANICS_LENNARDJONESTERM_H - -#include "InteractionTermBase.h" -#include "LennardJones.h" -#include -#include -#include - -namespace Scine { - -namespace Utils { -class FullSecondDerivativeCollection; -} // namespace Utils - -namespace MolecularMechanics { -/** - * @class LennardJonesTerm LennardJonesTerm.h - * @brief Class evaluating electrostatic interaction between two atoms. - */ -class LennardJonesTerm : public InteractionTermBase { - public: - using AtomIndex = int; - - /** - * @brief Constructor from two atom indices and an instance of the LennardJones class. - */ - LennardJonesTerm(AtomIndex firstAtom, AtomIndex secondAtom, const LennardJones& lj, std::shared_ptr cutoffRadius); - - /** - * @brief Evaluates energy contribution and adds the derivatives. - */ - double evaluateLennardJonesTerm(const Utils::PositionCollection& positions, - Utils::FullSecondDerivativeCollection& derivatives) const; - - /** @brief Getter for index of first atom. */ - int getFirstAtom() const; - /** @brief Getter for index of second atom. */ - int getSecondAtom() const; - - private: - AtomIndex firstAtom_, secondAtom_; - LennardJones lj_; - // Pointer so that it is updated whenever the settings are updated. - // Reference does not work because then the class loses its assignment operator. - std::shared_ptr cutoffRadius_; -}; - -} // namespace MolecularMechanics -} // namespace Scine - -#endif // MOLECULARMECHANICS_LENNARDJONESTERM_H diff --git a/src/Swoose/Swoose/MolecularMechanics/Parameters/LennardJonesParameters.cpp b/src/Swoose/Swoose/MolecularMechanics/Parameters/LennardJonesParameters.cpp index 482dbef..071a96a 100644 --- a/src/Swoose/Swoose/MolecularMechanics/Parameters/LennardJonesParameters.cpp +++ b/src/Swoose/Swoose/MolecularMechanics/Parameters/LennardJonesParameters.cpp @@ -21,6 +21,9 @@ LennardJonesParameters::LennardJonesParameters(double vdwRadius, double wellDept */ LennardJones LennardJonesParameters::toMMLennardJones(const LennardJonesParameters& otherLjParameters, double scalingFactor) const { + // TODO: It is extremely inefficient to do the unit conversion and scaling here. + // Either, we first convert all parameters or we do the conversion after summing over all + // terms. double rij = (vdwRadius_ + otherLjParameters.vdwRadius_) * Utils::Constants::bohr_per_angstrom; double eij = std::sqrt(wellDepth_ * otherLjParameters.wellDepth_) * Utils::Constants::hartree_per_kCalPerMol * scalingFactor; diff --git a/src/Swoose/Swoose/MolecularMechanics/PotentialTermsHelper.cpp b/src/Swoose/Swoose/MolecularMechanics/PotentialTermsHelper.cpp index 08f4e2e..8237caa 100644 --- a/src/Swoose/Swoose/MolecularMechanics/PotentialTermsHelper.cpp +++ b/src/Swoose/Swoose/MolecularMechanics/PotentialTermsHelper.cpp @@ -77,37 +77,6 @@ std::vector getAngleTerms(const IndexedStructuralTopology& topology, return angleList; } -std::vector -getElectrostaticTerms(bool applyCutoff, std::shared_ptr cutoffRadius, double scalingFactorForOneFourTerms, - const Eigen::MatrixXi& exclusionTypeMatrix, const Utils::PositionCollection& positions) { - std::vector electrostaticList; - assert(positions.rows() == exclusionTypeMatrix.rows()); - assert(positions.rows() == exclusionTypeMatrix.cols()); - - for (int a = 0; a < positions.rows(); ++a) { - for (int b = 0; b < a; ++b) { - if (applyCutoff) { - const auto& R = (positions.row(b) - positions.row(a)).norm(); - if (R > *cutoffRadius) - continue; - } - - auto type = exclusionTypeMatrix(a, b); - if (type == 0) - continue; - - double factor = 1; - if (type == -1) - factor = scalingFactorForOneFourTerms; - Electrostatic m(factor); - ElectrostaticTerm term(a, b, m, cutoffRadius); - electrostaticList.push_back(term); - } - } - - return electrostaticList; -} - } // namespace PotentialTermsHelper } // namespace MolecularMechanics } // namespace Scine diff --git a/src/Swoose/Swoose/MolecularMechanics/PotentialTermsHelper.h b/src/Swoose/Swoose/MolecularMechanics/PotentialTermsHelper.h index b893b9f..73123ee 100644 --- a/src/Swoose/Swoose/MolecularMechanics/PotentialTermsHelper.h +++ b/src/Swoose/Swoose/MolecularMechanics/PotentialTermsHelper.h @@ -10,7 +10,6 @@ #include "Interactions/AngleTerm.h" #include "Interactions/BondedTerm.h" -#include "Interactions/ElectrostaticTerm.h" #include #include @@ -51,22 +50,6 @@ std::vector getBondedTerms(const IndexedStructuralTopology& topology */ std::vector getAngleTerms(const IndexedStructuralTopology& topology, const MMParameters& parameters, const AtomTypesHolder& atomTypesHolder); - -/** - * @brief Getter for the vector of electrostatic terms. This is the same for GAFF and SFAM. - * @param applyCutoff Whether to apply a distance cutoff. - * @param cutoffRadius The distance cutoff in bohr. - * @param scalingFactorForOneFourTerms The factor by which the interaction of a - * 1,4-bonded pair of atoms should be scaled. - * @param exclusionTypeMatrix Exclusion-type matrix, that encodes for each atom pair how its - * non-covalent contribution is handled. - * @param positions The positions of the structure. - * @return Vector of electrostatic terms. - */ -std::vector -getElectrostaticTerms(bool applyCutoff, std::shared_ptr cutoffRadius, double scalingFactorForOneFourTerms, - const Eigen::MatrixXi& exclusionTypeMatrix, const Utils::PositionCollection& positions); - } // namespace PotentialTermsHelper } // namespace MolecularMechanics } // namespace Scine diff --git a/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamMolecularMechanicsCalculator.cpp b/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamMolecularMechanicsCalculator.cpp index 55e67fa..14374f6 100644 --- a/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamMolecularMechanicsCalculator.cpp +++ b/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamMolecularMechanicsCalculator.cpp @@ -311,7 +311,13 @@ void SfamMolecularMechanicsCalculator::generatePotentialTerms(const SfamParamete if (!onlyCalculateBondedContribution_) { dispersionEvaluator_->setDispersionTerms(generator.getDispersionTerms(applyCutoffDuringInitialization_)); repulsionEvaluator_->setRepulsionTerms(generator.getRepulsionTerms(applyCutoffDuringInitialization_)); - electrostaticEvaluator_->setElectrostaticTerms(generator.getElectrostaticTerms(applyCutoffDuringInitialization_)); + + electrostaticEvaluator_->setCutOffRadius(std::make_shared(nonCovalentCutoffRadius_)); + electrostaticEvaluator_->resetExclusions(structure_.size()); + electrostaticEvaluator_->addExclusions(topology); + // No Scaling of any interactions (e.g., 1-4 interactions). Therefore, we only resize the scaling information + // and leave it empty. + electrostaticEvaluator_->resetScaledInteractions(structure_.size()); } } diff --git a/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamMolecularMechanicsCalculator.h b/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamMolecularMechanicsCalculator.h index cee84a7..ff9e158 100644 --- a/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamMolecularMechanicsCalculator.h +++ b/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamMolecularMechanicsCalculator.h @@ -34,10 +34,10 @@ class InteractionTermEliminator; namespace MolecularMechanics { class SfamParameters; class AtomTypesHolder; -class BondType; -class AngleType; -class DihedralType; -class ImproperDihedralType; +struct BondType; +struct AngleType; +struct DihedralType; +struct ImproperDihedralType; /** * @class SfamMolecularMechanicsCalculator SfamMolecularMechanicsCalculator.h * @brief Calculator for the SFAM Molecular Mechanics method. diff --git a/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamParameters.h b/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamParameters.h index 1be9927..61d64cb 100644 --- a/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamParameters.h +++ b/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamParameters.h @@ -34,7 +34,7 @@ class AtomTypesHolder; * @brief Class containing the parameters for SFAM's MM model obtained after parsing a SFAM parameter file. * The angle and bond parameters are handled by the base class MMParameters. */ -class SfamParameters : public MMParameters { +class SfamParameters final : public MMParameters { public: /** @brief Checks whether the SFAM parameters are valid. */ bool sanityCheck(const AtomTypesHolder& atomTypes) const; diff --git a/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamPotentialTermsGenerator.cpp b/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamPotentialTermsGenerator.cpp index 34c51c3..47b59ca 100644 --- a/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamPotentialTermsGenerator.cpp +++ b/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamPotentialTermsGenerator.cpp @@ -164,13 +164,6 @@ std::vector SfamPotentialTermsGenerator::getRepulsionTerms(bool a return repulsionList; } -std::vector SfamPotentialTermsGenerator::getElectrostaticTerms(bool applyCutoff) { - // 0 = excluded, 1 = included, -1 = scaled - Eigen::MatrixXi exclusionType = PotentialTermsHelper::getExclusionTypeMatrix(topology_, nAtoms_); - return PotentialTermsHelper::getElectrostaticTerms(applyCutoff, cutoff_, scalingFactorForElectrostaticOneFourTerms_, - exclusionType, positions_); -} - std::vector SfamPotentialTermsGenerator::getHydrogenBondTerms() { std::vector hydrogenBondList; for (const auto& hydrogenBond : topology_.getHydrogenBondContainer()) { diff --git a/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamPotentialTermsGenerator.h b/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamPotentialTermsGenerator.h index 229b12b..f651c76 100644 --- a/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamPotentialTermsGenerator.h +++ b/src/Swoose/Swoose/MolecularMechanics/SFAM/SfamPotentialTermsGenerator.h @@ -12,7 +12,6 @@ #include "../Interactions/BondedTerm.h" #include "../Interactions/DihedralTerm.h" #include "../Interactions/DispersionTerm.h" -#include "../Interactions/ElectrostaticTerm.h" #include "../Interactions/HydrogenBondTerm.h" #include "../Interactions/ImproperDihedralTerm.h" #include "../Interactions/RepulsionTerm.h" @@ -21,7 +20,7 @@ namespace Scine { namespace Core { -class Log; +struct Log; } // namespace Core namespace MolecularMechanics { @@ -63,8 +62,6 @@ class SfamPotentialTermsGenerator { std::vector getDispersionTerms(bool applyCutoff); /** @brief Getter for Pauli repulsion terms with decision whether to use a cutoff radius. */ std::vector getRepulsionTerms(bool applyCutoff); - /** @brief Getter for electrostatic terms with decision whether to use a cutoff radius. */ - std::vector getElectrostaticTerms(bool applyCutoff); /** @brief Getter for hydrogen bond terms */ std::vector getHydrogenBondTerms(); diff --git a/src/Swoose/Swoose/MolecularMechanics/ScaledInteractions.cpp b/src/Swoose/Swoose/MolecularMechanics/ScaledInteractions.cpp new file mode 100644 index 0000000..abb56ae --- /dev/null +++ b/src/Swoose/Swoose/MolecularMechanics/ScaledInteractions.cpp @@ -0,0 +1,60 @@ +/** + * @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. + */ + +#include "Swoose/MolecularMechanics/ScaledInteractions.h" +#include "Topology/IndexedStructuralTopology.h" + +namespace Scine { +namespace MolecularMechanics { + +ScaledInteractions::ScaledInteractions(unsigned int nAtoms) + : scaledAtomPairs_(Eigen::SparseMatrix(nAtoms, nAtoms)) { +} + +void ScaledInteractions::setScaledInteractionPairs(const Eigen::SparseMatrix& pairs) { + if (pairs.cols() != scaledAtomPairs_.cols() || pairs.rows() != scaledAtomPairs_.rows()) { + throw std::runtime_error("Matrix dimensions are incorrect. The number of atoms in the system is" + " not allowed to change!"); + } + scaledAtomPairs_ = pairs; +} + +void ScaledInteractions::addScaledInteractionPairs(const Eigen::SparseMatrix& pairs) { + if (pairs.cols() != scaledAtomPairs_.cols() || pairs.rows() != scaledAtomPairs_.rows()) { + throw std::runtime_error("Matrix dimensions are incorrect. The added parameter scaling must refer to the same" + " atom indices as the existing ones!"); + } + scaledAtomPairs_ = scaledAtomPairs_ || pairs; +} + +const Eigen::SparseMatrix& ScaledInteractions::getScaledInteractionPairs() { + return scaledAtomPairs_; +} + +void ScaledInteractions::setInteractionScalingFactor(const double& factor) { + scalingFactor_ = factor; +} + +double ScaledInteractions::getInteractionScalingFactor() { + return scalingFactor_; +} + +void ScaledInteractions::addScaledInteractionPairs(const IndexedStructuralTopology& topology) { + std::vector> triplets; + for (const auto& scaled : topology.getScaledNonBondedContainer()) { + triplets.emplace_back(scaled.atom1, scaled.atom2, true); + triplets.emplace_back(scaled.atom2, scaled.atom1, true); + } + Eigen::SparseMatrix scaling(this->scaledAtomPairs_.cols(), this->scaledAtomPairs_.cols()); + scaling.setFromTriplets(triplets.begin(), triplets.end()); + this->addScaledInteractionPairs(scaling); +} +void ScaledInteractions::resetScaledInteractions(unsigned int nAtoms) { + scaledAtomPairs_ = Eigen::SparseMatrix(nAtoms, nAtoms); +} +} // namespace MolecularMechanics +} // namespace Scine diff --git a/src/Swoose/Swoose/MolecularMechanics/ScaledInteractions.h b/src/Swoose/Swoose/MolecularMechanics/ScaledInteractions.h new file mode 100644 index 0000000..b0b1e3b --- /dev/null +++ b/src/Swoose/Swoose/MolecularMechanics/ScaledInteractions.h @@ -0,0 +1,73 @@ +/** + * @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 SWOOSE_SCALEDINTERACTIONS_H +#define SWOOSE_SCALEDINTERACTIONS_H + +#include + +namespace Scine { +namespace MolecularMechanics { +class IndexedStructuralTopology; + +/** + * @class ScaledInteractions ScaledInteractions.h + * @brief This class keeps track which pair-wise interaction between atoms is scaled. + * The scaling information is stored as a sparse matrix atoms x atoms of booleans. + */ +class ScaledInteractions { + public: + /** + * @brief Constructor. + * @param nAtoms The number of atoms. + */ + ScaledInteractions(unsigned int nAtoms); + /** + * @brief Setter for the scaling information. + * @param pairs The scaled atom pairs. + */ + void setScaledInteractionPairs(const Eigen::SparseMatrix& pairs); + /** + * @brief Add new scaled interactions. + * @param pairs The additional scaled atom pairs. + */ + void addScaledInteractionPairs(const Eigen::SparseMatrix& pairs); + /** + * @brief Add new scaled interactions from the topology, e.g., 1-4 interactions. + * @param topology The topology. + */ + void addScaledInteractionPairs(const IndexedStructuralTopology& topology); + /** + * @brief Getter for scaling information. + * @return The scaling matrix. + */ + const Eigen::SparseMatrix& getScaledInteractionPairs(); + /** + * @brief Setter for the interaction scaling factor. + * @param factor The factor. + */ + void setInteractionScalingFactor(const double& factor); + /** + * @brief Getter for the interaction scaling factor. + * @return The interaction scaling factor. + */ + double getInteractionScalingFactor(); + /** + * @brief Resize the scaling information matrix to the nem number of atoms. + * @param nAtoms The number of atoms. + */ + void resetScaledInteractions(unsigned int nAtoms); + + private: + Eigen::SparseMatrix scaledAtomPairs_; + double scalingFactor_ = 0.5; +}; + +} // namespace MolecularMechanics +} // namespace Scine + +#endif // SWOOSE_SCALEDINTERACTIONS_H diff --git a/src/Swoose/Swoose/QMMM/InteractionTermEliminator.cpp b/src/Swoose/Swoose/QMMM/InteractionTermEliminator.cpp index 049550c..5ee308a 100644 --- a/src/Swoose/Swoose/QMMM/InteractionTermEliminator.cpp +++ b/src/Swoose/Swoose/QMMM/InteractionTermEliminator.cpp @@ -36,7 +36,7 @@ void InteractionTermEliminator::eliminateInteractionTerms(bool electrostaticEmbe throw std::runtime_error("GAFF Calculator could not be casted to derived class."); eliminateSharedInteractionTerms(electrostaticEmbedding, *calc); eliminateDihedralTerms(calc->improperDihedralsEvaluator_->dihedrals_, true); - eliminateLennardJonesTerms(calc->lennardJonesEvaluator_->ljTerms_); + eliminateLennardJonesTerms(*calc->lennardJonesEvaluator_); } else { throw std::runtime_error("The given MM model is not supported by the interaction term eliminator."); @@ -48,22 +48,27 @@ void InteractionTermEliminator::eliminateSharedInteractionTerms(bool electrostat eliminateBondedTerms(calculator.bondsEvaluator_->bonds_); eliminateAngleTerms(calculator.anglesEvaluator_->angles_); eliminateDihedralTerms(calculator.dihedralsEvaluator_->dihedrals_); - eliminateElectrostaticTerms(calculator.electrostaticEvaluator_->electrostaticTerms_, electrostaticEmbedding); + eliminateElectrostaticTerms(*calculator.electrostaticEvaluator_, electrostaticEmbedding); } void InteractionTermEliminator::eliminateTerm(MolecularMechanics::InteractionTermBase& term, const std::vector& atomsInTerm, int allowedInQmRegion) { + if (termToEliminate(atomsInTerm, allowedInQmRegion)) { + term.disable(); + } +} + +bool InteractionTermEliminator::termToEliminate(const std::vector& atomsInTerm, int allowedInQmRegion) { int inQmRegion = 0; for (const auto& atom : atomsInTerm) { if (isQmAtom(atom)) inQmRegion++; } if (inQmRegion > allowedInQmRegion) - term.disable(); - if (eliminateEnvironmentOnlyTerms_) { - if (inQmRegion == 0) - term.disable(); - } + return true; + if (eliminateEnvironmentOnlyTerms_ && inQmRegion == 0) + return true; + return false; } bool InteractionTermEliminator::isQmAtom(int index) { @@ -128,11 +133,43 @@ void InteractionTermEliminator::eliminateRepulsionTerms(std::vector& ljTerms) { - for (auto& ljTerm : ljTerms) { - std::vector atoms = {ljTerm.getFirstAtom(), ljTerm.getSecondAtom()}; - eliminateTerm(ljTerm, atoms, 1); +void InteractionTermEliminator::eliminateLennardJonesTerms(MolecularMechanics::LennardJonesEvaluator& lennardJonesEvaluator) { + const unsigned int nAtoms = calculator_->getStructure()->size(); + originalLJExclusions_ = lennardJonesEvaluator.getExclusions(); + std::vector> triplets; + for (unsigned int iAtom = 0; iAtom < nAtoms; ++iAtom) { + for (unsigned int jAtom = 0; jAtom < iAtom; ++jAtom) { + if (termToEliminate({int(iAtom), int(jAtom)}, 1)) { + triplets.emplace_back(iAtom, jAtom, true); + triplets.emplace_back(jAtom, iAtom, true); + } + } + } + Eigen::SparseMatrix exclusions(nAtoms, nAtoms); + exclusions.setFromTriplets(triplets.begin(), triplets.end()); + lennardJonesEvaluator.addExclusions(exclusions); +} + +void InteractionTermEliminator::eliminateElectrostaticTerms(MolecularMechanics::ElectrostaticEvaluator& electrostaticEvaluator, + bool electrostaticEmbedding) { + const unsigned int nAtoms = calculator_->getStructure()->size(); + originalElectrostaticExclusions_ = electrostaticEvaluator.getExclusions(); + std::vector> triplets; + // If electrostatic embedding is switched on: + // Eliminate if at least one atom is in the QM region, + // because QM-MM electrostatic interaction is covered by the electrostatic embedding + const int nQMAtomsAllowed = (electrostaticEmbedding) ? 0 : 1; + for (unsigned int iAtom = 0; iAtom < nAtoms; ++iAtom) { + for (unsigned int jAtom = 0; jAtom < iAtom; ++jAtom) { + if (termToEliminate({int(iAtom), int(jAtom)}, nQMAtomsAllowed)) { + triplets.emplace_back(iAtom, jAtom, true); + triplets.emplace_back(jAtom, iAtom, true); + } + } } + Eigen::SparseMatrix exclusions(nAtoms, nAtoms); + exclusions.setFromTriplets(triplets.begin(), triplets.end()); + electrostaticEvaluator.addExclusions(exclusions); } void InteractionTermEliminator::eliminateHydrogenBondTerms(std::vector& hydrogenBondTerms) { @@ -146,20 +183,6 @@ void InteractionTermEliminator::eliminateHydrogenBondTerms(std::vector& electrostaticTerms, - bool electrostaticEmbedding) { - for (auto& electrostaticTerm : electrostaticTerms) { - std::vector atoms = {electrostaticTerm.getFirstAtom(), electrostaticTerm.getSecondAtom()}; - // If electrostatic embedding is switched on: - // Eliminate if at least one atom is in the QM region, - // because QM-MM electrostatic interaction is covered by the electrostatic embedding - if (electrostaticEmbedding) - eliminateTerm(electrostaticTerm, atoms, 0); - else // For mechanical embedding: If one atom is in the QM region -> don't eliminate this term - eliminateTerm(electrostaticTerm, atoms, 1); - } -} - void InteractionTermEliminator::reset() { auto calculatorName = calculator_->name(); if (calculatorName == "SFAM") { @@ -178,7 +201,12 @@ void InteractionTermEliminator::reset() { throw std::runtime_error("GAFF Calculator could not be casted to derived class."); enableSharedInteractionTerms(*calc); enableTerms(calc->improperDihedralsEvaluator_->dihedrals_); - enableTerms(calc->lennardJonesEvaluator_->ljTerms_); + // The exclusion reset is only possible if we know the original exclusion state. + // If it is unknown at this state, this class never changed any exclusions. + // We do not need to do anything. + if (originalLJExclusions_.cols() != 0) { + calc->lennardJonesEvaluator_->setExclusions(originalLJExclusions_); + } } else { throw std::runtime_error("The given MM model is not supported by the interaction term eliminator."); @@ -197,7 +225,12 @@ void InteractionTermEliminator::enableSharedInteractionTerms(CalculatorType& cal enableTerms(calculator.bondsEvaluator_->bonds_); enableTerms(calculator.anglesEvaluator_->angles_); enableTerms(calculator.dihedralsEvaluator_->dihedrals_); - enableTerms(calculator.electrostaticEvaluator_->electrostaticTerms_); + // The exclusion reset is only possible if we know the original exclusion state. + // If it is unknown at this state, this class never changed any exclusions. + // We do not need to do anything. + if (originalElectrostaticExclusions_.cols() > 0) { + calculator.electrostaticEvaluator_->setExclusions(this->originalElectrostaticExclusions_); + } } } // namespace Qmmm diff --git a/src/Swoose/Swoose/QMMM/InteractionTermEliminator.h b/src/Swoose/Swoose/QMMM/InteractionTermEliminator.h index 39327c5..fc0b820 100644 --- a/src/Swoose/Swoose/QMMM/InteractionTermEliminator.h +++ b/src/Swoose/Swoose/QMMM/InteractionTermEliminator.h @@ -8,6 +8,7 @@ #ifndef SWOOSE_QMMM_INTERACTIONTERMELIMINATOR_H #define SWOOSE_QMMM_INTERACTIONTERMELIMINATOR_H +#include #include #include #include @@ -22,10 +23,10 @@ class DihedralTerm; class ImproperDihedralTerm; class DispersionTerm; class RepulsionTerm; -class ElectrostaticTerm; class HydrogenBondTerm; -class LennardJonesTerm; class InteractionTermBase; +class LennardJonesEvaluator; +class ElectrostaticEvaluator; } // namespace MolecularMechanics namespace Qmmm { @@ -80,7 +81,7 @@ class InteractionTermEliminator { void eliminateRepulsionTerms(std::vector& repulsionTerms); // @brief Eliminates those Lennard-Jones interactions which are already covered by the QM calculation in QM/MM. - void eliminateLennardJonesTerms(std::vector& ljTerms); + void eliminateLennardJonesTerms(MolecularMechanics::LennardJonesEvaluator& lennardJonesEvaluator); // @brief Eliminates those hydrogen bond interactions which are already covered by the QM calculation in QM/MM. void eliminateHydrogenBondTerms(std::vector& hydrogenBondTerms); @@ -90,7 +91,7 @@ class InteractionTermEliminator { * @param electrostaticTerms The electrostatic terms of the MM model. * @param electrostaticEmbedding Whether electrostatic embedding is switched on in the QM/MM calculator. */ - void eliminateElectrostaticTerms(std::vector& electrostaticTerms, + void eliminateElectrostaticTerms(MolecularMechanics::ElectrostaticEvaluator& electrostaticEvaluator, bool electrostaticEmbedding); /* @@ -98,6 +99,7 @@ class InteractionTermEliminator { */ void eliminateTerm(MolecularMechanics::InteractionTermBase& term, const std::vector& atomsInTerm, int allowedInQmRegion); + bool termToEliminate(const std::vector& atomsInTerm, int allowedInQmRegion); /* * @brief Returns whether an atom with the given index is in the QM region. */ @@ -107,11 +109,11 @@ class InteractionTermEliminator { template void enableTerms(std::vector& terms); - // Eliminates all of the shared interaction terms that are present in all of the possible MM calculators. + // Eliminates all the shared interaction terms that are present in all the possible MM calculators. template void eliminateSharedInteractionTerms(bool electrostaticEmbedding, CalculatorType& calculator); - // Enables all of the shared interaction terms that are present in all of the possible MM calculators. + // Enables all of the shared interaction terms that are present in all the possible MM calculators. template void enableSharedInteractionTerms(CalculatorType& calculator); @@ -123,6 +125,9 @@ class InteractionTermEliminator { // Boolean to keep track of whether we are in a reduced QM/MM energy calculation. bool eliminateEnvironmentOnlyTerms_ = false; + // The exclusion set before eliminating QM-QM internal terms. + Eigen::SparseMatrix originalLJExclusions_; + Eigen::SparseMatrix originalElectrostaticExclusions_; }; } // namespace Qmmm diff --git a/src/Swoose/Swoose/QMMM/QmRegionSelection/QmRegionCandidateGenerator.h b/src/Swoose/Swoose/QMMM/QmRegionSelection/QmRegionCandidateGenerator.h index 9ceaac1..7685886 100644 --- a/src/Swoose/Swoose/QMMM/QmRegionSelection/QmRegionCandidateGenerator.h +++ b/src/Swoose/Swoose/QMMM/QmRegionSelection/QmRegionCandidateGenerator.h @@ -14,7 +14,7 @@ namespace Scine { namespace Core { -class Log; +struct Log; } // namespace Core namespace Utils { diff --git a/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmDirectCalculationsHelper.cpp b/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmDirectCalculationsHelper.cpp index 490287b..0cb1006 100644 --- a/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmDirectCalculationsHelper.cpp +++ b/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmDirectCalculationsHelper.cpp @@ -41,8 +41,8 @@ QmmmDirectCalculationsHelper::QmmmDirectCalculationsHelper(std::shared_ptr QmmmDirectCalculationsHelper::calculateForces() { // Initialize some variables - int numCandidateModels = qmmmModelCandidates_.size(); - int numModels = numCandidateModels + qmmmReferenceModels_.size(); + size_t numCandidateModels = qmmmModelCandidates_.size(); + size_t numModels = numCandidateModels + qmmmReferenceModels_.size(); std::vector forces(numModels); // Get max allowed symmetry score @@ -57,8 +57,9 @@ std::vector QmmmDirectCalculationsHelper::calculateForces() { } #pragma omp parallel for - for (long unsigned int i = 0; i < numModels; ++i) { - int molecularCharge, spinMultiplicity; + for (size_t i = 0; i < numModels; ++i) { + int molecularCharge = 0; + int spinMultiplicity = 1; if (i < numCandidateModels) { molecularCharge = qmmmModelCandidates_.at(i).molecularCharge; spinMultiplicity = qmmmModelCandidates_.at(i).spinMultiplicity; diff --git a/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmDirectCalculationsHelper.h b/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmDirectCalculationsHelper.h index a875492..11f357b 100644 --- a/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmDirectCalculationsHelper.h +++ b/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmDirectCalculationsHelper.h @@ -15,7 +15,7 @@ namespace Scine { namespace Core { -class Log; +struct Log; } // namespace Core namespace Utils { diff --git a/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmModelAnalyzer.h b/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmModelAnalyzer.h index 4130a3e..ea0b6b8 100644 --- a/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmModelAnalyzer.h +++ b/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmModelAnalyzer.h @@ -15,7 +15,7 @@ namespace Scine { namespace Core { -class Log; +struct Log; } // namespace Core namespace Utils { diff --git a/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmReferenceDataManager.h b/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmReferenceDataManager.h index 7df32ae..bb1f3d2 100644 --- a/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmReferenceDataManager.h +++ b/src/Swoose/Swoose/QMMM/QmRegionSelection/QmmmReferenceDataManager.h @@ -15,7 +15,7 @@ namespace Scine { namespace Core { -class Log; +struct Log; } // namespace Core namespace Utils { diff --git a/src/Swoose/Swoose/QMMM/QmmmHelpers.cpp b/src/Swoose/Swoose/QMMM/QmmmHelpers.cpp index b683454..c7757d2 100644 --- a/src/Swoose/Swoose/QMMM/QmmmHelpers.cpp +++ b/src/Swoose/Swoose/QMMM/QmmmHelpers.cpp @@ -70,13 +70,23 @@ void writePointChargesFile(const Utils::PositionCollection& positions, // First write the atomic charges for (int i = 0; i < positions.rows(); ++i) { // Write charge if this atom is NOT a QM atom. + double charge = chargeRedistributionResult.atomicCharges.at(i); + // TODO: This is a quick fix to ensure that no point charges are written which are + // so small that Turbomole doesn't consider them. A cleaner solution would be + // desirable. + // Note that we set the charge to 2e-6 such that all point charges are properly + // read in: our Turbomole output parser checks for a value of 1e-6 to determine + // the number of point charges. + if (writeTurbomoleFormat && std::fabs(charge) < 2e-6) { + charge = 2e-6; + } if (std::find(listOfQmAtoms.begin(), listOfQmAtoms.end(), i) == listOfQmAtoms.end()) { if (writeTurbomoleFormat) { pcFile << positions.row(i) << " "; - pcFile << chargeRedistributionResult.atomicCharges.at(i) << "\n"; + pcFile << charge << "\n"; } else { - pcFile << chargeRedistributionResult.atomicCharges.at(i) << " "; + pcFile << charge << " "; pcFile << positions.row(i) * Utils::Constants::angstrom_per_bohr << "\n"; } } diff --git a/src/Swoose/Swoose/StructurePreparation/Protonation/ProtonationHandler.cpp b/src/Swoose/Swoose/StructurePreparation/Protonation/ProtonationHandler.cpp index 413c048..04e49bb 100644 --- a/src/Swoose/Swoose/StructurePreparation/Protonation/ProtonationHandler.cpp +++ b/src/Swoose/Swoose/StructurePreparation/Protonation/ProtonationHandler.cpp @@ -298,20 +298,20 @@ void ProtonationHandler::protonateTetrahedralGroups(const Utils::AtomCollection& if (nNeighbors_[index] == 1) { auto firstBondedAtom = listsOfNeighbors_[index].begin(); std::vector hPos(3); // v2, v3, v4; - Eigen::Vector3d v1 = (structure.at(*firstBondedAtom).getPosition() - pos).normalized(); + const Eigen::Vector3d v1 = (structure.at(*firstBondedAtom).getPosition() - pos).normalized(); Utils::StructuralCompletion::generate3TetrahedronCornersFrom1Other(v1, hPos.at(0), hPos.at(1), hPos.at(2)); - int counter = 0; + size_t counter = 0; for (auto i : hPos) { if (counter == hPos.size() - 1 && isPseudoTetrahedralAtom) { break; } auto bonded_element = structure.at(index).getElementType(); - double sumOfVdWRadii = Utils::ElementInfo::covalentRadius(bonded_element) + - Utils::ElementInfo::covalentRadius(Utils::ElementType::H); + const double sumOfVdWRadii = Utils::ElementInfo::covalentRadius(bonded_element) + + Utils::ElementInfo::covalentRadius(Utils::ElementType::H); i *= sumOfVdWRadii; i += pos; - Utils::Atom hydrogen(Utils::ElementType::H, i); + const Utils::Atom hydrogen(Utils::ElementType::H, i); hydrogenAtoms_.push_back(hydrogen); counter++; } @@ -321,21 +321,22 @@ void ProtonationHandler::protonateTetrahedralGroups(const Utils::AtomCollection& auto secondBondedAtom = std::next(listsOfNeighbors_[index].begin()); std::vector hPos(2); - Eigen::Vector3d v1 = (structure.at(*firstBondedAtom).getPosition() - pos).normalized(); - Eigen::Vector3d v2 = (structure.at(*secondBondedAtom).getPosition() - pos).normalized(); + const Eigen::Vector3d v1 = (structure.at(*firstBondedAtom).getPosition() - pos).normalized(); + const Eigen::Vector3d v2 = (structure.at(*secondBondedAtom).getPosition() - pos).normalized(); Utils::StructuralCompletion::generate2TetrahedronCornersFrom2Others(v1, v2, hPos.at(0), hPos.at(1)); - int counter = 0; + // TODO: Is there a reason why this counter is never incremented? + size_t counter = 0; for (auto i : hPos) { if (counter == hPos.size() - 1 && isPseudoTetrahedralAtom) { break; } auto bonded_element = structure.at(index).getElementType(); - double sumOfVdWRadii = Utils::ElementInfo::covalentRadius(bonded_element) + - Utils::ElementInfo::covalentRadius(Utils::ElementType::H); + const double sumOfVdWRadii = Utils::ElementInfo::covalentRadius(bonded_element) + + Utils::ElementInfo::covalentRadius(Utils::ElementType::H); i *= sumOfVdWRadii; i += pos; - Utils::Atom hydrogen(Utils::ElementType::H, i); + const Utils::Atom hydrogen(Utils::ElementType::H, i); hydrogenAtoms_.push_back(hydrogen); } diff --git a/src/Swoose/Swoose/StructurePreparation/SpecialCaseHandler.h b/src/Swoose/Swoose/StructurePreparation/SpecialCaseHandler.h index 1ef3401..a5844a7 100644 --- a/src/Swoose/Swoose/StructurePreparation/SpecialCaseHandler.h +++ b/src/Swoose/Swoose/StructurePreparation/SpecialCaseHandler.h @@ -13,7 +13,7 @@ namespace Scine { namespace Core { -class Log; +struct Log; } namespace Utils { @@ -22,7 +22,7 @@ enum class ElementType : unsigned; } // namespace Utils namespace StructurePreparation { -class StructurePreparationData; +struct StructurePreparationData; namespace SpecialCaseHandler { @@ -59,4 +59,4 @@ bool isCarboxylateC(const StructurePreparationData& data, int index); } // namespace StructurePreparation } // namespace Scine -#endif // PDBPREPARATION_SPECIALCASEHANDLER_H \ No newline at end of file +#endif // PDBPREPARATION_SPECIALCASEHANDLER_H diff --git a/src/Swoose/Swoose/StructurePreparation/StructurePreparationHelper.cpp b/src/Swoose/Swoose/StructurePreparation/StructurePreparationHelper.cpp index a3abf92..b57102d 100644 --- a/src/Swoose/Swoose/StructurePreparation/StructurePreparationHelper.cpp +++ b/src/Swoose/Swoose/StructurePreparation/StructurePreparationHelper.cpp @@ -56,6 +56,7 @@ std::string indexToAtomType(StructurePreparationData& data, int index) { if (proteinAtom.index == index) return proteinAtom.atomType; } + throw std::runtime_error("No atom type available for the given index!"); } void getSideChainNeighbors(StructurePreparationData& data, int atomIndex, std::list& atomsToAdd) { @@ -272,7 +273,7 @@ void mapSubsystemIndicesToFullStructure(const Utils::AtomCollection& fullStructu subsystemMapping.push_back(indicesInStructure); } -void handleBoundariesBetweenProteinAndNonRegContainer(StructurePreparationData& data, StructurePreparationFiles& files) { +void handleBoundariesBetweenProteinAndNonRegContainer(StructurePreparationData& data) { // detect boundary regions std::map boundaryAtoms; for (int i = 0; i < data.fullStructure.size(); ++i) { @@ -467,8 +468,8 @@ std::vector collectTitrableSites(StructurePreparationData& data) { void determineChargedSites(std::vector& listOfNegatives, std::vector& listOfPositives, const StructurePreparationData& data) { for (int index = 0; index < data.numberOfAtoms; ++index) { - bool negative = SpecialCaseHandler::isNegative(data, index, listOfNegatives); - bool positive = SpecialCaseHandler::isPositive(data, index, listOfPositives); + SpecialCaseHandler::isNegative(data, index, listOfNegatives); + SpecialCaseHandler::isPositive(data, index, listOfPositives); } } diff --git a/src/Swoose/Swoose/StructurePreparation/StructurePreparationHelper.h b/src/Swoose/Swoose/StructurePreparation/StructurePreparationHelper.h index b184a9f..ee8b172 100644 --- a/src/Swoose/Swoose/StructurePreparation/StructurePreparationHelper.h +++ b/src/Swoose/Swoose/StructurePreparation/StructurePreparationHelper.h @@ -56,7 +56,7 @@ void reevaluateConnectivityForAminoAcids(StructurePreparationData& data); * @brief After combination of (modified) substructures, the covalent bonds at the subsystem boundaries must be fixed * (superfluous hydrogens must be removed). */ -void handleBoundariesBetweenProteinAndNonRegContainer(StructurePreparationData& data, StructurePreparationFiles& files); +void handleBoundariesBetweenProteinAndNonRegContainer(StructurePreparationData& data); /** * @brief Solvates the structure with a number of solvent shells that can be set as a setting. */ @@ -110,4 +110,4 @@ bool compareByIndex(const ProteinAtom& a, const ProteinAtom& b); } // namespace StructurePreparation } // namespace Scine -#endif // PDBPREPARATION_PDBPREPARATIONHELPER_H \ No newline at end of file +#endif // PDBPREPARATION_PDBPREPARATIONHELPER_H diff --git a/src/Swoose/Swoose/StructurePreparation/StructurePreparationIO.cpp b/src/Swoose/Swoose/StructurePreparation/StructurePreparationIO.cpp index a3f2275..93cfcc7 100644 --- a/src/Swoose/Swoose/StructurePreparation/StructurePreparationIO.cpp +++ b/src/Swoose/Swoose/StructurePreparation/StructurePreparationIO.cpp @@ -25,7 +25,7 @@ namespace StructurePreparationIO { void xyzToPdb(const std::string& xyzFile, const std::string& pdbFile) { Utils::AtomCollection at; at = Utils::ChemicalFileHandler::read(xyzFile).first; - Utils::ChemicalFileHandler::write(pdbFile, at, "Dummy comment"); + Utils::ChemicalFileHandler::write(pdbFile, at); } void writeAtomicInfoFileForProtein(const StructurePreparationData& data, const std::string& atomicInfoFile) { @@ -79,6 +79,7 @@ std::string getSuffix(const bfs::path& filepath) { return suffix.substr(1); } +// TODO: This is already supported through the utils and could be removed in the future. void writePdbFileWithResidueSpecifier(const StructurePreparationData& data, const std::string& proteinFile, Core::Log& log) { std::ofstream pdbFile; if (!data.protein.empty()) { diff --git a/src/Swoose/Swoose/StructurePreparation/StructureProcessor.cpp b/src/Swoose/Swoose/StructurePreparation/StructureProcessor.cpp index ae06631..c9e9899 100644 --- a/src/Swoose/Swoose/StructurePreparation/StructureProcessor.cpp +++ b/src/Swoose/Swoose/StructurePreparation/StructureProcessor.cpp @@ -71,11 +71,13 @@ void StructureProcessor::analyzeStructure(const std::string& structureFile) { log_.output << "Analyzing the input structure ..." << Core::Log::nl; StructurePreparationHelper::performInitialCheck(files_, 1, settings_); auto structure = getAtomCollectionFromInput(structureFile, false); + structure.removeAtomsByResidueLabel({"HOH"}); StructurePreparationHelper::updatePdbPreparationData(data_, structure); StructurePreparationHelper::performGraphAnalysisOnStructure(data_, data_.fullStructure); StructurePreparationHelper::updateNonRegContainerVector(data_); StructurePreparationHelper::reevaluateConnectivityForAminoAcids(data_); - if ((data_.vectorOfNonRegContainerIndices.size() + data_.vectorOfProteinIndices.size()) != data_.fullStructure.size()) { + if (size_t(data_.vectorOfNonRegContainerIndices.size() + data_.vectorOfProteinIndices.size()) != + size_t(data_.fullStructure.size())) { log_.output << "Structure analysis failed!" << Core::Log::nl; } StructurePreparationIO::writePdbFileWithResidueSpecifier(data_, files_.proteinFile, log_); @@ -107,7 +109,7 @@ void StructureProcessor::finalize() { // merge substructures StructurePreparationHelper::mergeProteinAndNonRegContainer(data_, files_); // handle the boundaries - StructurePreparationHelper::handleBoundariesBetweenProteinAndNonRegContainer(data_, files_); + StructurePreparationHelper::handleBoundariesBetweenProteinAndNonRegContainer(data_); // write the atomic info file for the protein first StructurePreparationIO::writeAtomicInfoFileForProtein(data_, files_.atomicInfoFile); StructurePreparationIO::addAtomicInformationForNonRegContainer(files_, data_.subsystemMapping); @@ -164,7 +166,6 @@ Utils::AtomCollection StructureProcessor::getAtomCollectionFromInput(const std:: Utils::PdbStreamHandler handler; handler.setReadH(includeH); // first read in all solvent molecules - handler.parseOnlySolvent(false); std::ifstream input; input.open(structureFile); auto structures = handler.read(input); @@ -174,8 +175,10 @@ Utils::AtomCollection StructureProcessor::getAtomCollectionFromInput(const std:: std::vector solvents; try { - handler.parseOnlySolvent(true); solvents = handler.read(input); + for (auto& solvent : solvents) { + solvent.keepAtomsByResidueLabel({"HOH"}); + } } catch (...) { solvents.resize(0); diff --git a/src/Swoose/Swoose/Utilities/SettingsPopulator.h b/src/Swoose/Swoose/Utilities/SettingsPopulator.h index 89a95c9..ba409c4 100644 --- a/src/Swoose/Swoose/Utilities/SettingsPopulator.h +++ b/src/Swoose/Swoose/Utilities/SettingsPopulator.h @@ -88,6 +88,7 @@ class SettingsPopulator { static void addMaxNumRefModelsOption(Utils::UniversalSettings::DescriptorCollection& settings); static void addTolerancesForQmRegionSelection(Utils::UniversalSettings::DescriptorCollection& settings); static void addQmRegionSelectionRandomSeed(Utils::UniversalSettings::DescriptorCollection& settings); + static void addPerformSfamParametrizations(Utils::UniversalSettings::DescriptorCollection& settings); // These functions used for QM region selection via a Database static void addMethodFamily(Utils::UniversalSettings::DescriptorCollection& settings); diff --git a/src/Swoose/Swoose/Utilities/SubsystemGenerator.h b/src/Swoose/Swoose/Utilities/SubsystemGenerator.h index 6ae1cc4..aef0b6f 100644 --- a/src/Swoose/Swoose/Utilities/SubsystemGenerator.h +++ b/src/Swoose/Swoose/Utilities/SubsystemGenerator.h @@ -16,7 +16,7 @@ namespace Scine { namespace Core { -class Log; +struct Log; } // namespace Core namespace Utils { diff --git a/src/Swoose/Tests/Files/tests_file_location.h.in b/src/Swoose/Tests/Files/tests_file_location.h.in index 12f54b1..eca8098 100644 --- a/src/Swoose/Tests/Files/tests_file_location.h.in +++ b/src/Swoose/Tests/Files/tests_file_location.h.in @@ -41,4 +41,4 @@ const char* const gaff_parameter_file = "@Tests_dir@gaff.dat"; const char* const gaff_dummy_atom_types_file = "@Tests_dir@gaff_dummy_atom_types.dat"; const char* const qm_region_selection_yaml_file = "@Tests_dir@settings_qm_region_selection_direct.yaml"; -#endif // TESTS_FILE_LOCATION_H \ No newline at end of file +#endif // TESTS_FILE_LOCATION_H diff --git a/src/Swoose/Tests/StructurePreparationTest.cpp b/src/Swoose/Tests/StructurePreparationTest.cpp index cd466eb..7eee061 100644 --- a/src/Swoose/Tests/StructurePreparationTest.cpp +++ b/src/Swoose/Tests/StructurePreparationTest.cpp @@ -79,8 +79,8 @@ TEST_F(AStructurePreparationTest, ProteinAndNonRegContainerAreSeparatedCorrectly processor.setFiles(testFiles); processor.analyzeStructure(plastocyanin_pdb_file); - Utils::AtomCollection nonRegContainer = Utils::ChemicalFileHandler::read(testFiles.nonRegContainerFile).first; - Utils::PdbStreamHandler handler; + const Utils::AtomCollection nonRegContainer = Utils::ChemicalFileHandler::read(testFiles.nonRegContainerFile).first; + const Utils::PdbStreamHandler handler; std::ifstream input; input.open(testFiles.proteinFile); auto protein = handler.read(input); @@ -184,7 +184,7 @@ TEST_F(AStructurePreparationTest, CustomInformationForNonRegContainerCanBeAdded) Utils::NativeFilenames::combinePathSegments(structure_preparation_dir, systemName, "system.xyz")); StructurePreparationHelper::mergeProteinAndNonRegContainer(data, testFiles); - StructurePreparationHelper::handleBoundariesBetweenProteinAndNonRegContainer(data, testFiles); + StructurePreparationHelper::handleBoundariesBetweenProteinAndNonRegContainer(data); // write the atomic info file for the protein first StructurePreparationIO::writeAtomicInfoFileForProtein(data, testFiles.atomicInfoFile); StructurePreparationIO::addAtomicInformationForNonRegContainer(testFiles, data.subsystemMapping); diff --git a/src/Swoose/Tests/TestUtilities/MockQmCalculator.h b/src/Swoose/Tests/TestUtilities/MockQmCalculator.h index d9ebd91..3da67fa 100644 --- a/src/Swoose/Tests/TestUtilities/MockQmCalculator.h +++ b/src/Swoose/Tests/TestUtilities/MockQmCalculator.h @@ -53,7 +53,7 @@ class MockQmCalculatorSettings : public Scine::Utils::Settings { }; /// @brief Mock class for a QM calculator to use in QM/MM tests. -class MockQmCalculator : public Utils::CloneInterface { +class MockQmCalculator final : public Utils::CloneInterface { public: static constexpr const char* model = "MOCK-QM"; /// @brief Default Constructor. diff --git a/src/Swoose/config.cmake.in b/src/Swoose/config.cmake.in index 6673b83..b8b5642 100644 --- a/src/Swoose/config.cmake.in +++ b/src/Swoose/config.cmake.in @@ -9,7 +9,7 @@ find_dependency(ScineMolassembler REQUIRED) get_filename_component(SWOOSE_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" PATH) # Generate Target -include("${CMAKE_CURRENT_LIST_DIR}/SwooseTargets.cmake") +include("${CMAKE_CURRENT_LIST_DIR}/ScineSwooseTargets.cmake") @PACKAGE_INIT@