A simple script to run torsiondrives via qcengine

simple-td.py

# A script to quickly run a torsion drive using QCEngine torsiondrive/geometric openff toolkit and rdkit uff
from openff.toolkit.topology import Molecule, Atom, Bond
from qcelemental.models.common_models import Model
from qcelemental.models import DriverEnum
from qcelemental.models.procedures import TDKeywords, TorsionDriveInput, QCInputSpecification, OptimizationSpecification, TorsionDriveResult
import qcengine as qcng


# Make our molecule
print("Building molecule")
ethane: Molecule = Molecule.from_smiles("CC")
ethane.generate_conformers(n_conformers=1)


# Now we can find the indices of the torsion atoms
bond: Bond = ethane.find_rotatable_bonds()[0]
atom2: Atom = bond.atom1
atom3: Atom = bond.atom2
atom1 = [atom.molecule_particle_index for atom in atom2.bonded_atoms if atom != atom3][0]
atom4 = [atom.molecule_particle_index for atom in atom3.bonded_atoms if atom != atom2][0]
dihedrals = [(atom1, atom2.molecule_particle_index, atom3.molecule_particle_index, atom4)]
print("Dihedral found with indices ", dihedrals[0])

#Build the input schema for qcengine, here we will run with uff from rdkit
keywords = TDKeywords(dihedrals=dihedrals, grid_spacing=[15])
model = Model(method="uff", basis=None)
input_specification = QCInputSpecification(driver=DriverEnum.gradient, model=model)
input_molecule = ethane.to_qcschema()
optimisation_specification = OptimizationSpecification(
    procedure="geometric",
    # set the geometric keywords following the torsiondrivce CLI defaults
    keywords={
        "coordsys": "dlc",
        "maxiter": 300,
        "program": "rdkit",
        "enforce": 0.1,
        "reset": True,
        "qccnv": True,
        "epsilon": 0.0
    }
)
input_data = TorsionDriveInput(
    keywords=keywords,
    input_specification=input_specification,
    initial_molecule=input_molecule,
    optimization_spec=optimisation_specification
)

# run the torsiondrive
print("Running torsiondrive ...")
# Note this only uses a single input conformation, and runs all opts in serial 
# edit the local options to pass more resources to qcengine
result = qcng.compute_procedure(input_data=input_data, procedure="torsiondrive", raise_error=True, local_options={"ncores": 1, "memory": 1})
print("Torsiondrive complete, writing results to result.json")
with open("result.json", "w") as output:
    output.write(result.json())