Coding with Theory Objects

AaronTools.theory.Theory() is used to designate computational options when building input files for different QM packages. It is designed to keep much of the code for writing quantum chemistry input files similar across different software packages to avoid having to dig into the nuances of the corresponding input file formats or worry about specific names for DFT functionals (e.g. PBE0 vs PBE1PBE or M06-2X vs M062X) or basis sets (def2-TZVP vs def2tzvp) in different packages.

Here, we will cover the basics of creating and using a Theory object to write input files for popular QM packages.

To have a functional Theory object, at a minimum you need to define the method, basis (unless using a semi-empirical method), and job_type.

Optionally, you can also specify:

• charge - overall charge

• multiplicity - multiplicity

• processors - allocated cores

• memory - allocated RAM

• empirical_dispersion - Grimme D2, D3, etc.

• solvent - continuum solvent model and solvent

• grid - integration grid

When writing an input file, additional keywords can be passed to AaronTools.geometry.Geometry.write() that specify any other options (often program-specific). These keywords are described in Theory parameters. See Additional Keywords.

For method, basis, job_type, empirical_dispersion, and grid you can either explicitly construct the corresponding object or provide a keyword and let Theory automatically construct the required object.

The latter approach is much simpler, but the former provides more control. This control is most often needed for basis (e.g. when using mixed basis sets or ECPs, etc).

For specifying solvent and continuum solvent model, you must use the ImplicitSolvent class object.

Building Basic Input Files

As a first example, the following will build a Gaussian input file called benzene.com that does a B3LYP/def2-SVP optimization of the coordinates in benzene.xyz followed by the computation of vibrational frequencies

from AaronTools.geometry import Geometry
from AaronTools.theory import *

geom = Geometry('benzene.xyz')

method = Theory(
    method="B3LYP",
    basis="def2-SVP",
    job_type=[OptimizationJob(), FrequencyJob()]
)
outfile = "benzene.com"
geom.write(outfile=outfile, theory=method)

Note that the resulting Gaussian input file will use def2svp for the basis set, even though we specified def2-SVP.

Equivalently, we could have done the following, where we explicitly build Method and BasisSet objects:

from AaronTools.geometry import Geometry
from AaronTools.theory import *

geom = Geometry('benzene.xyz')

method = Theory(
    method=Method("B3LYP"),
    basis=BasisSet("def2-SVP"),
    job_type=[OptimizationJob(), FrequencyJob()]
)
outfile = "benzene.com"
geom.write(outfile=outfile, theory=method)

If a molecule has a charge and multiplicity other than 0 and 1, we need to pass that to Theory:

method = Theory(
    method="B3LYP",
    charge=1,
    multiplicity=2,
    basis="def2-SVP",
    job_type=[OptimizationJob(), FrequencyJob()]
)

Similarly, if we want to use B3LYP-D3, instead of B3LYP, we can specify empirical_dispersion="D3"

method = Theory(
    method="B3LYP",
    charge=1,
    multiplicity=2,
    empirical_dispersion="D3",
    basis="def2-SVP",
    job_type=[OptimizationJob(), FrequencyJob()]
)

By changing the extension of the file being written, the corresponding format and keyword changes (def2-svp vs def2svp, etc) are automatically handled.

For example, the example below will write (essentially) equivalent input files for Gaussian, ORCA, and Psi4.

from AaronTools.geometry import Geometry
from AaronTools.theory import *

geom = Geometry('benzene.xyz')

method = Theory(
    method="B3LYP",
    basis="def2-SVP",
    job_type=[OptimizationJob(), FrequencyJob()]
)
for outfile in ["gaussian.com", "ORCA.inp", "psi4.in"]:
    geom.write(outfile=outfile, theory=method)

The above examples will all be in the gas phase. To instead use a continuum solvent model (e.g. PCM, SMD, etc.) we create and then use an ImplicitSolvent object (here requesting the SMD model and water as solvent):

from AaronTools.geometry import Geometry
from AaronTools.theory import *

geom = Geometry('benzene.xyz')

method = Theory(
    method="B3LYP",
    basis="def2-SVP",
    solvent=ImplicitSolvent("smd", "water"),
    job_type=[OptimizationJob(), FrequencyJob()]
)

for outfile in ["gaussian.com", "ORCA.inp", "psi4.in"]:
    geom.write(outfile=outfile, theory=method)

Job Types

There are currently seven job types in the theory package:

• AaronTools.theory.job_types.OptimizationJob()

• AaronTools.theory.job_types.FrequencyJob()

• AaronTools.theory.job_types.SinglePointJob()

• AaronTools.theory.job_types.ForceJob()

• AaronTools.theory.job_types.ConformerSearchJob()

• AaronTools.theory.job_types.TDDFTJob()

• AaronTools.theory.job_types.NMRJob()

A single JobType can be given to a Theory. If multiple JobType instances are given as list, the job-related information will appear in the order it appears in the list. For example, above we used job_type=[OptimizationJob(), FrequencyJob()] to specify a geometry optimization followed by vibrational frequencies.

However, if we instead did

job_type = [FrequencyJob(), OptimizationJob()]

then any Psi4 input file constructed using the corresponding Theory object will request frequencies before the optimization. Other programs are not sensitive to the order these jobs will appear in the input file.

Many of these job types take additional arguments (click the links above to see the options). For example, for a transition state optimization you need to specify OptimizationJob(transition_state=True).

JobType Abbreviations

There are shorthand versions of many common job types. For instance, to request numerical frequencies you could do

job_type = FrequencyJob(numerical=True)

or, equivalently,

job_type = "freq.num"

Similarly,

job_type = OptimizationJob(transition_state=True)

is equivalent to

job_type = "opt.ts"

See AaronTools.theory.job_types.job_from_string() for more information.

Constrained Optimizations

If we want to do a constrained optimization, we need to do a little more work. For example, suppose we have an AaronTools Geometry (probably not benzene) called geom and we want to write an input file for an optimization with a constraint on the distance between atoms 1 and 4. Constraints are passed to OptimizationJob() as a dictionary, with the keys corresponding to the types of constraints (bonds, angles, torsions, etc). Each entry in the dictionary is a list of lists of AaronTools Atoms. In our case, we are constraining a distance (‘bond’) so need to supply a list of a list of two atoms, which is most easily built using Geometry.find():

constraints = {}
constraints["bonds"] = [geom.find("1,4")]

Now we can pass this constraint dictionary to OptimizationJob():

method = Theory(
    method="B3LYP",
    basis="def2-SVP",
    job_type=OptimizationJob(constraints=constraints)
)

An input file written using this Theory object will include this geometric constraint, formatted properly for the correspinding QM package.

To add more constraints we simply append more pairs (or triples for an angle, quadruples for a torsion, etc) to the corresponding entry in the constraints dictionary. The following (silly) example will constrain distances 1-4 and 7-11, angle 2-3-5, and torsion 1-2-3-4:

constraints = {}
constraints["bonds"] = [geom.find("1,4"), geom.find("7,11")]
constraints["angles"] = [geom.find("2,3,5")]
constraints["torsions"] = [geom.find("1,2,3,4")]

method = Theory(
    method="B3LYP",
    basis="def2-SVP",
    job_type=OptimizationJob(constraints=constraints)
)

Alternatively, we can use Finders to build our constraints dictionary. For example, suppose we want to optimize the structure of benzene with all C-H bonds frozen. We could look up the atom numbers for each carbon and the associated H, or we can loop over all H atoms and use AaronTools.finders.BondedTo() to figure out the connected C atom:

from AaronTools.finders import BondedTo

constraints = {}
constraints["bonds"] = []

for H in geom.find('H'):
    bondedC = geom.find(BondedTo(H))[0]
    constraints["bonds"].append([bondedC, H])

method = Theory(
    method="B3LYP",
    basis="def2-SVP",
    job_type=OptimizationJob(constraints=constraints)
)

Finer Control

If you need more control over one or more of these objects you can explicitly define various objects and pass these to Theory. This is most likely to occur for BasisSet, for example, when working with mixed basis sets and/or ECPs.

The various objects that can be passed to Theory are discussed below.

Method Class

AaronTools.theory.Method() is used to keep method keywords the same across different formats. As an example:

from AaronTools.theory import Method

pbe0 = Method("PBE0")

When used to write a Gaussian input file, this Method will use the Gaussian keyword for PBE0 (PBE1PBE).

Method also takes a is_semiempirical argument:

rm1 = Method("RM1", is_semiempirical=True)

For Gaussian and ORCA input files, using a semi-empirical method will cause basis set information to be omitted.

SAPTMethod

AaronTools.theory.SAPTMethod() is a subclass of Method that is specific for SAPT jobs. When used to make a Psi4 input file, the molecule will be split into monomers, which are specified by the components attribute of the Geometry instance.

sapt0 = SAPTMethod("sapt0")

See SAPT Calculations for an example.

Basis Sets

The AaronTools.theory.BasisSet() object is a collection of AaronTools.theory.Basis() and (optionally) AaronTools.theory.ECP() objects.

The second argument given to each Basis determines which elements that basis applies to. By default, a Basis applies to all elements while an ECP applies to any transition metal.

For example, suppose we have some Pt carbonyl complex. To build a BasisSet object for a calculation in which we use LANL2DZ basis set and ECP on Pt and 6-31G(d) on C and O, we could do

from AaronTools.theory import Basis, ECP, BasisSet
basis = BasisSet(
    [
        Basis("6-31G(d)", ["C", "O"]),
        Basis("LANL2DZ", "Pt")
    ],
    [ECP("LANL2DZ")]
)

Alternatively, we can use Finders to automatically build lists of elements:

from AaronTools.theory import Basis, ECP, BasisSet
from AaronTools.finders import AnyTransitionMetal, AnyNonTransitionMetal

basis = BasisSet(
    [
        Basis("6-31G(d)", AnyNonTransitionMetal()),
        Basis("LANL2DZ", AnyTransitionMetal()),
    ],
    [ECP("LANL2DZ")]
)

Finally, the aux_type keyword is used for ORCA and Psi4 input files to specify auxiliary basis sets.

from AaronTools.theory import Basis, ECP, BasisSet
from AaronTools.finders import AnyTransitionMetal, AnyNonTransitionMetal
basis = BasisSet(
    [
        Basis("cc-pVTZ", AnyNonTransitionMetal()),
        Basis("cc-pVTZ", AnyNonTransitionMetal(), aux_type='C'),
        Basis("cc-pVTZ-PP", AnyTransitionMetal()),
        Basis("cc-pVTZ-PP", AnyTransitionMetal(), aux_type='C')
    ],
    [ECP("SK-MCDHF-RSC")]
)

Any of these BasiSet objects can then be passed to a Theory object. For example, the following will write a Gaussian input file TM_complex.com for an optimization + frequency job at the M06-2X/6-31G(d)/LANL2DZ level of theory for any transition metal complex in TM_complex.xyz:

from AaronTools.geometry import Geometry
from AaronTools.theory import *
from AaronTools.finders import AnyTransitionMetal, AnyNonTransitionMetal

geom = Geometry('TM_complex.xyz')

basis = BasisSet(
    [
        Basis("6-31G(d)", AnyNonTransitionMetal()),
        Basis("LANL2DZ", AnyTransitionMetal()),
    ],
    [ECP("LANL2DZ")]

method = Theory(
    method="M062X",
    basis=basis,
    job_type=[OptimizationJob(), FrequencyJob()]
)
outfile = "TM_complex.com"
geom.write(outfile=outfile, theory=method)

Empirical Dispersion

AaronTools.theory.emp_dispersion.EmpiricalDispersion() keeps specifying dispersion corrections consistent across different input file formats.

from AaronTools.theory import EmpiricalDispersion

disp = EmpiricalDispersion("Grimme D2")

# The following are equivalent:
disp = EmpiricalDispersion("Grimme D2")
disp = EmpiricalDispersion("GD2")
disp = EmpiricalDispersion("D2")
disp = EmpiricalDispersion("-D2")

Some dispersion methods are not available in all QM software programs. Check the get_gaussian, get_orca, etc. methods of the EmpiricalDispersion class (or the respective manuals) for acceptable dispersion methods.

Continuum Solvent Model

The AaronTools.theory.ImplicitSolvent() object is how you request a continuum solvent model in a Theory object:

from AaronTools.theory import ImplicitSolvent

solvent = ImplicitSolvent("PCM", "toluene")

Integration Grid

As with other objects in the AaronTools.theory package, the AaronTools.theory.IntegrationGrid() object is a way to specify grids in a similar manner across different file formats.

It’s important to note that different programs use different types of grids. This, combined with varied grid pruning algorithms, mean that getting exactly equivalent grids in two QM programs is nearly impossible. If you use a keyword from one program to make an input file for a different program, IntegrationGrid will at least try to specify an equivalent grid.

from AaronTools.theory import IntegrationGrid

grid = IntegrationGrid("SuperFineGrid")

Gaussian, ORCA, and Psi4 all have different ways of specifying integration grids. Gaussian and ORCA have grid keywords. When using an ORCA grid keyword to write a Gaussian input file, IntegrationGrid will try to approximate the ORCA grid’s density. Psi4 specifies grid density by supplying a number of radial and angular points. Gaussian allows a similar specification. These can be specified as a string of the format "(radial, angular)". As an example,

grid = IntegrationGrid("(99, 590)")

This grid can be used with Gaussian and Psi4, and should give similar results (down to grid pruning and other algorithmic differences). If you’re going to write an ORCA input file with this grid, the number of radial points is set indirectly with the IntAcc option. IntAcc will be set for the number of radial points in the 2nd row of the periodic table.

Additional Keywords

Additional program options are often program-specific and are passed to AaronTools.geometry.Geometry.write() differently depending on the QM package and the location where the additional options are required. These keywords are described in Theory parameters.

For example, in Constrained Optimizations we added constraints using the constraints option in OptimizationJob(). Alternatively, we can directly write data to the constraints section of a Gaussian input file using GAUSSIAN_CONSTRAINTS. For instance, we can write the constraints from Constrained Optimizations by modifying the geom.write line:

geom.write(outfile=outfile, theory=method, GAUSSIAN_CONSTRAINTS = "B 1 4 F\nB 7 11 F\nA 1 2 3 F\nD 1 2 3 4 F")

The advantage of building a constraints dictionary and passing that to OptimizationJob() is that you can more easily switch to a different QM package.

Some of these additional keywords take a dictionary. GAUSSIAN_ROUTE provides a nice example. As noted above, you can requiest a TS optimization by passing transition_state=True to OptimizationJob(). However, what if you also want to include noeigen as an option to opt? In other words, by using transition_state=True the route section will include opt=(ts,CalcFC), but we want to add noeigen to the list of opt options. We can do this by defining a dictionary with key opt and value noeigen and noeigen will automatically be added to the list of options under opt:

route = {"opt": "noeigen"}
geom.write(outfile=outfile, theory=method, GAUSSIAN_ROUTE=route)

For route entries with no options (e.g. nosym) you simply provide the key but an empty value:

route = {"nosym": ""}
geom.write(outfile=outfile, theory=method, GAUSSIAN_ROUTE=route)