From 73b0b26ad19bc8bfea94b89e1c3aa35e1925bd8c Mon Sep 17 00:00:00 2001 From: Reiher Research Group Date: Thu, 23 Jun 2022 09:03:39 +0200 Subject: [PATCH] Release 4.0.0 --- .gitignore | 7 + CHANGELOG.rst | 9 ++ CMakeLists.txt | 10 +- README.rst | 6 +- conanfile.py | 8 +- dev | 2 +- manual/readuct_manual.tex | 152 ++++++++++++++---- src/Readuct/App/Tasks/AfirOptimizationTask.h | 4 + .../App/Tasks/BSplineInterpolationTask.h | 5 +- src/Readuct/App/Tasks/BondOrderTask.h | 7 +- .../App/Tasks/GeometryOptimizationTask.h | 111 ++++++++++--- src/Readuct/App/Tasks/HessianTask.h | 2 + src/Readuct/App/Tasks/IrcTask.h | 6 + src/Readuct/App/Tasks/NtOptimization2Task.h | 147 +++++++++++++++++ src/Readuct/App/Tasks/NtOptimizationTask.h | 27 ++-- src/Readuct/App/Tasks/SinglePointTask.h | 77 ++++++++- src/Readuct/App/Tasks/Task.h | 11 +- src/Readuct/App/Tasks/TaskFactory.h | 7 +- src/Readuct/App/Tasks/TsOptimizationTask.h | 11 ++ src/Readuct/App/main.cpp | 10 +- src/Readuct/CMakeLists.txt | 2 +- src/Readuct/Files.cmake | 2 + src/Readuct/Python/TasksPython.cpp | 4 + src/Readuct/Python/sphinx/index.rst | 3 +- .../ElementaryStepOptimizer.h | 6 +- src/Readuct/Tests/AppTests/test_readuct.py | 111 ++++++++++--- 26 files changed, 616 insertions(+), 131 deletions(-) create mode 100644 src/Readuct/App/Tasks/NtOptimization2Task.h diff --git a/.gitignore b/.gitignore index b534b4a..6c687d6 100644 --- a/.gitignore +++ b/.gitignore @@ -14,6 +14,8 @@ tmp_clang_* compile_commands.json .lvimrc *.vim +.DS_Store +*~ # manual output manual/readuct_manual.aux @@ -30,3 +32,8 @@ manual/readuct_manual.ilg # conan ci files conan/* + +# python +*/__pycache__/ +*.pyc +.pytest_cache/ diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 4ea0101..975a2c9 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,6 +1,15 @@ Changelog ========= +Release 4.0.0 +------------- + +- Add 2nd Newton trajectory scan algorithm (NT2) +- Update automatic TS mode picking to respect frequencies to some degree +- Deprecate BondOrderTask and allow bond order calculation in SinglePointTask +- Add option for spin propensity check +- Add option to optimize periodic boundaries in geometry optimization + Release 3.0.0 ------------- diff --git a/CMakeLists.txt b/CMakeLists.txt index 936d94a..6a836dc 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(Readuct - VERSION 3.0.0 + VERSION 4.0.0 DESCRIPTION "This is the SCINE module Readuct." ) @@ -21,12 +21,18 @@ if(SCINE_BUILD_TESTS) endif() # Enable testing -option(BUILD_SPARROW "Will download and build Sparrow (the Scine semi-empirical module)." OFF) +option(BUILD_SPARROW "Will download and build Sparrow (the SCINE semi-empirical module)." OFF) if(BUILD_SPARROW) include(ImportSparrow) import_sparrow() endif() +option(BUILD_XTB "Will download and build SCINE Xtb (the SCINE Wrapper around the xtb API by the Grimme group)." OFF) +if(BUILD_XTB) + include(ImportXtb) + import_xtb() +endif() + option(SCINE_USE_MKL "Use the optimized MKL library for linear algebra operations of Eigen" OFF) # Subdirectories diff --git a/README.rst b/README.rst index 564a5a7..3b8b15a 100644 --- a/README.rst +++ b/README.rst @@ -18,7 +18,7 @@ SCINE ReaDuct is a command-line tool that allows you to carry out For these calculations, it relies on a backend program to provide the necessary quantum chemical properties (such as nuclear gradients). Currently, SCINE Sparrow, -SCINE XTB, CP2K, Gaussian, ORCA, and Turbomole are supported as backend programs. +XTB, CP2K, Gaussian, ORCA, and Turbomole are supported as backend programs. License and Copyright Information --------------------------------- @@ -33,8 +33,8 @@ Installation and Usage For instructions on how to install and use ReaDuct as well as for a detailed documentation of the entire functionality of ReaDuct, please consult the user manual found in the ``manual`` directory in in the repository. -Alternatively the manual can also be found on the official GitHub website, -on the SCINE website and in the hosted documentation. +Alternatively the manual can also be found on the official GitHub website +and on the SCINE website. How to Cite ----------- diff --git a/conanfile.py b/conanfile.py index 0af057c..56f222e 100644 --- a/conanfile.py +++ b/conanfile.py @@ -3,7 +3,7 @@ class ScineReaductConan(ScineConan): name = "scine_readuct" - version = "3.0.0" + version = "4.0.0" url = "https://github.com/qcscine/readuct" description = """ SCINE ReaDuct is a command-line tool that allows to carry out: @@ -18,7 +18,7 @@ class ScineReaductConan(ScineConan): - Newton trajectory scans searching for transition state guesses. For these calculations, it relies on a backend program to provide the necessary quantum chemical properties (such as nuclear gradients). Currently, SCINE Sparrow -SCINE XTB, CP2K, Gaussian, ORCA, and Turbomole are supported as backend programs.""" +XTB, CP2K, Gaussian, ORCA, and Turbomole are supported as backend programs.""" options = { "shared": [True, False], "python": [True, False], @@ -38,7 +38,7 @@ class ScineReaductConan(ScineConan): "dev/cmake/*", "src/*", "CMakeLists.txt", "README.rst", "LICENSE.txt", "dev/conan/hook.cmake", "dev/conan/glue/*" ] - requires = ["scine_utilities/[~=4.0.0]", + requires = ["scine_utilities/[=5.0.0]", "boost/[>1.65.0]", "yaml-cpp/0.6.3"] cmake_name = "Readuct" @@ -55,6 +55,6 @@ def configure(self): def build_requirements(self): if self.options.tests: - self.build_requires("scine_sparrow/[~=3.0.0]") + self.build_requires("scine_sparrow/[=3.0.1]") super().build_requirements() diff --git a/dev b/dev index 413ce8a..abb0acc 160000 --- a/dev +++ b/dev @@ -1 +1 @@ -Subproject commit 413ce8aa2ccd1cdcdcd49366f493c90627fad5c7 +Subproject commit abb0acc9e1bb05da2a8dd4c5804d415b7f596e1a diff --git a/manual/readuct_manual.tex b/manual/readuct_manual.tex index 923e994..fc69a0e 100644 --- a/manual/readuct_manual.tex +++ b/manual/readuct_manual.tex @@ -1,5 +1,8 @@ \documentclass[]{tufte-book} +\titleclass{\subsubsection}{straight} +\titleformat{\subsubsection}[hang]{\normalfont\itshape}{\thesubsubsection}{1em}{}[] + \hypersetup{colorlinks}% uncomment this line if you prefer colored hyperlinks (e.g., for onscreen viewing) %% @@ -12,7 +15,7 @@ %% % Book metadata -\title[SCINE ReaDuct manual]{User Manual \vskip 0.5em {\setlength{\parindent}{0pt} \Huge SCINE ReaDuct 3.0.0}} +\title[SCINE ReaDuct manual]{User Manual \vskip 0.5em {\setlength{\parindent}{0pt} \Huge SCINE ReaDuct 4.0.0}} \author[The SCINE ReaDuct Developers]{The SCINE ReaDuct Developers: \newline \noindent Christoph Brunken, Katja-Sophia Csizi, Stephanie Grimmel, Stefan Gugler, Jan-Grimo Sobez, Miguel Steiner, Paul T\"urtscher, Jan Unsleber, Alain C. Vaucher, Thomas Weymuth, and Markus Reiher} \publisher{ETH Z\"urich} @@ -253,7 +256,7 @@ \chapter{Installation}\label{ch:installation} \textsc{ReaDuct} is distributed as an open source code. In order to compile \textsc{ReaDuct} from this source code, you need \begin{itemize} - \item a C++ compiler supporting the C++14 standard (GCC at least 7.3.0), + \item a C++ compiler supporting the C++17 standard (GCC at least 7.3.0), \item cmake (at least 3.9.0), \item the Boost library (at least 1.65.0), and \item the Eigen3 library (at least 3.3.2). @@ -371,7 +374,7 @@ \subsection{SCINE \textsc{Sparrow}} \subsection{External Programs} -SCINE \textsc{ReaDuct} supports calculations with \textsc{ORCA}\cite{orca} (version 4.2.0), \textsc{Turbomole}\cite{turbomole} (version 7.4.1), \textsc{CP2K} (8.1)\cite{cp2k}, \textsc{Gaussian}\cite{gaussian09} (version 09 Rev. D01), and xtb (6.4.1)\cite{xtb}. Support is currently not fully tested. For each program, there might be specific calculation types and/or settings which do not work. Also, we cannot guarantee compatibility with any version different from the one's mentioned above since we have no control over the output format of an external program. If you encounter any problems when +SCINE \textsc{ReaDuct} supports calculations with \textsc{ORCA}\cite{orca} (version 4.2.0), \textsc{Turbomole}\cite{turbomole} (version 7.4.1), \textsc{CP2K} (8.1)\cite{cp2k}, \textsc{Gaussian}\cite{gaussian09} (version 09 Rev. D01), and xtb (6.4.1)\cite{xtb}. Support is currently not fully tested. For each program, there might be specific calculation types and/or settings which do not work. Also, we cannot guarantee compatibility with any version different from the one's mentioned above since we have no control over the output format of an external program. If you encounter any problems when using one of these software packages together with \textsc{ReaDuct}, please write a short message to \href{scine@phys.chem.ethz.ch}{scine@phys.chem.ethz.ch}. For these programs, the following settings can be applied. @@ -391,9 +394,10 @@ \subsection{External Programs} hartree). The default is \texttt{1.0e-7} (\textit{i.e.}, 10\textsuperscript{$-$7}\,hartree) for all programs. \item \texttt{max\_scf\_iterations}: The maximum number of SCF iterations allowed by the program. By default, it is 100. Note that this keyword is currently not available in conjunction with \textsc{Gaussian}. \item \texttt{electronic\_temperature}: The electronic temperature used for the calculation. By default, it is 0\,K, except for xtb, in which case it is 300\,K. +\item \texttt{temperature}: The temperature used for the calculation of thermochemical data. By default, it is 298.15\,K. Note that the calculation of thermochemistry data is currently not supported for \textsc{Gaussian}. \end{itemize} -Additional program-specific settings are compiled in the subsequent sections. +A form of dispersion correction must be specified within the method following a hypen, e.g. \texttt{PBE-D3BJ}. The supported dispersion correction between programs may vary. Additional program-specific settings are compiled in the subsequent sections. \subsection{\textsc{ORCA}} @@ -427,8 +431,8 @@ \subsection{\textsc{ORCA}} in no point charges considered. \item \texttt{delete\_tmp\_files}: Whether temporary files (i.e., all files with a ".tmp" extension) should be deleted in case the calculation fails. By default, this is set to \texttt{true}. -\item \texttt{temperature}: The temperature used for the calculation of thermochemical data. By default, it is 298.15\,K. \item \texttt{hessian\_calculation\_type}: The way to calculate a Hessian, i.e., either ``analytical`` or ``numerical``. By default, it is set to ``analytical``, unless it is already specified in our code that \textsc{ORCA} does not provide an analytical Hessian for the specified method. However, this list is not complete and may vary from version to version. +\item \texttt{special\_option}: A possibility to add an input that is not covered by our settings, only recommended for experts. \end{itemize} \subsection{\textsc{Turbomole}} @@ -446,22 +450,29 @@ \subsection{\textsc{Turbomole}} You can specify the following settings in the settings block: \begin{itemize} \item \texttt{method}: -This specifies the calculation method to be used. By default, it is \texttt{'PBE'}. Dispersion functional can be added to this method string, e.g: \texttt{'PBE D3BJ'}. Current supported dispersion corrections are D3\cite{grimmeD3} and D3BJ.\cite{bjDamping} +This specifies the calculation method to be used. By default, it is \texttt{'PBE'}. Dispersion correction can be added to this method string, e.g: \texttt{'PBE-D3BJ'}. Current supported dispersion corrections are D3\cite{grimmeD3} and D3BJ.\cite{bjDamping} \item \texttt{basis\_set}: This specifies the basis set string. By default, Turbomole assigns the \texttt{'def-SV(P)'} basis set. You can specify any valid \textsc{Turbomole} basis set string from the Karlsruhe basis sets\cite{karlsruheBasisSets}, Dunning's correlation-consistent basis sets\cite{dunningBasisSets} or Pople-style basis sets\cite{popleBasisSets} (see the \textsc{Turbomole} manual for a complete list). \item \texttt{solvation}: This specifies the implicit solvation model. By default, it is empty. Currently, you can specify \texttt{cosmo}. \item \texttt{solvent}: This specifies the solvent for implicit solvation. By default, it is empty. Current supported solvents are: \texttt{water}, \texttt{acetone}, \texttt{benzene} , \texttt{dmso}, \texttt{ethanol}, \texttt{methanol}, \texttt{hexane}, \texttt{toluene}, \texttt{ammonia}, \texttt{chloroform}, \texttt{nitrobenzene}, and \texttt{thf}. -\item \texttt{scf\_damping}: This specifies whether SCF damping is switched on to aid SCF convergence. By default, it is set to \texttt{false}. -\item \texttt{scf\_orbitalshift}: Shift all closed shells to lower energies to aid convergence. The default is set to \texttt{0.1}\,hartree. +\item \texttt{scf\_damping}: This specifies the SCF damping to be applied to aid SCF convergence. Available settings are \texttt{default} (default value of \textsc{Turbomole}), \texttt{low}, \texttt{medium}, and \texttt{high}. Note that the default damping values are lower than any of the other settings. +\item \texttt{scf\_orbitalshift}: Shifts either the energies of unoccupied MOs to higher energies or, in open-shell cases, the energies of closed-shell MOs to lower energies to aid convergence. The default is set to \texttt{0.1}\,hartree. \item \texttt{steer\_orbitals}: This specifies whether the orbitals should be perturbed following a randomized scheme after a converged single point calculation\cite{orbitalperturbation}. By default, this is set to \texttt{false}. -\item \texttt{temperature}: The temperature used for the calculation of thermochemical data. By default, it is 298.15\,K. \end{itemize} \subsection{\textsc{CP2K}} In order to use \textsc{CP2K} with \textsc{ReaDuct}, specify \texttt{program: 'CP2K'} in the respective system block and the desired calculation method family (e.g., \texttt{'DFT'}) in the \texttt{method\_family} key. -The method (e.g., \texttt{'PBE'}) can be specified in the system settings. +The method (e.g., \texttt{'PBE'}) can be specified in the system settings. Dispersion correction can be added to this method string, e.g: \texttt{'PBE-D3BJ'}. Current supported dispersion corrections are: +\begin{itemize} + \item D2 + \item D3 + \item D3BJ + \item DRSLL + \item LMKLL + \item RVV10 +\end{itemize} Implicit solvation is currently not supported. @@ -477,15 +488,6 @@ \subsection{\textsc{CP2K}} \item \texttt{plane\_wave\_cutoff}: Sets the plane wave cutoff of the finest grid in Ry. The default is 300. \item \texttt{relative\_multi\_grid\_cutoff}: Determines the grid at which a Gaussian is mapped, giving the cutoff in Ry used for a Gaussian with exponent (denoted alpha in the \textsc{CP2K} manual) of one. The default is 60. \item \texttt{n\_grids}: Sets the desired number of grids. The default is set to 5. - \item \texttt{vdw\_functional}: Sets the dispersion correction. By default, dispersion correction is disabled. The possible options are: - \begin{itemize} - \item DFTD2 - \item DFTD3 - \item DFTD3(BJ) - \item DRSLL - \item LMKLL - \item RVV10 - \end{itemize} \item \texttt{additional\_mos}: Specify the number of additional molecular orbitals. By default, it is set to zero. \item \texttt{orbital\_transformation}: Specify an orbital transformation minimizer, such as ``cg``. By default, it is set to ``none``, which deactivates orbital transformation. The available options are: \begin{itemize} @@ -530,14 +532,14 @@ \subsection{\textsc{CP2K}} By default it is set to ``additional\_output'', which results in a generated file named ``additional\_output-1.0.Log''. If this setting is identical to ``cp2k\_filename\_base``, no separate files will be generated. \item \texttt{delete\_tmp\_files}: Whether temporary files (i.e., all files with a ".bak" extension) should be deleted in case the calculation fails. By default, this is set to \texttt{true}. - \item \texttt{temperature}: The temperature used for the calculation of thermochemical data. By default, it is 298.15\,K. + \item \texttt{electronic\_temperature}: The temperature used for Fermi smearing. By default, it is 0.0\,K and therefore deactivated. Setting it to a value also requires to set the number of \texttt{additional\_mos}. \end{itemize} \subsection{\textsc{Gaussian}} In order to use \textsc{Gaussian} with \textsc{ReaDuct}, specify \texttt{program: 'GAUSSIAN'} in the respective system block and the desired calculation method family (e.g., \texttt{'DFT'}) in the \texttt{method\_family} key. -The method (e.g., \texttt{'PBEPBE'}) can be specified in the system settings. +The method (e.g., \texttt{'PBEPBE'}) can be specified in the system settings. Dispersion correction can be added to this method string, e.g: \texttt{'PBE-D3BJ'}. Current supported dispersion corrections are D2, D3\cite{grimmeD3} and D3BJ.\cite{bjDamping} Implicit solvation can be activated specifying the \texttt{solvent} key. To model implicit solvation, the Conductor-like Polarizable Continuum Model\cite{cpcm} (C-PCM) with the default solvent dependent settings implemented in Gaussian is used. @@ -568,7 +570,7 @@ \subsection{\textsc{xtb}} In order to use \textsc{xtb} with \textsc{ReaDuct}, specify \texttt{program: 'XTB'} in the respective system block and the desired calculation method family (e.g., \texttt{'GFN2'}) in the \texttt{method\_family} key. -You have to separately install the SCINE wrapper for xtb and make sure that the corresponding module is available in \texttt{SCINE\_MODULE\_PATH}. +It can be downloaded and integrated into \textsc{ReaDuct} automatically at compilation time with the CMake option \texttt{-DBUILD\_XTB=ON} (\textit{cf.}~also section~\nameref{ch:installation}). You can specify the following settings in the settings block: \begin{itemize} @@ -577,7 +579,6 @@ \subsection{\textsc{xtb}} \texttt{benzene}, \texttt{ch2cl2}, \texttt{chcl3}, \texttt{cs2}, \texttt{dmso}, \texttt{ether}, \texttt{methanol}, \texttt{toluene}, \texttt{thf}, \texttt{water}, \texttt{h2o}. For GFN-FF, they are \texttt{acetone}, \texttt{acetonitrile}, \texttt{benzene}, \texttt{ch2cl2}, \texttt{chcl3}, \texttt{cs2}, \texttt{dmf}, \texttt{dmso}, \texttt{ether}, \texttt{toluene}, \texttt{thf}, \texttt{water}, \texttt{h2o} -\item \texttt{temperature}: This specifies the temperature used for thermochemical calculations. The default is 298.15\,K. \item \texttt{symmetry\_number}: This specifies the symmetry number used thermochemical calculations. The default is 1. \item \texttt{print\_level}: This specifies the verbosity of the output. Possible values are 0, 1, and 2. By default, it is set to 0 (least verbose). \end{itemize} @@ -592,17 +593,25 @@ \subsection{Single Point Calculation} \texttt{type: 'sp'}, or \texttt{type: 'energy'}. The single point task will print partial atomic charges if the given model provides any. If charges are not required the keyword \texttt{'require\_charges': false} can be given in the tasks settings. In this case the program will not abort if a method that does not provide charges is requested. -Furthermore, the gradients with respect to the nuclear coordinates can be requested by setting the keyword \texttt{'require\_gradients': true} +Furthermore, the gradients with respect to the nuclear coordinates can be requested by setting the keyword \texttt{'require\_gradients': true}, bond orders can be requested by setting the keyword \texttt{'require\_bond\_orders': true}, and the orbital energies can be requested by setting the keyword \texttt{'orbital\_energies': true}. +This task also allows to carry out test calculations with other spin multiplicities to ensure that the ground state is calculated. +This option can be activated by setting \texttt{spin\_propensity\_check} to an integer value, which specifies the range to check above and below the current spin multiplicity. +For example, setting it to 1 will check a multiplicity of 1 and 5 for a system with a spin multiplicity of 3. +This setting is defaulted to 0, which deactivates the check.\\ Like in any task, the setting \texttt{stop\_on\_error} can be given to control whether the program throws an exception for a failed energy calculation or simply returns \texttt{false} for the task and proceeds with the remaining tasks. The default value is \texttt{true}. +Like in any task, the setting \texttt{silent\_stdout\_calculator} can be given to control whether any standard output of the calculation itself should be printed. The default value is \texttt{true}. \subsection{Bond Order Analysis} +This task has been deprecated; it will be removed in future versions. Please use the Single Point Task instead. + Depending on the chosen method it is possible to generate Mayer bond orders for a given system. In order to carry out this task, specify any of the following in the respective task block: \texttt{type: 'bond\_orders'}, \texttt{type: 'bondorders'}, \texttt{type: 'bonds'}, \texttt{type: 'bos'}, or \texttt{type: 'bo'}. This task also generates and states the electronic energy. Like in any task, the setting \texttt{stop\_on\_error} can be given to control whether the program throws an exception for a failed energy calculation or simply returns \texttt{false} for the task and proceeds with the remaining tasks. The default value is \texttt{true}. +Like in any task, the setting \texttt{silent\_stdout\_calculator} can be given to control whether any standard output of the calculation itself should be printed. The default value is \texttt{true}. \subsection{Hessian Calculation} @@ -616,6 +625,7 @@ \subsection{Hessian Calculation} as the backend program, the molecular symmetry number $\sigma$ has to be specified there if it differs from the default value of one. See the \textsc{Sparrow} manual for details. Like in any task, the setting \texttt{stop\_on\_error} can be given to control whether the program throws an exception for a failed energy calculation or simply returns \texttt{false} for the task and proceeds with the remaining tasks. The default value is \texttt{true}. +Like in any task, the setting \texttt{silent\_stdout\_calculator} can be given to control whether any standard output of the calculation itself should be printed. The default value is \texttt{true}. \subsection{Structure Optimization} @@ -648,6 +658,21 @@ \subsection{Structure Optimization} It is given as a list containing the corresponding atom indices (e.g., \texttt{[0, 12, 32, 42]}). This setting can only be set for a true Cartesian coordinate system. By default, this list is empty. \item \texttt{stop\_on\_error}: Determine whether the program throws an exception for failed energy calculations or optimizations or simply returns \texttt{false} for the failed task and proceeds with the remaining tasks. In case of optimizations, the resulting structure after the maximum number of optimization steps has been reached, is still safed, when this is set to \texttt{false}. The default value is \texttt{true}. +\item \texttt{silent\_stdout\_calculator}: Determine whether any standard output of the calculation itself should be printed. The default value is \texttt{true}. +\item \texttt{unitcelloptimizer}: This sets the desired optimization algorithm for the unit cell. You can set \texttt{'bfgs'} +for the BFGS algorithm including G-DIIS, \texttt{'lbfgs'} for the L-BFGS algorithm, and \texttt{'steepestdescent'} or \texttt{'sd'} +for a steepest descent algorithm. By default, it is deactivated. A unit cell optimization requires a program that support the calculation +of a stress tensor. +\end{itemize} +If you specified a \texttt{unitcelloptimizer}, you can also set the following options: +\begin{itemize} +\item \texttt{cellopt\_optimize\_angles"}: Specify whether the angles of the unit cell shall be fixed (true) or constrained (false). By default set to \texttt{true}. +\item \texttt{cellopt\_optimize\_a"}: Specify whether the \textit{a} vector of the unit cell (defined along the positive x-axis) shall be optimized. By default set to \texttt{'true'}. +\item \texttt{cellopt\_optimize\_b"}: Specify whether the \textit{b} vector of the unit cell (defined x-y-plane) shall be optimized. By default set to \texttt{true}. +\item \texttt{cellopt\_optimize\_c"}: Specify whether the \textit{c} vector of the unit cell (defined in the positive z-direction) shall be optimized. By default set to \texttt{true}. +\item \texttt{cellopt\_geoopt\_max\_convergence\_iterations"}: Define the maximum number of microiterations in the optimization of the geometry. By default set to \texttt{100}. +\item \texttt{cellopt\_cellopt\_max\_convergence\_iterations"}: Define the maximum number of microiterations in the optimization of the unitcell. By default set to \texttt{100}. +The total number of iterations is still controlled by \texttt{geoopt\_max\_convergence\_iterations}. \end{itemize} If you specified \texttt{optimizer: 'bfgs'}, the default coordinate system is 'internal' @@ -656,8 +681,12 @@ \subsection{Structure Optimization} \item \texttt{bfgs\_use\_gdiis}: Switch to enable the use of a G-DIIS possibly accelerating convergence. By default set to \texttt{true}. \item \texttt{bfgs\_gdiis\_max\_store}: The maximum number of old steps used in the G-DIIS. By default set to \texttt{5}. -\item \texttt{bfgs\_use\_trust\_radius}: Whether to use the trust radius. By default set to \texttt{false}. -\item \texttt{bfgs\_trust\_radius}: The maximum movement in any cartesian direction by any atom. By default set to \texttt{0.1}. +\item \texttt{bfgs\_use\_trust\_radius}: Whether to use the trust radius. By default set to \texttt{true}. +\item \texttt{bfgs\_trust\_radius}: The maximum movement in any cartesian direction by any atom. By default set to \texttt{0.3}. +\item \texttt{bfgs\_min\_iterations}: The minimum number of cycles before the convergence criteria are checked. +By default set to \texttt{1}.\\ +Setting this option to values larger than \texttt{1} enforces one or multiple updates of the B matrix. +This might help for structures with small absolute values of imaginary frequencies. \end{itemize} If you specified \texttt{optimizer: 'lbfgs'}, the default coordinate system is 'cartesianWithoutRotTrans' @@ -702,7 +731,9 @@ \subsection{Transition State Optimization} \item \texttt{optimizer}: This sets the desired optimization algorithm. You can set \texttt{'bofill'} for Bofill's algorithm\cite{bofill1, bofill2}, or any of \texttt{'eigenvector\_following'}, \texttt{'eigenvectorfollowing'}, \texttt{'ef'}, \texttt{'evf'}, or \texttt{'ev'} for a eigenvector following algorithm, or \texttt{'dimer'} for the Dimer algorithm~\cite{dimer1,dimer2,dimer3}. By default, it is set to \texttt{'bofill'}. -\item \texttt{automatic\_mode\_selection}: Specify the indices of atoms that are essential for the reaction in a list to automatically select the imaginary frequency mode with the highest movement of these atoms. +\item \texttt{automatic\_mode\_selection}: Specify the indices of atoms that are essential for the reaction in a list to automatically + select the imaginary frequency mode to follow. Out of all modes that have high contributions by movement of the selected atoms + (all within 10\,\% of the highest contributing mode) the one with the biggest imaginary frequency is picked. \item \texttt{convergence\_step\_max\_coefficient}: The convergence threshold for the maximum absolute element of the last step taken. By default set to \texttt{1.0e-4}. \item \texttt{convergence\_step\_rms}: The convergence threshold for the root mean square of the last step taken. By default set to @@ -721,6 +752,7 @@ \subsection{Transition State Optimization} It is given as a list containing the corresponding atom indices (e.g., \texttt{[0, 12, 32, 42]}). This setting can only be set for a true Cartesian coordinate system. By default, this list is empty. \item \texttt{stop\_on\_error}: Determine whether the program throws an exception for failed energy calculations or optimizations or simply returns \texttt{false} for the failed task and proceeds with the remaining tasks. In case of optimizations, the resulting structure after the maximum number of optimization steps has been reached, is still safed, when this is set to \texttt{false}. The default value is \texttt{true}. +\item \texttt{silent\_stdout\_calculator}: Determine whether any standard output of the calculation itself should be printed. The default value is \texttt{true}. \end{itemize} If you specified \texttt{optimizer: 'bofill'}, the default coordinate system is 'cartesianWithoutRotTrans' @@ -822,6 +854,7 @@ \subsection{Intrinsic Reaction Coordinate Calculation} The default for the IRC task is 'cartesianWithoutRotTrans'. \item \texttt{irc\_initial\_step\_size}: Maximum displacement of one coordinate of one atom along the given mode. All other coordinates are scaled down accordingly; by default it is set to \texttt{0.3}. \item \texttt{stop\_on\_error}: Determine whether the program throws an exception for failed energy calculations or optimizations or simply returns \texttt{false} for the failed task and proceeds with the remaining tasks. In case of optimizations, the resulting structure after the maximum number of optimization steps has been reached, is still safed, when this is set to \texttt{false}. Especially, when using a SD type optimizer this option can be helpful. The default value is \texttt{true}. +\item \texttt{silent\_stdout\_calculator}: Determine whether any standard output of the calculation itself should be printed. The default value is \texttt{true}. \end{itemize} If you specified \texttt{optimizer: 'bfgs'}, you can also set the following options: @@ -903,6 +936,7 @@ \subsection{Artificial Force Induced Reaction Calculation} \item \texttt{convergence\_requirement}: The number of criteria that have to converge besides the value criterion. This has to be between \texttt{0} and \texttt{4}; by default it is set to \texttt{3}. \item \texttt{stop\_on\_error}: Determine whether the program throws an exception for failed energy calculations or optimizations or simply returns \texttt{false} for the failed task and proceeds with the remaining tasks. In case of optimizations, the resulting structure after the maximum number of optimization steps has been reached, is still safed, when this is set to \texttt{false}. Especially, when using a SD type optimizer this option can be helpful. The default value is \texttt{true}. +\item \texttt{silent\_stdout\_calculator}: Determine whether any standard output of the calculation itself should be printed. The default value is \texttt{true}. \end{itemize} If you specified \texttt{optimizer: 'bfgs'}, you can also set the following options: @@ -970,6 +1004,7 @@ \subsection{B-Spline Interpolation and Optimization}\label{sec: b-spline} \item \texttt{num\_structures}: Sets the number of structures into which a reaction path is discretized for the final output (\textit{i.e.,} when writing it to a XYZ trajectory file). By default set to \texttt{10}. \item \texttt{stop\_on\_error}: Determine whether the program throws an exception for failed energy calculations or optimizations or simply returns \texttt{false} for the failed task and proceeds with the remaining tasks. In case of optimizations, the resulting structure after the maximum number of optimization steps has been reached, is still safed, when this is set to \texttt{false}. Especially, when using a SD type optimizer this option can be helpful. The default value is \texttt{true}. +\item \texttt{silent\_stdout\_calculator}: Determine whether any standard output of the calculation itself should be printed. The default value is \texttt{true}. \end{itemize} If you specified \texttt{optimizer: 'lbfgs'}, you can also set the following options: @@ -991,12 +1026,20 @@ \subsection{B-Spline Interpolation and Optimization}\label{sec: b-spline} \item \texttt{sd\_factor}: The scaling factor to be used in the steepest descent algorithm. By default set to \texttt{0.1}. \end{itemize} -\subsection{Newton Trajectory Calculation} +\subsection{Newton Trajectory Calculations} + +These tasks are used to generate transition state guesses. +To this end, sets of atoms are force onto one another, or pushed apart. +There are currently two algorithms available for these types of calculations: +Newton Trajectory 1 (\texttt{type: 'nt'}, \texttt{type: 'nt1'}) and Newton Trajectory 2 (\texttt{type: 'nt2'}). + +\subsubsection{Newton Trajectory Algorithm 1} This task is used to generate transition state guesses. Two groups of atoms are forced onto one another. In order to carry out this task, specify any of the following in the respective task block: -\texttt{type: 'newtontrajectory'}, \texttt{type: 'ntoptimization'}, \texttt{type: 'ntopt'}, or \texttt{type: 'nt'}. +\texttt{type: 'newtontrajectory'}, \texttt{type: 'ntoptimization'}, \texttt{type: 'ntopt'}, \texttt{type: 'nt'} or \texttt{type: 'nt1'}. +In this version of the algorithm two groups of atoms (LHS and RHS) are pused together or pulled appart. You usually want to set the following settings: \begin{itemize} \item \texttt{nt\_rhs\_list}: This specifies list of indices of atoms to be forced onto or away from those @@ -1048,6 +1091,57 @@ \subsection{Newton Trajectory Calculation} covalent radii sums, or if the geometric centers of the two groups are further than the value in a.u. apart. By default the value is set to \texttt{4.0}. \item \texttt{stop\_on\_error}: Determine whether the program throws an exception for failed energy calculations or optimizations or simply returns \texttt{false} for the failed task and proceeds with the remaining tasks. In case of optimizations, the resulting structure after the maximum number of optimization steps has been reached, is still safed, when this is set to \texttt{false}. Especially, when using a SD type optimizer this option can be helpful. The default value is \texttt{true}. +\item \texttt{silent\_stdout\_calculator}: Determine whether any standard output of the calculation itself should be printed. The default value is \texttt{true}. +\end{itemize} + +The Newton trajectory task uses a steepest descent algorithm internally; for this reason the following setting is available: +\begin{itemize} +\item \texttt{sd\_factor}: The scaling factor to be used in the steepest descent algorithm. By default set to \texttt{1.0}. +\end{itemize} + +\subsubsection{Newton Trajectory Algorithm 2} +In order to carry out this task, specify any of the following in the respective task block: +\texttt{type: 'newtontrajectory2'}, \texttt{type: 'ntoptimization2'}, \texttt{type: 'ntopt2'}, or \texttt{type: 'nt2'}. + +In this version of the task any number of pairs of atoms are pushed together or pulled apart. +You usually want to set the following settings: +\begin{itemize} +\item \texttt{nt\_associations}: This specifies list of indices of atoms pairs to be forced onto another. + Pairs are given in a flat list such that indices $2i$ and $2i+1$ form a pair. +\item \texttt{nt\_dissociations}: This specifies list of indices of atoms pairs to be forced away from another. + Pairs are given in a flat list such that indices $2i$ and $2i+1$ form a pair. +\end{itemize} + +The task works without the specification of any additional settings; the default settings work usually fine. However, +if desired, the following settings can always be set: +\begin{itemize} +\item \texttt{nt\_total\_force\_norm}: The maximum norm of the force applied to an atom listed in any of the pairs. + Note that some atoms may have forces applied that are less than the given value, all forces applied will be scaled + such that the 'fastest' atom will be receiving at maximum a force with the norm given in this setting. + By default it is set to \texttt{0.1}. +\item \texttt{nt\_coordinate\_system}: Select the representation of the coordinate system, either + 'internal'\cite{libirc}, 'cartesianWithoutRotTrans', or 'cartesian'. +Please be aware that this only affects the SD steps, while the BFGS micro cycles will always be performed with true + Cartesian coordinates due to constraints. +\item \texttt{nt\_constrained\_atoms}: A list of atoms for which the Cartesian coordinates are constrained during the optimization. + It is given as a list containing the corresponding atom indices (e.g., \texttt{[0, 12, 32, 42]}). This setting can only be set for a true Cartesian coordinate system. +\item \texttt{nt\_use\_micro\_cycles}: Use micro cycles inbetween forced steps that move the constrained atoms. +In these micro cycles a BFGS/GDIIS will be used to optimize the geometry with the lhs/rhs atoms fixed in place. +By default it is set to \texttt{true}. +\item \texttt{nt\_fixed\_number\_of\_micro\_cycles}: Use a fixed number of micro cycles per macro cycle. If set to +\texttt{false}, the number of micro cycles will grow with the number of macro iterations. It will grow by one micro +cycle per macro iteration performed until \texttt{nt\_number\_of\_micro\_cycles} is reached. By default it is set to \texttt{true}. +\item \texttt{nt\_number\_of\_micro\_cycles}: The maximum number of micro cycles used. By default it is set to \texttt{10}. +\item \texttt{nt\_filter\_passes}: This task uses a Savitzky--Golay filter before analyzing the reaction curve. The number of passes +through this filter is adjusted by the \texttt{nt\_filter\_passes} option. By default the number of passes is \texttt{10}. +\item \texttt{convergence\_max\_iterations}: The maximum number of Newton trajectory macro iterations. By default set to \texttt{500}. +\item \texttt{convergence\_attractive\_stop}: The distance between the forced atoms at which the algorithm stops and tries to identify +a transition state guess. The given value is applied as a factor to the covalent radii of the compared atoms. +In the case of an attractive force between two atoms, this criterion set to 1.0 would stop the NT procedure as soon as the two atoms +are closer than 1.0 times the sum of their covalent radii. +By default the value it is set to \texttt{0.9}. +\item \texttt{stop\_on\_error}: Determine whether the program throws an exception for failed energy calculations or optimizations or simply returns \texttt{false} for the failed task and proceeds with the remaining tasks. In case of optimizations, the resulting structure after the maximum number of optimization steps has been reached, is still safed, when this is set to \texttt{false}. Especially, when using a SD type optimizer this option can be helpful. The default value is \texttt{true}. +\item \texttt{silent\_stdout\_calculator}: Determine whether any standard output of the calculation itself should be printed. The default value is \texttt{true}. \end{itemize} The newton trajectory task uses a steepest descent algorithm internally; for this reason the following setting is available: diff --git a/src/Readuct/App/Tasks/AfirOptimizationTask.h b/src/Readuct/App/Tasks/AfirOptimizationTask.h index f7d49d9..5d93aad 100644 --- a/src/Readuct/App/Tasks/AfirOptimizationTask.h +++ b/src/Readuct/App/Tasks/AfirOptimizationTask.h @@ -48,10 +48,12 @@ class AfirOptimizationTask : public Task { warningIfMultipleInputsGiven(); warningIfMultipleOutputsGiven(); + bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true); std::shared_ptr calc; if (!testRunOnly) { // leave out in case of task chaining --> attention calc is NULL // Note: _input is guaranteed not to be empty by Task constructor calc = copyCalculator(systems, _input.front(), name()); + Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator); // Check system size if (calc->getStructure()->size() == 1) { @@ -66,11 +68,13 @@ class AfirOptimizationTask : public Task { if (optimizertype == "LBFGS") { auto tmp = std::make_shared>(*calc); tmp->optimizer.useTrustRadius = true; + tmp->optimizer.trustRadius = 0.1; optimizer = std::move(tmp); } else if (optimizertype == "BFGS") { auto tmp = std::make_shared>(*calc); tmp->optimizer.useTrustRadius = true; + tmp->optimizer.trustRadius = 0.1; optimizer = std::move(tmp); } else if (optimizertype == "SD" || optimizertype == "STEEPESTDESCENT") { diff --git a/src/Readuct/App/Tasks/BSplineInterpolationTask.h b/src/Readuct/App/Tasks/BSplineInterpolationTask.h index 31f8e97..9e1e688 100644 --- a/src/Readuct/App/Tasks/BSplineInterpolationTask.h +++ b/src/Readuct/App/Tasks/BSplineInterpolationTask.h @@ -63,6 +63,7 @@ class BSplineInterpolationTask : public Task { tangentFileName_ = taskSettings.extract("tangent_file", tangentFileName_); coordinateThresholdForMaximumExtraction_ = taskSettings.extract("extract_threshold", coordinateThresholdForMaximumExtraction_); + bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true); // If no errors encountered until here, the basic settings should be alright if (testRunOnly) { @@ -72,6 +73,8 @@ class BSplineInterpolationTask : public Task { // Note: _input is guaranteed not to be empty by Task constructor auto calc = copyCalculator(systems, _input.front(), name()); auto secondCalculator = systems.at(_input.back()); + Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator); + Utils::CalculationRoutines::setLog(*secondCalculator, true, true, !silentCalculator); if (calc->settings() != secondCalculator->settings()) { _logger->warning << " Warning: The given systems have different settings. Only taking first and ignoring second.\n"; @@ -296,7 +299,7 @@ class BSplineInterpolationTask : public Task { "this optimization Task"); } - const int cycles = optimizer->optimize(); + const int cycles = optimizer->optimize(*_logger); const int maxiter = settings.getInt("convergence_max_iterations"); bool converged = cycles < maxiter; diff --git a/src/Readuct/App/Tasks/BondOrderTask.h b/src/Readuct/App/Tasks/BondOrderTask.h index 0d2606a..91713e5 100644 --- a/src/Readuct/App/Tasks/BondOrderTask.h +++ b/src/Readuct/App/Tasks/BondOrderTask.h @@ -11,7 +11,7 @@ #include "Tasks/Task.h" /* Scine */ #include -#include +#include /* External */ #include "boost/exception/diagnostic_information.hpp" #include @@ -27,6 +27,9 @@ class BondOrderTask : public Task { public: /** * @brief Construct a new BondOrderTask. + * + * @deprecated Please use the SinglePointTask with the `require_bond_orders` setting. + * * @param input The input system names for the task. * @param output The output system names for the task. * @param logger The logger to/through which all text output will be handled. @@ -45,6 +48,7 @@ class BondOrderTask : public Task { // Read and delete special settings bool stopOnError = stopOnErrorExtraction(taskSettings); + bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true); if (!taskSettings.empty()) { throw std::logic_error(falseTaskSettingsErrorMessage(name())); } @@ -55,6 +59,7 @@ class BondOrderTask : public Task { // Note: _input is guaranteed not to be empty by Task constructor auto calc = copyCalculator(systems, _input.front(), name()); + Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator); // Calculate bond orders and energy if not present in the results yet if (!calc->results().has()) { diff --git a/src/Readuct/App/Tasks/GeometryOptimizationTask.h b/src/Readuct/App/Tasks/GeometryOptimizationTask.h index 6b2e926..21f268b 100644 --- a/src/Readuct/App/Tasks/GeometryOptimizationTask.h +++ b/src/Readuct/App/Tasks/GeometryOptimizationTask.h @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -48,41 +49,19 @@ class GeometryOptimizationTask : public Task { bool run(SystemsMap& systems, Utils::UniversalSettings::ValueCollection taskSettings, bool testRunOnly = false) const final { warningIfMultipleInputsGiven(); warningIfMultipleOutputsGiven(); + bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true); // Get/Copy Calculator std::shared_ptr calc; if (!testRunOnly) { // leave out in case of task chaining --> attention calc is NULL // Note: _input is guaranteed not to be empty by Task constructor calc = copyCalculator(systems, _input.front(), name()); + Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator); } // Generate optimizer auto optimizertype = taskSettings.extract("optimizer", std::string{"BFGS"}); - std::transform(optimizertype.begin(), optimizertype.end(), optimizertype.begin(), ::toupper); - std::shared_ptr optimizer; - if (optimizertype == "LBFGS") { - auto tmp = std::make_shared>(*calc); - // Default convergence options - optimizer = std::move(tmp); - } - else if (optimizertype == "BFGS") { - auto tmp = std::make_shared>(*calc); - // Default convergence options - optimizer = std::move(tmp); - } - else if (optimizertype == "SD" || optimizertype == "STEEPESTDESCENT") { - auto tmp = std::make_shared>(*calc); - // Default convergence options - optimizer = std::move(tmp); - } - else if (optimizertype == "NR" || optimizertype == "NEWTONRAPHSON") { - auto tmp = std::make_shared>(*calc); - // Default convergence options - optimizer = std::move(tmp); - } - else { - throw std::runtime_error( - "Unknown Optimizer requested for a geometry optimization, available are: SD, NR, BFGS and LBFGS!"); - } + auto unitcelloptimizertype = taskSettings.extract("unitcelloptimizer", std::string{""}); + auto optimizer = constructOptimizer(*calc, optimizertype, unitcelloptimizertype); // Read and delete special settings bool stopOnError = stopOnErrorExtraction(taskSettings); @@ -186,6 +165,86 @@ class GeometryOptimizationTask : public Task { return cycles < maxiter; } + + inline std::shared_ptr constructOptimizer(Core::Calculator& calc, std::string type, + std::string cellType) const { + std::transform(type.begin(), type.end(), type.begin(), ::toupper); + std::transform(cellType.begin(), cellType.end(), cellType.begin(), ::toupper); + if (type == "LBFGS") { + if (cellType.empty()) { + return std::make_shared>(calc); + } + else { + if (cellType == "LBFGS") { + return std::make_shared>(calc); + } + else if (cellType == "BFGS") { + return std::make_shared>(calc); + } + else if (cellType == "SD" || cellType == "STEEPESTDESCENT") { + return std::make_shared>(calc); + } + throw std::runtime_error( + "Unknown CellOptimizer requested for a geometry optimization, available are: SD, BFGS and LBFGS!"); + } + } + else if (type == "BFGS") { + if (cellType.empty()) { + return std::make_shared>(calc); + } + else { + if (cellType == "LBFGS") { + return std::make_shared>(calc); + } + else if (cellType == "BFGS") { + return std::make_shared>(calc); + } + else if (cellType == "SD" || cellType == "STEEPESTDESCENT") { + return std::make_shared>(calc); + } + throw std::runtime_error( + "Unknown CellOptimizer requested for a geometry optimization, available are: SD, BFGS and LBFGS!"); + } + } + else if (type == "SD" || type == "STEEPESTDESCENT") { + if (cellType.empty()) { + return std::make_shared>(calc); + } + else { + if (cellType == "LBFGS") { + return std::make_shared>(calc); + } + else if (cellType == "BFGS") { + return std::make_shared>(calc); + } + else if (cellType == "SD" || cellType == "STEEPESTDESCENT") { + return std::make_shared>(calc); + } + throw std::runtime_error( + "Unknown CellOptimizer requested for a geometry optimization, available are: SD, BFGS and LBFGS!"); + } + } + else if (type == "NR" || type == "NEWTONRAPHSON") { + if (cellType.empty()) { + return std::make_shared>(calc); + } + else { + if (cellType == "LBFGS") { + return std::make_shared>(calc); + } + else if (cellType == "BFGS") { + return std::make_shared>(calc); + } + else if (cellType == "SD" || cellType == "STEEPESTDESCENT") { + return std::make_shared>(calc); + } + throw std::runtime_error( + "Unknown CellOptimizer requested for a geometry optimization, available are: SD, BFGS and LBFGS!"); + } + } + throw std::runtime_error( + "Unknown Optimizer requested for a geometry optimization, available are: SD, NR, BFGS and LBFGS!"); + }; }; } // namespace Readuct diff --git a/src/Readuct/App/Tasks/HessianTask.h b/src/Readuct/App/Tasks/HessianTask.h index ecbc23a..cbb5fe3 100644 --- a/src/Readuct/App/Tasks/HessianTask.h +++ b/src/Readuct/App/Tasks/HessianTask.h @@ -50,6 +50,7 @@ class HessianTask : public Task { // Read and delete special settings bool stopOnError = stopOnErrorExtraction(taskSettings); + bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true); if (!taskSettings.empty()) { throw std::logic_error(falseTaskSettingsErrorMessage(name())); } @@ -61,6 +62,7 @@ class HessianTask : public Task { // Note: _input is guaranteed not to be empty by Task constructor auto calc = copyCalculator(systems, _input.front(), name()); + Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator); // Check system size if (calc->getStructure()->size() == 1) { diff --git a/src/Readuct/App/Tasks/IrcTask.h b/src/Readuct/App/Tasks/IrcTask.h index 72ce0ea..609e2d5 100644 --- a/src/Readuct/App/Tasks/IrcTask.h +++ b/src/Readuct/App/Tasks/IrcTask.h @@ -53,10 +53,12 @@ class IrcTask : public Task { << " Warning: To store IRC results in a new location two output systems have to be specified.\n"; } + bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true); std::shared_ptr calc; if (!testRunOnly) { // leave out in case of task chaining --> attention calc is NULL // Note: _input is guaranteed not to be empty by Task constructor calc = copyCalculator(systems, _input.front(), name()); + Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator); // Check system size if (calc->getStructure()->size() == 1) { @@ -216,6 +218,9 @@ class IrcTask : public Task { // Add observer oldEnergy = 0.0; oldParams.resize(0); + // Reset optimizer + optimizer->setSettings(settings); + optimizer->reset(); // Trajectory stream boost::filesystem::path dirB(((_output.size() > 1) ? _output[1] : _input[0])); boost::filesystem::create_directory(dirB); @@ -241,6 +246,7 @@ class IrcTask : public Task { // Run optimization structure = systems.at(_input[0])->getStructure(); + cycles = 0; try { cycles = optimizer->optimize(*structure, *_logger, ircModeVector, false); } diff --git a/src/Readuct/App/Tasks/NtOptimization2Task.h b/src/Readuct/App/Tasks/NtOptimization2Task.h new file mode 100644 index 0000000..89667ad --- /dev/null +++ b/src/Readuct/App/Tasks/NtOptimization2Task.h @@ -0,0 +1,147 @@ +/** + * @file + * @copyright This code is licensed under the 3-clause BSD license.\n + * Copyright ETH Zurich, Laboratory of Physical Chemistry, Reiher Group.\n + * See LICENSE.txt for details. + */ +#ifndef READUCT_NTOPTIMIZATION2TASK_H_ +#define READUCT_NTOPTIMIZATION2TASK_H_ + +/* Readuct */ +#include "Tasks/Task.h" +/* Scine */ +#include +#include +#include +/* External */ +#include +#include +/* std c++ */ +#include +#include +#include + +namespace Scine { +namespace Readuct { + +class NtOptimization2Task : public Task { + public: + /** + * @brief Construct a new NtOptimization2Task. + * @param input The input system names for the task. + * @param output The output system names for the task. + * @param logger The logger to/through which all text output will be handled. + */ + NtOptimization2Task(std::vector input, std::vector output, std::shared_ptr logger = nullptr) + : Task(std::move(input), std::move(output), std::move(logger)) { + } + + std::string name() const override { + return "NT2 Optimization"; + } + + bool run(SystemsMap& systems, Utils::UniversalSettings::ValueCollection taskSettings, bool testRunOnly = false) const final { + warningIfMultipleInputsGiven(); + warningIfMultipleOutputsGiven(); + + bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true); + std::shared_ptr calc; + if (!testRunOnly) { // leave out in case of task chaining --> attention calc is NULL + // Note: _input is guaranteed not to be empty by Task constructor + calc = copyCalculator(systems, _input.front(), name()); + Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator); + + // Check system size + if (calc->getStructure()->size() == 1) { + throw std::runtime_error("Cannot calculate NT2 optimization for monoatomic systems."); + } + } + + // Generate optimizer + auto optimizer = std::make_unique(*calc); + // Read and delete special settings + bool stopOnError = stopOnErrorExtraction(taskSettings); + // Apply user settings + auto settings = optimizer->getSettings(); + settings.merge(taskSettings); + if (!settings.valid()) { + settings.throwIncorrectSettings(); + } + optimizer->setSettings(settings); + + // If no errors encountered until here, the basic settings should be alright + if (testRunOnly) { + return true; + } + + // Add observer + // Trajectory stream + using Writer = Utils::XyzStreamHandler; + const std::string& outputSystem = ((!_output.empty()) ? _output[0] : _input[0]); + boost::filesystem::path dir(outputSystem); + boost::filesystem::create_directory(dir); + boost::filesystem::path trjfile(outputSystem + ".nt.trj.xyz"); + std::ofstream trajectory((dir / trjfile).string(), std::ofstream::out); + double oldEnergy = 0.0; + Eigen::VectorXd oldParams; + auto func = [&](const int& cycle, const double& energy, const Eigen::VectorXd& params) { + if (oldParams.size() != params.size()) { + oldParams = params; + } + if (cycle == 1) { + printf("%7s %16s %16s %16s %16s\n", "Cycle", "Energy", "Energy Diff.", "Step RMS", "Max. Step"); + } + auto diff = (params - oldParams).eval(); + printf("%7d %+16.9f %+16.9f %+16.9f %+16.9f\n", cycle, energy, energy - oldEnergy, + sqrt(diff.squaredNorm() / diff.size()), diff.cwiseAbs().maxCoeff()); + oldEnergy = energy; + oldParams = params; + auto structure = calc->getStructure(); + Writer::write(trajectory, *structure, std::to_string(energy)); + }; + optimizer->addObserver(func); + auto cout = _logger->output; + + // Run optimization + auto structure = calc->getStructure(); + int cycles = 0; + try { + cycles = optimizer->optimize(*structure, *_logger); + } + catch (...) { + trajectory.close(); + if (stopOnError) { + throw; + } + // check if thrown exception corresponds to no actual calculation failures but simply no guess found + // --> no thrown exception intended in this case if stopOnError equals false + std::string noGuess = "No transition state guess was found in Newton Trajectory scan"; + size_t found = boost::current_exception_diagnostic_information().find(noGuess); + if (found == std::string::npos) { // did not find harmless exception -> throw + throw; + } + cout << Core::Log::endl << " No TS guess found in NT2 scan." << Core::Log::endl << Core::Log::endl; + return false; + } + + trajectory.close(); + + cout << Core::Log::endl + << " Found TS guess after " << cycles << " iterations." << Core::Log::endl + << Core::Log::endl; + + // Print/Store results + systems[outputSystem] = calc; + boost::filesystem::path xyzfile(outputSystem + ".xyz"); + std::ofstream xyz((dir / xyzfile).string(), std::ofstream::out); + Writer::write(xyz, *(calc->getStructure())); + xyz.close(); + + return true; + } +}; + +} // namespace Readuct +} // namespace Scine + +#endif // READUCT_NTOPTIMIZATION2TASK_H_ diff --git a/src/Readuct/App/Tasks/NtOptimizationTask.h b/src/Readuct/App/Tasks/NtOptimizationTask.h index 6c05c3c..9efea5c 100644 --- a/src/Readuct/App/Tasks/NtOptimizationTask.h +++ b/src/Readuct/App/Tasks/NtOptimizationTask.h @@ -37,21 +37,23 @@ class NtOptimizationTask : public Task { } std::string name() const override { - return "NT Optimization"; + return "NT1 Optimization"; } bool run(SystemsMap& systems, Utils::UniversalSettings::ValueCollection taskSettings, bool testRunOnly = false) const final { warningIfMultipleInputsGiven(); warningIfMultipleOutputsGiven(); + bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true); std::shared_ptr calc; if (!testRunOnly) { // leave out in case of task chaining --> attention calc is NULL // Note: _input is guaranteed not to be empty by Task constructor calc = copyCalculator(systems, _input.front(), name()); + Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator); // Check system size if (calc->getStructure()->size() == 1) { - throw std::runtime_error("Cannot calculate NT optimization for monoatomic systems."); + throw std::runtime_error("Cannot calculate NT1 optimization for monoatomic systems."); } } @@ -118,26 +120,15 @@ class NtOptimizationTask : public Task { if (found == std::string::npos) { // did not find harmless exception -> throw throw; } - cout << Core::Log::endl << " No TS guess found in NT scan." << Core::Log::endl << Core::Log::endl; + cout << Core::Log::endl << " No TS guess found in NT1 scan." << Core::Log::endl << Core::Log::endl; return false; } trajectory.close(); - int maxiter = settings.getInt("convergence_max_iterations"); - if (cycles < maxiter) { - cout << Core::Log::endl - << " Converged after " << cycles << " iterations." << Core::Log::endl - << Core::Log::endl; - } - else { - cout << Core::Log::endl - << " Stopped after " << maxiter << " iterations." << Core::Log::endl - << Core::Log::endl; - if (stopOnError) { - throw std::runtime_error("Problem: NT optimization did not converge."); - } - } + cout << Core::Log::endl + << " Found TS guess after " << cycles << " iterations." << Core::Log::endl + << Core::Log::endl; // Print/Store results systems[outputSystem] = calc; @@ -146,7 +137,7 @@ class NtOptimizationTask : public Task { Writer::write(xyz, *(calc->getStructure())); xyz.close(); - return cycles < maxiter; + return true; } }; diff --git a/src/Readuct/App/Tasks/SinglePointTask.h b/src/Readuct/App/Tasks/SinglePointTask.h index 88b714e..de67b6f 100644 --- a/src/Readuct/App/Tasks/SinglePointTask.h +++ b/src/Readuct/App/Tasks/SinglePointTask.h @@ -11,7 +11,7 @@ #include "Tasks/Task.h" /* Scine */ #include -#include +#include #include #include #include @@ -56,9 +56,18 @@ class SinglePointTask : public Task { bool stopOnError = stopOnErrorExtraction(taskSettings); bool requireCharges = taskSettings.extract("require_charges", true); bool requireGradients = taskSettings.extract("require_gradients", false); + bool requireStressTensor = taskSettings.extract("require_stress_tensor", false); + bool requireBondOrders = taskSettings.extract("require_bond_orders", false); bool requireOrbitalEnergies = taskSettings.extract("orbital_energies", false); + bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true); + int spinPropensityCheck = taskSettings.extract("spin_propensity_check", 0); if (!taskSettings.empty()) { - throw std::logic_error("Specified a task setting that is not available for this task"); + std::string keyListing = "\n"; + auto keys = taskSettings.getKeys(); + for (const auto& key : keys) { + keyListing += "'" + key + "'\n"; + } + throw std::logic_error("Specified one or more task settings that are not available for this task:" + keyListing); } if (testRunOnly) { @@ -67,11 +76,14 @@ class SinglePointTask : public Task { // Note: _input is guaranteed not to be empty by Task constructor auto calc = copyCalculator(systems, _input.front(), name()); + Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator); // Check for available properties bool chargesAvailable = calc->possibleProperties().containsSubSet(Utils::Property::AtomicCharges); bool gradientsAvailable = calc->possibleProperties().containsSubSet(Utils::Property::Gradients); + bool stressTensorAvailable = calc->possibleProperties().containsSubSet(Utils::Property::StressTensor); bool orbitalEnergiesAvailable = calc->possibleProperties().containsSubSet(Utils::Property::OrbitalEnergies); + bool bondOrdersAvailable = calc->possibleProperties().containsSubSet(Utils::Property::BondOrderMatrix); if (requireCharges && !chargesAvailable) { throw std::logic_error("Charges required, but chosen calculator does not provide them.\n" "If you do not need charges, set 'require_charges' to 'false' in the task settings"); @@ -79,22 +91,32 @@ class SinglePointTask : public Task { if (requireGradients && !gradientsAvailable) { throw std::logic_error("Gradients required, but chosen calculator does not provide them."); } + if (requireStressTensor && !stressTensorAvailable) { + throw std::logic_error("Stress tensor required, but chosen calculator does not provide it."); + } if (requireOrbitalEnergies && !orbitalEnergiesAvailable) { throw std::logic_error("Orbital energies required, but chosen calculator does not provide them."); } - + if (requireBondOrders && !bondOrdersAvailable) { + throw std::logic_error("Bond orders required, but chosen calculator does not provide them."); + } // Calculate energy, and possibly atomic charges and/or gradients Utils::PropertyList requiredProperties = Utils::Property::Energy; if (requireCharges) { requiredProperties.addProperty(Utils::Property::AtomicCharges); } - if (requireGradients) { requiredProperties.addProperty(Utils::Property::Gradients); } + if (requireStressTensor) { + requiredProperties.addProperty(Utils::Property::StressTensor); + } if (requireOrbitalEnergies) { requiredProperties.addProperty(Utils::Property::OrbitalEnergies); } + if (requireBondOrders) { + requiredProperties.addProperty(Utils::Property::BondOrderMatrix); + } try { calc->setRequiredProperties(requiredProperties); @@ -107,12 +129,35 @@ class SinglePointTask : public Task { if (stopOnError) { throw; } - _logger->warning + _logger->error << " " + name() + " was not successful with error:\n " + boost::current_exception_diagnostic_information() << Core::Log::endl; return false; } + if (spinPropensityCheck != 0) { + calc = Utils::CalculationRoutines::spinPropensity(*calc, *_logger, spinPropensityCheck); + // some calculators are lazy and don't copy properties, hence we recalculate if necessary + if (!calc->getRequiredProperties().containsSubSet(requiredProperties)) { + try { + calc->setRequiredProperties(requiredProperties); + calc->calculate(name()); + if (!calc->results().get()) { + throw std::runtime_error(name() + " was not successful"); + } + } + catch (...) { + if (stopOnError) { + throw; + } + _logger->error + << " " + name() + " was not successful with error:\n " + boost::current_exception_diagnostic_information() + << Core::Log::endl; + return false; + } + } + } + // Store result if (!_output.empty()) { systems[_output[0]] = calc; @@ -147,6 +192,28 @@ class SinglePointTask : public Task { auto gradients = calc->results().get(); cout << " Gradients (hartree / bohr):\n\n"; cout << [&gradients](std::ostream& os) { Utils::matrixPrettyPrint(os, gradients); }; + cout << Core::Log::nl; + } + // Print stress tensor + if (requireStressTensor) { + Eigen::Matrix3d stressTensor = calc->results().get(); + cout << " Stress tensor (hartree / bohr^3):\n\n"; + cout << [&stressTensor](std::ostream& os) { Utils::matrixPrettyPrint(os, stressTensor); }; + cout << Core::Log::nl; + } + // Print bond orders + if (requireBondOrders) { + cout << " Atom#1 Atom#2 Bond Order\n\n"; + auto bos = calc->results().get(); + const auto& mat = bos.getMatrix(); + for (int i = 0; i < mat.outerSize(); i++) { + for (typename Eigen::SparseMatrix::InnerIterator it(mat, i); it; ++it) { + if (it.value() > 0.3 && it.row() < it.col()) { + cout.printf(" %6d %6d %+16.9f\n", it.row(), it.col(), it.value()); + } + } + } + cout << Core::Log::nl << Core::Log::endl; } // Print orbital energies if (requireOrbitalEnergies) { diff --git a/src/Readuct/App/Tasks/Task.h b/src/Readuct/App/Tasks/Task.h index 0c51391..4a2dc00 100644 --- a/src/Readuct/App/Tasks/Task.h +++ b/src/Readuct/App/Tasks/Task.h @@ -115,10 +115,13 @@ class Task { }; static std::string falseTaskSettingsErrorMessage(const std::string& name) { - return " You gave Task settings for the " + name + ",\n" + " but the only possible setting for this task, is the\n" + - " 'stop_on_error' option to control whether Readuct fails\n" + - " with a failed calculation or simply returns false.\n" + - " You might want to specify the settings you put into the task settings\n" + " in the systems section."; + return " You gave Task settings for the " + name + ",\n" + " but the only possible setting for this task, are the\n" + + " 'stop_on_error' option to control whether ReaDuct fails\n" + + " with a failed calculation or simply returns false\n" + + " and the 'silent_stdout_calculator' option to control whether\n" + " the standard output of the calculator should be printed.\n" + " You might want to specify the settings you put into the task settings\n" + + " in the systems section."; } protected: diff --git a/src/Readuct/App/Tasks/TaskFactory.h b/src/Readuct/App/Tasks/TaskFactory.h index 2ff8980..5c31b57 100644 --- a/src/Readuct/App/Tasks/TaskFactory.h +++ b/src/Readuct/App/Tasks/TaskFactory.h @@ -13,6 +13,7 @@ #include "Tasks/GeometryOptimizationTask.h" #include "Tasks/HessianTask.h" #include "Tasks/IrcTask.h" +#include "Tasks/NtOptimization2Task.h" #include "Tasks/NtOptimizationTask.h" #include "Tasks/SinglePointTask.h" #include "Tasks/Task.h" @@ -62,9 +63,13 @@ class TaskFactory { task = std::make_unique(input, output, logger); } else if (name == "NT" || name == "NEWTONTRAJECTORY" || name == "NTOPTIMIZATION" || - name == "NEWTONTRAJECTORYOPTIMIZATION" || name == "NTOPT") { + name == "NEWTONTRAJECTORYOPTIMIZATION" || name == "NTOPT" || name == "NT1") { task = std::make_unique(input, output); } + else if (name == "NT2" || name == "NEWTONTRAJECTORY2" || name == "NTOPTIMIZATION2" || + name == "NEWTONTRAJECTORYOPTIMIZATION2" || name == "NTOPT2") { + task = std::make_unique(input, output); + } else if (name == "HESSIAN" || name == "FREQUENCY_ANALYSIS" || name == "FREQUENCYANALYSIS" || name == "FREQ" || name == "FREQUENCY" || name == "FREQUENCIES") { task = std::make_unique(input, output, logger); diff --git a/src/Readuct/App/Tasks/TsOptimizationTask.h b/src/Readuct/App/Tasks/TsOptimizationTask.h index 2a205da..1e26108 100644 --- a/src/Readuct/App/Tasks/TsOptimizationTask.h +++ b/src/Readuct/App/Tasks/TsOptimizationTask.h @@ -52,10 +52,12 @@ class TsOptimizationTask : public Task { warningIfMultipleInputsGiven(); warningIfMultipleOutputsGiven(); + bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true); std::shared_ptr calc; if (!testRunOnly) { // leave out in case of task chaining --> attention calc is NULL // Note: _input is guaranteed not to be empty by Task constructor calc = copyCalculator(systems, _input.front(), name()); + Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator); // Check system size if (calc->getStructure()->size() == 1) { @@ -447,6 +449,15 @@ class TsOptimizationTask : public Task { contributions.push_back(contribution); } int selection = std::distance(contributions.begin(), std::max_element(contributions.begin(), contributions.end())); + double tolerance = contributions[selection] * 0.1; + // Contributions are sorted by wave number + // -> pick the first one that fits into the tolerance. + for (unsigned int i = 0; i < contributions.size(); i++) { + if ((contributions[selection] - contributions[i]) < tolerance) { + selection = i; + break; + } + } _logger->output << " Automatically selected mode " << std::to_string(selection) << " with wave number " << std::to_string(waveNumbers[selection]) << " cm^-1" << Core::Log::endl; return selection; diff --git a/src/Readuct/App/main.cpp b/src/Readuct/App/main.cpp index 5195cf0..892b94a 100644 --- a/src/Readuct/App/main.cpp +++ b/src/Readuct/App/main.cpp @@ -59,7 +59,7 @@ int main(int argc, char* argv[]) { // Handle help if (optionsMap.count("help") > 0 || optionsMap.count("config") == 0) { cout << optionsDescription << Core::Log::endl; - return 1; + return 0; } // Load input file @@ -233,14 +233,14 @@ int main(int argc, char* argv[]) { const auto& outpt = tasks[i]->output(); cout << " Input System(s): " + tasks[i]->input()[0] << Core::Log::endl; - for (unsigned int i = 1; i < inpt.size(); i++) { - cout << " " + inpt[i] << Core::Log::endl; + for (unsigned int j = 1; j < inpt.size(); j++) { + cout << " " + inpt[j] << Core::Log::endl; } cout << Core::Log::endl; if (!outpt.empty()) { cout << " Output System(s): " + outpt[0] << Core::Log::endl; - for (unsigned int i = 1; i < outpt.size(); i++) { - cout << " " + outpt[i] << Core::Log::endl; + for (unsigned int j = 1; j < outpt.size(); j++) { + cout << " " + outpt[j] << Core::Log::endl; } cout << Core::Log::endl; } diff --git a/src/Readuct/CMakeLists.txt b/src/Readuct/CMakeLists.txt index 69c1ba7..e0b41fc 100644 --- a/src/Readuct/CMakeLists.txt +++ b/src/Readuct/CMakeLists.txt @@ -17,7 +17,7 @@ function(target_set_rpath) set(options "") set(oneValueArgs "") set(multiValueArgs "TARGETS;UNIX_RPATH") - cmake_parse_arguments(SCINE_RPATH "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) + cmake_parse_arguments(SCINE_RPATH_TARGETS "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) if(APPLE) string(REPLACE $ORIGIN @loader_path osx_rpath "${SCINE_RPATH_UNIX_RPATH}") set_target_properties(${SCINE_RPATH_TARGETS} PROPERTIES diff --git a/src/Readuct/Files.cmake b/src/Readuct/Files.cmake index a8075ee..a99bc6c 100644 --- a/src/Readuct/Files.cmake +++ b/src/Readuct/Files.cmake @@ -6,6 +6,8 @@ set(READUCT_APP_FILES ${CMAKE_CURRENT_SOURCE_DIR}/App/Tasks/BondOrderTask.h ${CMAKE_CURRENT_SOURCE_DIR}/App/Tasks/AfirOptimizationTask.h ${CMAKE_CURRENT_SOURCE_DIR}/App/Tasks/TsOptimizationTask.h + ${CMAKE_CURRENT_SOURCE_DIR}/App/Tasks/NtOptimizationTask.h + ${CMAKE_CURRENT_SOURCE_DIR}/App/Tasks/NtOptimization2Task.h ${CMAKE_CURRENT_SOURCE_DIR}/App/Tasks/HessianTask.h ${CMAKE_CURRENT_SOURCE_DIR}/App/Tasks/IrcTask.h ${CMAKE_CURRENT_SOURCE_DIR}/App/Tasks/SinglePointTask.h diff --git a/src/Readuct/Python/TasksPython.cpp b/src/Readuct/Python/TasksPython.cpp index 4227d51..3a04c2f 100644 --- a/src/Readuct/Python/TasksPython.cpp +++ b/src/Readuct/Python/TasksPython.cpp @@ -10,6 +10,7 @@ #include "Tasks/GeometryOptimizationTask.h" #include "Tasks/HessianTask.h" #include "Tasks/IrcTask.h" +#include "Tasks/NtOptimization2Task.h" #include "Tasks/NtOptimizationTask.h" #include "Tasks/SinglePointTask.h" #include "Tasks/Task.h" @@ -128,6 +129,9 @@ void init_tasks(pybind11::module& m) { m.def("run_nt_task", &run, pybind11::arg("systems"), pybind11::arg("names_of_used_systems"), pybind11::arg("test_run_only") = false); + m.def("run_nt2_task", &run, pybind11::arg("systems"), + pybind11::arg("names_of_used_systems"), pybind11::arg("test_run_only") = false); + m.def("run_afir_task", &run, pybind11::arg("systems"), pybind11::arg("names_of_used_systems"), pybind11::arg("test_run_only") = false); diff --git a/src/Readuct/Python/sphinx/index.rst b/src/Readuct/Python/sphinx/index.rst index c5af5e1..b7c814d 100644 --- a/src/Readuct/Python/sphinx/index.rst +++ b/src/Readuct/Python/sphinx/index.rst @@ -6,7 +6,8 @@ User Manual The user manuals are available for download as PDF files: -- :download:`Download Latest User Manual PDF ` +- :download:`Download Latest User Manual PDF ` +- :download:`Download User Manual PDF v4.0.0 ` - :download:`Download User Manual PDF v3.0.0 ` - :download:`Download User Manual PDF v2.0.0 ` - :download:`Download User Manual PDF v1.0.0 ` diff --git a/src/Readuct/Readuct/ElementaryStepOptimization/ElementaryStepOptimizer.h b/src/Readuct/Readuct/ElementaryStepOptimization/ElementaryStepOptimizer.h index a7eb6b6..f30cea6 100644 --- a/src/Readuct/Readuct/ElementaryStepOptimization/ElementaryStepOptimizer.h +++ b/src/Readuct/Readuct/ElementaryStepOptimization/ElementaryStepOptimizer.h @@ -43,7 +43,7 @@ class ElementaryStepOptimizerBase { * * @return int The final number of optimization cycles carried out. */ - virtual int optimize() = 0; + virtual int optimize(Core::Log& log) = 0; /** * @brief Get the current reaction profile (at the end of an optimization, this is the final reaction profile). @@ -126,7 +126,7 @@ class ElementaryStepOptimizer : public ElementaryStepOptimizerBase { * * @return int The final number of optimization cycles carried out. */ - int optimize() override { + int optimize(Core::Log& log) override { // Preparations EnergiesAndGradientsAlongSpline valuesAlongSpline; const auto& elements = _profile.getMolecularSpline().getElements(); @@ -171,7 +171,7 @@ class ElementaryStepOptimizer : public ElementaryStepOptimizerBase { }; // Optimize - auto cycles = optimizer.optimize(variablesVector, update, check); + auto cycles = optimizer.optimize(variablesVector, update, check, log); return cycles; } diff --git a/src/Readuct/Tests/AppTests/test_readuct.py b/src/Readuct/Tests/AppTests/test_readuct.py index b1e7052..ad28bd7 100644 --- a/src/Readuct/Tests/AppTests/test_readuct.py +++ b/src/Readuct/Tests/AppTests/test_readuct.py @@ -63,7 +63,7 @@ def __parse_output(self): self.energy = float(splitted_output[index].split()[4]) if line == ' # cm^-1': self.lowest_frequency = float(splitted_output[index + 1].split()[1]) - if self._input['tasks'][0]['type'] in ['afir', 'geoopt', 'ts', 'irc', 'nt']: + if self._input['tasks'][0]['type'] in ['afir', 'geoopt', 'ts', 'irc', 'nt', 'nt2']: self.optimized_structure = self.__load_xyz_file('default_system_output/default_system_output.xyz') if self._input['tasks'][0]['type'] in ['bspline']: self.interpolated_trajectory = self.__load_trajectory( @@ -145,7 +145,7 @@ def test_wrong_input(self): 'name': 'default_system', 'method_family': 'PM6', 'settings': {'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7} + 'density_rmsd_criterion': 1e-7} } ], 'tasks': [ @@ -167,7 +167,7 @@ def test_wrong_input(self): 'method_family': 'PM6', 'wrong_keyword': 'something', 'settings': {'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7} + 'density_rmsd_criterion': 1e-7} } ], 'tasks': [ @@ -187,7 +187,7 @@ def test_wrong_input(self): 'name': 'default_system', 'method_family': 'PM6', 'settings': {'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7} + 'density_rmsd_criterion': 1e-7} } ], 'tasks': [ @@ -214,7 +214,7 @@ def test_wrong_input(self): {'input': ['default_system'], 'type': 'sp', } - ] + ] } calculation = Calculation(wrong_system_input) @@ -234,25 +234,48 @@ def test_wrong_input(self): 'type': 'sp', 'settings': {'stop_on_error': False} } - ] + ] } calculation = Calculation(wrong_system_input) success = calculation.run() assert success + def test_cellopt_with_unavailable_stress(self): + inp = {'systems': [ + {'path': os.path.join(os.path.dirname(os.path.realpath(__file__)), 'h2o2.xyz'), + 'name': 'default_system', + 'method_family': 'PM6', + 'settings': {'self_consistence_criterion': 1e-7, + 'density_rmsd_criterion': 1e-7} + } + ], + 'tasks': [ + {'input': ['default_system'], + 'type': 'opt', + 'settings': {'unitcelloptimizer': 'bfgs'} + } + ], + } + + calculation = Calculation(inp) + with self.assertRaises(subprocess.CalledProcessError) as context: + success = calculation.run() + def test_single_point_task(self): default_input = {'systems': [ {'path': os.path.join(os.path.dirname(os.path.realpath(__file__)), 'h2o2.xyz'), 'name': 'default_system', 'method_family': 'PM6', 'settings': {'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7} + 'density_rmsd_criterion': 1e-7} } ], 'tasks': [ {'input': ['default_system'], 'type': 'sp', + 'settings': {'silent_stdout_calculator': False, + 'spin_propensity_check': 1} } ] } @@ -271,7 +294,7 @@ def test_single_point_task(self): 'spin_multiplicity': 2, 'spin_mode': 'unrestricted', 'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7 + 'density_rmsd_criterion': 1e-7 } } ], @@ -294,7 +317,7 @@ def test_hessian_task(self): 'method_family': 'PM6', 'settings': { 'self_consistence_criterion': 1e-8, - 'density_rmsd_criterion' : 1e-8 + 'density_rmsd_criterion': 1e-8 } } ], @@ -340,7 +363,7 @@ def test_structure_optimization_task(self): 'name': 'default_system', 'method_family': 'PM6', 'settings': {'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7} + 'density_rmsd_criterion': 1e-7} } ], 'tasks': [ @@ -369,7 +392,7 @@ def test_single_atom_opt_task(self): 'name': 'default_system', 'method_family': 'PM6', 'settings': {'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7, + 'density_rmsd_criterion': 1e-7, 'molecular_charge': -1} } ], @@ -379,7 +402,7 @@ def test_single_atom_opt_task(self): 'type': 'geoopt', 'settings': {'bfgs_use_gdiis': False} } - ] + ] } default_calculation = Calculation(default_input) success = default_calculation.run() @@ -394,7 +417,7 @@ def test_transition_state_optimization_task(self): 'method_family': 'PM6', 'settings': { 'self_consistence_criterion': 1e-8, - 'density_rmsd_criterion' : 1e-8 + 'density_rmsd_criterion': 1e-8 } } ], @@ -431,7 +454,7 @@ def test_transition_state_optimization_task_dimer(self): 'method_family': 'PM6', 'settings': { 'self_consistence_criterion': 1e-8, - 'density_rmsd_criterion' : 1e-8 + 'density_rmsd_criterion': 1e-8 } } ], @@ -493,7 +516,7 @@ def test_afir_task(self): 'method_family': 'PM6', 'settings': {'molecular_charge': -2, 'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7} + 'density_rmsd_criterion': 1e-7} } ], 'tasks': [ @@ -526,7 +549,7 @@ def test_afir_task_2(self): 'settings': { 'molecular_charge': -1, 'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7 + 'density_rmsd_criterion': 1e-7 } } ], @@ -586,7 +609,7 @@ def test_nt_task(self): 'settings': { 'molecular_charge': -1, 'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7 + 'density_rmsd_criterion': 1e-7 } }], 'tasks': [{ @@ -617,6 +640,42 @@ def test_nt_task(self): # The following check yields 2.3 in release mode, but 2.5 in debug mode self.assertAlmostEqual(actual, 2.4, places=0) + def test_nt2_task(self): + default_input = {'systems': [{ + 'path': os.path.join(os.path.dirname(os.path.realpath(__file__)), 'sn2_start.xyz'), + 'name': 'default_system', + 'method_family': 'dftb3', + 'settings': { + 'molecular_charge': -1, + 'self_consistence_criterion': 1e-7, + 'density_rmsd_criterion': 1e-7 + } + }], + 'tasks': [{ + 'input': ['default_system'], + 'output': ['default_system_output'], + 'type': 'nt2', + 'settings': { + 'nt_associations': [0, 1], + 'nt_total_force_norm': 0.02, + 'nt_coordinate_system': 'cartesian', + 'sd_factor': 0.5, + 'convergence_max_iterations': 1000, + 'convergence_attractive_stop': 0.9, + }, + }] + } + + default_calculation = Calculation(default_input) + success = default_calculation.run() + self.assertTrue(success) + loaded = default_calculation.optimized_structure + actual = np.linalg.norm(np.array(loaded[0][1:]) - np.array(loaded[1][1:])) + self.assertAlmostEqual(round(actual, 2), 2.40, places=2) + actual = np.linalg.norm(np.array(loaded[5][1:]) - np.array(loaded[1][1:])) + # The following check yields 2.3 in release mode, but 2.5 in debug mode + self.assertAlmostEqual(actual, 2.4, places=0) + def test_irc_task(self): default_input = {'systems': [ {'path': os.path.join(os.path.dirname(os.path.realpath(__file__)), 'h2o2_ts.xyz'), @@ -624,7 +683,7 @@ def test_irc_task(self): 'method_family': 'PM6', 'settings': { 'self_consistence_criterion': 1e-8, - 'density_rmsd_criterion' : 1e-8 + 'density_rmsd_criterion': 1e-8 } } ], @@ -663,7 +722,7 @@ def test_bsplines_task(self): 'settings': { 'molecular_charge': -1, 'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7 + 'density_rmsd_criterion': 1e-7 } }, {'path': os.path.join(os.path.dirname(os.path.realpath(__file__)), 'end.xyz'), @@ -672,7 +731,7 @@ def test_bsplines_task(self): 'settings': { 'molecular_charge': -1, 'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7 + 'density_rmsd_criterion': 1e-7 } } ], @@ -878,12 +937,12 @@ def test_bsplines_interpolation_task(self): {'path': os.path.join(os.path.dirname(os.path.realpath(__file__)), 'h2_1.xyz'), 'name': 'start', 'method_family': 'PM6', - 'settings': {'self_consistence_criterion': 1e-7, 'density_rmsd_criterion' : 1e-7} + 'settings': {'self_consistence_criterion': 1e-7, 'density_rmsd_criterion': 1e-7} }, {'path': os.path.join(os.path.dirname(os.path.realpath(__file__)), 'h2_2.xyz'), 'name': 'end', 'method_family': 'PM6', - 'settings': {'self_consistence_criterion': 1e-7, 'density_rmsd_criterion' : 1e-7} + 'settings': {'self_consistence_criterion': 1e-7, 'density_rmsd_criterion': 1e-7} } ], 'tasks': [ @@ -916,12 +975,12 @@ def test_bsplines_interpolation_task(self): {'path': os.path.join(os.path.dirname(os.path.realpath(__file__)), 'h2_1.xyz'), 'name': 'start', 'method_family': 'PM6', - 'settings': {'self_consistence_criterion': 1e-7, 'density_rmsd_criterion' : 1e-7} + 'settings': {'self_consistence_criterion': 1e-7, 'density_rmsd_criterion': 1e-7} }, {'path': os.path.join(os.path.dirname(os.path.realpath(__file__)), 'h2_3.xyz'), 'name': 'end', 'method_family': 'PM6', - 'settings': {'self_consistence_criterion': 1e-7, 'density_rmsd_criterion' : 1e-7} + 'settings': {'self_consistence_criterion': 1e-7, 'density_rmsd_criterion': 1e-7} } ], 'tasks': [ @@ -961,7 +1020,7 @@ def test_bsplines_task_2(self): 'settings': { 'molecular_charge': -2, 'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7 + 'density_rmsd_criterion': 1e-7 } }, {'path': os.path.join(os.path.dirname(os.path.realpath(__file__)), 'end2.xyz'), @@ -970,7 +1029,7 @@ def test_bsplines_task_2(self): 'settings': { 'molecular_charge': -2, 'self_consistence_criterion': 1e-7, - 'density_rmsd_criterion' : 1e-7 + 'density_rmsd_criterion': 1e-7 } } ],