Example BASH Scripts for Common Tasks

AaronTools command line scripts (CLSs, see List of Command Line Scripts) make it simple to write BASH scripts to automate many routine tasks. A few examples are provided here to get you started! More examples will be added over time to make this a sort of ‘cookbook’ of AaronTools BASH scripts If you have a BASH script that uses AaronTools that you’d like us to add send it our way (similarly, if there is something you are trying to do using AaronTools contact us and we can (possibly) help).

Note that these BASH 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

Perhaps the most common task in routine quantum chemistry applications is optimizing geometries for a set of molecules and computing harmonic vibrational frequencies. This can be done using AaronTools CLSs using any quantum chemistry software package. Here we use Gaussian as an example, but trivial changes to the scripts below would instead use Q-Chem, ORCA, or Psi4.

For instance, suppose we have a set of XYZ files containing initial geometries that we would like to optimize. These geometries could be for a set of molecules from some other source or conformations of a single molecule generated by CREST for example. The following script will submit wB97XD/def2-TZVP optimization and frequency jobs (requesting 8 cores/processors and 12 GB of memory) on all XYZ files in the current directory. It uses makeInput.py to make a Gaussian input file from each XYZ file and then calls jobSubmit.py to submit a job to the queue:

#!/bin/bash
# Perform wB97XD/def2TZVP optimizations/frequencies on all XYZ files 

for XYZ in `/bin/ls *.xyz`; do
    name=`basename $XYZ .xyz`  # get filename without .xyz
    makeInput.py $XYZ -m wB97XD -b "def2tzvp" -opt -freq -o $name.com -np 8 -mem 12
    jobSubmit.py $name.com -p 8 -m 12
done

If these XYZ files were different conformations of the same molecule, after all optimizations were complete we could then identify the unique ones using unique.py, for example.

If initial XYZ files are not available, we can instead build all molecules from a list of SMILES or IUPAC names (using fetchMolecule.py) and pipe this into makeInput.py to build input files for the optimization + frequency computations. For instance, suppose we have a file molecules that contains a list of IUPAC names:

#!/bin/bash
# Perform wB97XD/def2TZVP optimizations/frequencies on all molecules listed in 'molecules'

for mol in `cat molecules`; do
    safe_name=`echo $mol | tr -dc [:alnum:]` # strip non alphanumeric characters from name
    fetchMolecule.py -i "$mol" | makeInput.py -m wB97XD -b "def2tzvp" -opt -freq -o $safe_name.com -np 8 -mem 12
    jobSubmit.py $safe_name.com -p 8 -m 12
done

Because many IUPAC names contain spaces and/or other non-alphanumeric characters (e.g. 1,3-butadiene) we remove those to form the name of the input file using the BASH command tr -dc [:alnum:].

The above script can be easily extended to run optimizations/frequency jobs using different combinations of DFT functionals and basis sets:

#!/bin/bash
# Perform optimizations/frequencies on all molecules listed in 'molecules' using different methods

for func in wb97xd b3lyp m062x; do
    mkdir $func
    cd $func
    for basis in "def2tzvp" "6-31+G(d)" "cc-pVTZ"; do
        basis_name=`echo $basis | tr -dc [:alnum:]`
        mkdir $basis_name
        cd $basis_name
        for mol in `cat ../../molecules`; do
            safe_name=`echo $mol | tr -dc [:alnum:]` # strip non alphanumeric characters from name
            fetchMolecule.py -i "$mol" | makeInput.py -m $func -b $basis -opt -freq -o $safe_name.com -np 8 -mem 12
            jobSubmit.py $safe_name.com -p 8 -m 12
        done
        cd ..
    done
    cd ..
done

Checking and Rerunning Jobs

After running the above optimizations, it is necessary to check two things:

1. Did the computations all finish correctly?

2. Are all optimized structures energy minima (no imaginary frequencies)

The following script will check these. For jobs that did not generate an output file, the script just resubmits the original input file. For jobs that ran but did not finish, a new job is submitted using the last geometry from the previous attempt using the same level of theory used in the corresponding output file. For optimizations that ended with one or more imaginary vibrational frequency, the script displaces along the (lowest magnitude) imaginary vibrational mode using follow.py and re-optimizes. We’ll assume that we are working with a set of XYZ files, as in the first example above. The modifications to instead use a list of IUPAC names should be obvious. Changing the level of theory and/or quantum chemistry package should also be straightforward.

#!/bin/bash
# check all jobs and resubmit if needed, following imaginary frequency if present
method="wb97xd"
basis="def2tzvp"

for XYZ in `/bin/ls *.xyz`; do
    name=`basename $XYZ .xyz`
    # Check if output file exists
    if [ ! -e $name.log ]; then
        echo $name did not run! Trying again.
        jobSubmit.py $name.com -p 8 -m 12
    else
        # check if job finished
        if [ `printInfo.py -i FINISHED $name.log | grep -c True` == 1 ]; then
            # check if geometry optimized to a minimum
            if [ `printFreq.py -t neg $name.log | wc -l` == 1 ]; then
                echo $name finished!
            else
                echo $name optimized to a saddle point--following imaginary mode
                follow.py $name.log | makeInput.py -m $method -b "$basis" -opt -freq -o $name.com -np 8 -mem 12
                jobSubmit.py $name.com -p 8 -m 12
            fi
        else
            echo $name not finished! Restarting
            makeInput.py $name.log -m $method -b "$basis" -opt -freq -o $name.com -np 8 -mem 12
            jobSubmit.py $name.com -p 8 -m 12
        fi
    fi
done

Potential Energy Scans for Stacked Dimers

Density Functional Theory

AaronTools CLSs make it trivial to run scans over different coordinates for non-bonded dimers (e.g. see this recent paper). For instance, the script below will run wB97X-D/def2-TZVP single point energies on the benzene dimer as a function of x- and z-coordinates (with a fixed monomer geometry) from x = 0 to 5 A and z = 3 - 4 A. It assumes that there is a local file called benzene.xyz that contains the coordinates of benzene in the xy-plane that is oriented with a vertex along the x-axis.

#!/bin/bash
# scan over x, z coordinates for benzene dimer and calculate wB97X-D/def2-TZVP single point for each
# assumes there is a local file called 'benzene.xyz' containing benzene coordinates in xy-plane

for x in `seq 0 0.5 5`; do
    for z in `seq 3 0.2 4`; do
        name=x.$x.z.$z
        translate.py benzene.xyz -v $x 0 $z -o temp.xyz
        combineMonomers.py benzene.xyz temp.xyz | makeInput.py -m wb97xd -b def2tzvp -o $name.com -np 8 -mem 12
        jobSubmit.py $name.com -p 8 -m 12
    done
done

rm -f temp.xyz

Note that combineMonomers.py can only accept geometries from STDIN or from a list of files, so we save the shifted benzene geometry to temp.xyz before combining with benzene.xyz to build each dimer.

If you have an aversion to writing temporary files, we could instead combine the monomers first and then use translate.py to shift atoms 13-24:

#!/bin/bash
# scan over x, z coordinates for benzene dimer and calculate wB97X-D/def2-TZVP single point for each
# assumes there is a local file called 'benzene.xyz' containing benzene coordinates in xy-plane

for x in `seq 0 0.5 5`; do
    for z in `seq 3 0.2 4`; do
        name=x.$x.z.$z
        combineMonomers.py benzene.xyz benzene.xyz | translate.py -v $x 0 $z -t 13-24 | makeInput.py -m wb97xd -b def2tzvp -o $name.com -np 8 -mem 12
        jobSubmit.py $name.com -p 8 -m 12
    done
done

We could modify either of these to instead run the T-shaped benzene dimer by first rotating the second benzene by 90 degrees around the y-axis (e.g. rotate.py benzene.xyz -x y -a 90) and then piping this into translate.py -v $x 0 $z ....

We can also use AaronTools to gather the energies after making sure all of the single points ran to completion:

#!/bin/bash
# gather energies from dimer scan

for x in `seq 0 0.5 5`; do
    for z in `seq 3 0.2 4`; do
        name=x.$x.z.$z
        E=`printInfo.py $name.log -i scf_energy | tail -n1 | awk '{print $3}'`
        echo $x $z $E
    done
done

Symmetry Adapted Perturbation Theory (SAPT)

To instead run SAPT computations (using Psi4), we need to make some minor modifications. SAPT computations provide interaction energies between non-bonded molecules. As such, these computations require a dimer (for example) to be partitioned into individual components. In order to use makeInput.py to build a Psi4 SAPT input file we will exploit the fact that Psi4 can automatically detect components for SAPT computations using activate(auto_fragments()). Because this needs to appear in the Psi4 input file before the job, we pass this to makeInput.py using the --before-job option. Note that with Psi4 we need to explicitly specify that Psi4 calculate the energy (by passing -e to makeInput.py).

#!/bin/bash
# scan over x, z coordinates for benzene dimer and calculate SAPT0/jun-cc-pVDZ interaction energies
# assumes there is a local file called 'benzene.xyz' containing benzene coordinates in xy-plane

for x in `seq 0 0.5 5`; do
    for z in `seq 3 0.2 4`; do
        # name job by x and z value (e.g. x.0.0.z.3.0.com)
        name=x.$x.z.$z

        # move one benzene to (x, 0, z)
        translate.py benzene.xyz -v $x 0 $z -o temp.xyz
        
        # combine benzene.xyz and temp.xyz to form dimer then build Psi4 input file
        combineMonomers.py benzene.xyz temp.xyz | makeInput.py -e -m sapt0 -b jun-cc-pvdz --before-job 'activate(auto_fragments())' -o $name.in  -np 8 -mem 12

        # submit job
        jobSubmit.py $name.in -p 8 -m 12
    done
done

If Psi4 does not properly recognize your monomers, or you want more control over how a given complex is divided, you will need to use the Python API. See SAPT Calculations.