Example Python Scripts for Common Tasks

Many routine tasks can be handled using BASH scripts based on the AaronTools command line scripts (CLSs). See Example BASH Scripts for Common Tasks. However, for more complex tasks you’ll need to use the AaronTools Python API. A few examples are provided here to get you started (you can also peruse the CLSs, although some are overly complicated for generality). More examples will be added over time to make this a sort of ‘cookbook’ of AaronTools Python scripts

Note that these scripts assume you have correctly installed AaronTools (see Installation Guide) and configured the Job Submission Templates to allow for job submission to a local queue.

Optimization and Frequencies

Although the task of submitting optimization + frequency jobs can be readily handled using a BASH script (see Example BASH Scripts for Common Tasks), we include a Python version here for completeness. Here we use ORCA as an example, but the script below would instead use Psi4 or Gaussian by simply changing the filename extension to .in or .com, respectfully!

The following script will submit wB97XD/def2-TZVP optimization and frequency jobs (requesting 8 cores/processors, a 12 hour walltime, and 12 GB of memory) on all XYZ files in the current directory:

# build ORCA input files for optimization + frequency calculation on all XYZ files in current directory
from os.path import splitext
from AaronTools.geometry import Geometry
from AaronTools.theory import *
from AaronTools.job_control import SubmitProcess
import glob

# build theory object
method = Theory(
    method="wb97xd",
    basis="def2-tzvp",
    job_type=[OptimizationJob(), FrequencyJob()],
    processors=8,
    memory=12,
)

# loop over XYZ files in current directory
for XYZ in glob.glob("*.xyz"):
    name, _ = splitext(XYZ)   # grab name without .xyz ending
    geom = Geometry(XYZ)      # build AaronTools geometry
    outfile = f"{name}.inp"   # filename with .inp extension for ORCA
    geom.write(outfile=outfile, theory=method) # make input file

    # create SubmitProcess object
    submit_process = SubmitProcess(outfile, 12, 8, 12)

    # submit job
    try:
        submit_process.submit()
    except Exception as e:
        warn("failed to submit %s: %s" % (outfile, str(e)))

This could be easily modified to use IUPAC names or SMILES.

Transition Metal Complexes

Alternatively, if we have a set of transition metal complexes, we need to build a AaronTools.theory.BasisSet() object to specify different basis sets/ECPs for the transition metal (See Coding with Theory Objects).

The following will build Gaussian input files to run optimizations and frequencies at the B3LYP-D3/6-31G(d)/LANL2DZ level of theory for all XYZ files in the current directory:

# build Gaussian input files for optimization + frequency calculation
# on all XYZ files in current directory 
# use LANL2DZ basis and ECP on any transition metal and 6-31G(d) on everything else

from os.path import splitext
from AaronTools.geometry import Geometry
from AaronTools.theory import *
from AaronTools.job_control import SubmitProcess
import glob

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

# build theory object
method = Theory(
    method="b3lyp",
    empirical_dispersion="D3",
    basis=basis,
    job_type=[OptimizationJob(), FrequencyJob()]
    processors=8,
    memory=12,
)

# loop over XYZ files in current directory
for XYZ in glob.glob("*.xyz"):
    name, _ = splitext(XYZ)   # grab name without .xyz ending
    geom = Geometry(XYZ)      # build AaronTools geometry
    outfile = f"{name}.com"   # filename with .com extension for Gaussian
    geom.write(outfile=outfile, theory=method) # make input file

    # create SubmitProcess object
    submit_process = SubmitProcess(outfile, 12, 8, 12)

    # submit job
    try:
        submit_process.submit()
    except Exception as e:
        warn("failed to submit %s: %s" % (outfile, str(e)))

After these are run to completion, we could submit B3LYP-D3/6-311+G(d,p)/LANL2DZ single point energies using the following:

# build Gaussian input files for Single Point calculations on all LOG files in 
# current directory 
# use LANL2DZ basis and ECP on any transition metal and 6-311+G(d,p) on everything else

from os.path import splitext
from AaronTools.geometry import Geometry
from AaronTools.theory import *
from AaronTools.job_control import SubmitProcess
import glob

# build basis object
basis = BasisSet(
    [
        Basis("6-311+G(d,p)", AnyNonTransitionMetal()),
        Basis("LANL2DZ", AnyTransitionMetal()),
    ],
    [ECP("LANL2DZ")]
)

# build theory object
method = Theory(
    method="b3lyp",
    empirical_dispersion="D3",
    basis=basis,
    job_type=SinglePointJob()
    processors=8,
    memory=12,
)

# loop over LOG files in current directory
for LOG in glob.glob("*.log"):
    name, _ = splitext(LOG)    # grab name without .log ending
    geom = Geometry(LOG)       # build AaronTools geometry
    outfile = f"{name}.sp.com" # build filename with .com extension for Gaussian
    geom.write(outfile=outfile, theory=method) # make input file

    # create SubmitProcess object
    submit_process = SubmitProcess(outfile, 12, 8, 12)

    # submit job
    try:
        submit_process.submit()
    except Exception as e:
        warn("failed to submit %s: %s" % (outfile, str(e)))

SAPT Calculations

In Symmetry Adapted Perturbation Theory (SAPT) we built a script to run SAPT calculations on the parallel stacked benzene dimer as a function of horizontal and vertical displacements. In that case, because we were using makeInput.py, we had to rely on Psi4 automatically partitioning the dimer into fragments. If more control is required over how the supermolecule is fragmented for SAPT computations, we need to build a Psi4 input file explicitly separating the molecule into componets.

This requires some small changes to how we build the geometry and theory objects. First, we need to define the components of Geometry to be a list of the individual monomers. Second, in the theory object we need to use SAPTMethod instead of Method. Finally, the charge and multiplicity need to be lists consisting of the charge/multiplicity for the whole system and each monomer.

For instance, the following will submit a Psi4 job to calculate the SAPT0/jun-cc-pVDZ energy on a dimer (read from dimer.xyz) but use the AaronTools function AaronTools.geometry.Geometry.get_monomers() to separate the monomers:

from AaronTools.component import Component
from AaronTools.geometry import Geometry
from AaronTools.theory import *
from AaronTools.job_control import SubmitProcess

# make Geometry object from XYZ file containing dimer
dimer = Geometry("dimer.xyz")

# for SAPT, the monomers need to be in the geometry's components attribute
dimer.components = []
for mon in dimer.get_monomers():
    dimer.components.append(Component(mon))

# Modified Theory object for SAPT computations
theory = Theory(
    charge=[0, 0, 0],
    multiplicity=[1, 1, 1],
    method=SAPTMethod("sapt0"),
    basis="jun-cc-pvdz",
    processors=4,
    memory=8,
    job_type="energy",
)

outfile='dimer.in'

# write Psi4 input then submit job
dimer.write(outfile=outfile, theory=theory)
submit_process = SubmitProcess(outfile, 12, 4, 8)
try:
    submit_process.submit()
except Exception as e:
    warn("failed to submit %s: %s" % (outfile, str(e)))

Splitting Up Multi-structure XYZ Files

Often, one needs to take an XYZ file containing multiple structures (e.g. the conformers generated by CREST) and split this into individual XYZ files for each structure. This can be done easily in AaronTools:

#!/usr/bin/env python3
# split a multistructure XYZ file (e.g. from Crest) into individual XYZ files

from AaronTools.geometry import Geometry
from AaronTools.fileIO import FileReader

fr = FileReader('crest_conformers.xyz', get_all=True)

# loop over conformers
for confnum, conf in enumerate(fr.all_geom):
    geom = Geometry(conf["atoms"])
    outfile = f"conf{confnum}.xyz"
    geom.write(outfile=outfile)

Analyzing XYZ Trajectories

AaronTools AaronTools.fileIO.FileReader() objects can be used to read all geometries from multi-structure XYZ files. This can be used, for example, to analyze distances, angles, etc. from an XYZ trajectory file.

The script below does this, reading traj.xyz and printing the distance between atoms 1 and 2 and the 2-1-3 angle for each step. This could be easily modified to instead only print every N steps, etc. or other geometric parameters.

#!/usr/bin/env python3
# print a table of selected distances and angles from an XYZ trajectory

from AaronTools.geometry import Geometry
from AaronTools.fileIO import FileReader
import numpy as np

fr = FileReader('traj.xyz', get_all=True)

# print header
print("  Step     R(1,2)   A(2,1,3)")

# loop over all XYZ files
for step, struc in enumerate(fr.all_geom):
    geom = Geometry(struc["atoms"])
    # get atoms 1, 2, 3
    a1, a2, a3 = geom.find(['1','2','3'])
    # get distance and angle
    dist = a1.dist(a2)
    angle = np.degrees(geom.angle(a2, a1, a3))
    # print data
    print(f"{step:5} {dist:10.3f} {angle:10.2f}")