6. Calculators

Exploring potential energy surfaces (PESs) requires calculation of energies and its derivatives (gradients, Hessian matrices).

For testing purposes and method development, pysisyphus implements various analytical 2d potentials, as they allow fast evaluations of the aforementioned quantities. Actual (production) calculations are carried out by wrapping existing quantum chemistry codes (ORCA, TURBOMOLE, Gaussian etc.) and delegating the required calculations to them. Pysisyphus generates all necessary inputs, executes the QC code and parses their output for the requested quantities.

Furthermore pysisyphus provides several "meta"-calculators which wrap other (actual) calculators, to modify calculated energies and forces. Examples for this are the Dimer calculator, used for carrying out transition state searches or the ONIOM calculator, allowing multi-level calculations comprising different levels of theory.

External forces, e.g. a restraining spherical potential or harmonic restraints on primitive internal coordinates (stretches, bends, torsion) can be applied with ExternalPotential.

6.1. YAML input

Possible keywords for the YAML input can be derived from inspecting the possible arguments of the Calculator base-class (see below) and the possible arguments of the respective calculator, e.g. ORCA or XTB.

The most commonly used keywords, derived from the Calculator baseclass are mem, handling the requested memory per core in MB, pal, handling the number of requested CPU cores, charge, the total charge of the system and mult, the systems multiplicity.

For (excited state, ES) calculations carried out with calculators derived from OverlapCalculator additional keywords are possible. The most common keywords controlling ES calculations are track, activating ES tracking, ovlp_type, selecting the tracking algorithm and ovlp_with, handling the selection of the reference cycle.

An example input highlighting the most important keywords is shown below.

geom:
 [... omitted ...]
calc:
 type: orca                     # Calculator type, e.g. g09/g16/openmolcas/
                                # orca/pyscf/turbomole/dftb+/mopac/psi4/xtb

 pal: 1                         # Number of CPU cores
 mem: 1000                      # Memory per core
 charge: 0                      # Charge
 mult: 1                        # Multiplicity
 # Keywords for ES calculations
 track: False                   # Activate ES tracking
 ovlp_type: tden                # Tracking algorithm
 ovlp_with: previous            # Reference cycle selection
 # Additional calculator specific keywords
 [... omitted ...]
opt:
 [... omitted ...]

6.2. Calculator base classes

6.3. OverlapCalculator base class

6.4. Calculators with Excited state capabilities

6.4.1. Gaussian09

6.4.2. Gaussian16

6.4.3. OpenMolcas

Pysisyphus currently supports energy and gradient calculations utilizing the &rasscf and/or the &mcpdft sections. Neither analytical nor numerical Hessians are yet implemented for the OpenMolcas-calculator.

Two keywords are always required: inporb and basis, with the former pointing to a .RasOrb file and the latter containing the selected atomic orbital basis, e.g., ano-rcc-vdzp. Additional input for the &gateway, &rasscf and &mcpdft sections can be given via the respective keyword(s).

Due to restrictions of the current design, simple keywords that don't take further arguments as cmsi in &rasscf or grad and mspdft in &mcpdft still must be given with a trailing colon. See below for an example.

Listing 6.1 Optimization using compressed-multi-state PDFT.
geom:
 type: redund
 fn: |
  4
  
  C	-1.0398336639	 0.0	    0.0
  S	 0.6002429216	 0.0	    0.0
  H -1.6321592382	-0.94139864	0.0
  H	-1.6321592382	 0.94139864	0.0
calc:
 type: openmolcas
 basis: cc-pvdz
 charge: 0
 mult: 1
 inporb: /home/johannes/tmp/359_cmspdft/359_cmspdft.RasOrb
 rasscf:
  ciroot: 3 3 1
  mdrlxroot: 2
  cmsi:
 mcpdft:
  ksdft: t:pbe
  grad:
  mspdft:
opt:
 thresh: gau

6.4.4. ORCA 4.2.1 / 5.0.1

6.4.5. PySCF 1.7.6

6.4.6. Turbomole 7.x

Pysisyphus does not implement a wrapper for define, so the user has to manually prepare a directory containing a valid control file. An automated define wrapper, restricted to ground state functionality, is available via the QCEngine project, to which I contributed the Turbomole harness.

Care should be taken to include only the minimum amount of necessary files in the control_path directory, e.g., (auxbasis, basis, control, coord, mos) for a closed-shell calculation using RI. A gradient file must not be present in control_path, as well as other subdirectories and files with .out extension. The coord file, while not strictly required, should be kept too, to facilitate testing of the setup with standalone Turbomole.

It may be a good idea to pre-converge the calculation in control_path, to see if the setup is correct and actually works. Resulting files like energy, statistics can be deleted; mos should be kept, as the converged MOs are reused in pysisyphus.

If an excited-state optimization is desired, care has to be taken, to include $exopt [n] for TD-DFT/TDA or the geoopt state=([n]) (ricc2)! Tracking of excited states is currently possible for closed shell egrad and ricc2 calculations.

The current implementation was tested against Turbomole 7.4.1 and QCEngine 0.19.0. Please see examples/complex/11_turbomole_gs_tsopt for a full example where Turbmole is utilized in a growing string calculation. The same example, using QCEngine, is found in examples/complex/12_qcengine_turbomole_gs_tsopt. MOs are not reused with the QCEngine calculator, so the native pysisyphus calculator is probably faster.

6.4.7. DFTB+ 20.x

6.5. Calculators with Ground state capabilities

6.5.1. MOPAC 2016

6.5.2. Psi4

6.5.3. QCEngine

6.5.4. XTB 6.x

6.5.5. Dalton

6.5.6. OpenBabel

6.5.7. CFOUR

6.6. Meta (wrapping) Calculators

6.6.1. ExternalPotential

6.6.2. Restraint

# General input structure for restraints
calc:
 type: ext
 # Multiple potentials could be specified here as a list
 potentials:
   # Primitive internal coordinate restraint
   - type: restraint
     # List of restraints; could also be multiple restraints. Each restraint is given as
     # list of 2 or 3 items.
     #
     # The first item always specifies an internal coordinate,
     # whereas the second argument is a force constant (in atomic units; actual units
     # depend on the coordinate). Optionally a reference value (third argument) can be
     # given. If omitted, the initial coordinate value is used as reference value.
     restraints: [[[BOND, 0, 1], 10, 3.0]]
     # The commented out input below would restrain the bond at its initial value.
     #restraints: [[[BOND, 0, 1], 10]]
     # Multiple restraints are specified as given below.
     #restraints: [[[BOND, 0, 1], 10], [[BEND, 0, 1, 2], 1.0]]
calc:
 type: [actual calculator that is wrapped by ExternalPotential]

6.6.3. HarmonicSphere

6.6.4. LogFermi

6.6.5. RMSD

6.6.6. DFT-D3

Method to add DFT-D3 dispersion corrections as an external potential via the program developed by the Grimme group <https://www.chemie.uni-bonn.de/grimme/de/software/dft-d3/get_dft-d3>.

This is for use with calculators that do not natively provide D3 corrections (e.g. OpenMolcas). Usage mirrors that of other external potentials, with an example given below.

# General input structure for restraints
calc:
 type: ext
 # Multiple potentials could be specified here as a list
 potentials:
   # Add atom-pairwise D3 dispersion correction as a differentiable, external potential
   - type: d3
     # Functional is specified in TURBOMOLE format, all lower case.
     functional: pbe
     # Optional Becke-Johnson damping, default false, recommended true
     bjdamping: true

calc:
 type: [actual calculator that is wrapped by ExternalPotential]

6.6.7. AFIR

6.6.8. ONIOM

6.6.9. Dimer

6.7. Pure Python calculators

6.7.1. Sympy 2D Potentials

6.7.2. Lennard-Jones

6.7.3. FakeASE

6.7.4. TIP3P