You can run this notebook online in a Binder session or view it on Github.

QCEngine

Full QCEngine documentation is available.

QCEngine is a quantum chemistry abstraction layer where many different quantum chemistry (or quantum-chemistry-like!) programs can be run with identical input and output abstractions that match the MolSSI QCSchema.

Begin by importing qcengine.

import qcelemental as qcel
import qcengine as qcng

We can list all programs that QCEngine currently supports. It should be noted that there are many programs which provide force field or machine learning potential evaluation (e.g. rdkit and torchani) in addition to the traditional quantum chemistry programs.

qcng.list_all_programs()

{'dftd3',
 'entos',
 'molpro',
 'mopac',
 'mp2d',
 'psi4',
 'rdkit',
 'terachem',
 'torchani'}

We can then list all programs that QCEngine has detected on the current resource. This list will vary depending on installed packages. As a note, QCEngine does not install programs by default, and these must be installed separately.

qcng.list_available_programs()

{'dftd3', 'mopac', 'psi4', 'rdkit'}

Single Computations

QCEngine makes the distinction between a “single” evaluation which corresponds to a single molecular geometry and a “procedure” which involves multiple geometries or multiple molecules. “Single” evaluations include energy, gradient, Hessian, and property quantities. “Procedures” include geometry optimization and other complex multi-step procedures.

First, we can build a Molecule object using the QCElemental molecule builder:

mol = qcel.models.Molecule(geometry=[[0, 0, 0], [0, 1.5, 0], [0, 0, 1.5]],
                           symbols=["O", "H", "H"],
                           connectivity=[[0, 1, 1], [0, 2, 1]])
mol

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

<Molecule(name='H2O' formula='H2O' hash='b41d0c5')>

We can then provide minimal input for a quantum chemistry job which specifies the molecule, driver, and model that the computation should be run under:

computation = {
    "molecule": mol,
    "driver": "energy",
    "model": {"method": "B3LYP", "basis": "6-31g"}
}
ret = qcng.compute(computation, "psi4")

The result contains many attributes that hold relevant data. We can access the return_result which contains the desired value as determined by the driver input field. In this case, it is the B3LYP/6-31g energy (in Hartree):

ret.return_result

-76.2741297206346

QCEngine automatically parses additional data about the state of the computation and pulls several other component fields. Here we can see the energy breakdown as well as the basis information:

ret.properties.dict()

{'calcinfo_nbasis': 13,
 'calcinfo_nmo': 13,
 'calcinfo_nalpha': 5,
 'calcinfo_nbeta': 5,
 'calcinfo_natom': 3,
 'nuclear_repulsion_energy': 11.138071187457696,
 'return_energy': -76.2741297206346,
 'scf_one_electron_energy': -126.25159666378747,
 'scf_two_electron_energy': 46.556895182916136,
 'scf_xc_energy': -7.717499427220989,
 'scf_dipole_moment': [0.0, 0.0, 2.660795024634264],
 'scf_total_energy': -76.2741297206346,
 'scf_iterations': 6}

Finally, QCEngine records much of the state of the computation such as the hardware it was run on, the program it was run with, and the versions of programs used:

ret.provenance.dict()

{'creator': 'Psi4',
 'version': '1.3.2',
 'routine': 'psi4.json.run_json',
 'memory': 2.266,
 'nthreads': 2,
 'qcengine_version': 'v0.9.0',
 'wall_time': 2.847166061401367,
 'hostname': 'Daniels-MacBook-Pro.local',
 'cpu': 'Intel(R) Core(TM) i7-7820HQ CPU @ 2.90GHz',
 'username': 'daniel'}

Procedures

Since we created a pretty poor molecule to start with, we should optimize it first under a force field method to have a reasonable geometry. Here, we will use the standalone geomeTRIC package coupled with the rdkit force field evaluator.

oh_bond, hoh_angle = mol.measure([[0, 1], [1, 0, 2]])
print(f"O-H Bond length (Bohr): {oh_bond}")
print(f"H-O-H Angle (degrees):  {hoh_angle}")

O-H Bond length (Bohr): 1.5
H-O-H Angle (degrees):  90.0

opt_input = {
    "keywords": {
        "program": "rdkit"
    },
    "input_specification": {
        "driver": "gradient",
        "model": {"method": "UFF"},
    },
    "initial_molecule": mol
}
opt = qcng.compute_procedure(opt_input, "geometric")
opt

<Optimization(model='{'method': 'UFF', 'basis': None}' molecule_hash='b41d0c5')>

We can first check the geometry of the final molecule and see that it is something much more reasonable:

opt_mol = opt.final_molecule
oh_bond, hoh_angle = opt_mol.measure([[0, 1], [1, 0, 2]])
print(f"O-H Bond length (Bohr): {oh_bond}")
print(f"H-O-H Angle (degrees):  {hoh_angle}")

O-H Bond length (Bohr): 1.8713297962085038
H-O-H Angle (degrees):  104.5102429904286

We can explore additional data generated with this geometry optimization including details of every gradient evaluation performed:

opt.trajectory

[<Result(driver='gradient' model='{'method': 'UFF', 'basis': None}' molecule_hash='b41d0c5')>,
 <Result(driver='gradient' model='{'method': 'UFF', 'basis': None}' molecule_hash='1ad5fe3')>,
 <Result(driver='gradient' model='{'method': 'UFF', 'basis': None}' molecule_hash='04ec4cf')>,
 <Result(driver='gradient' model='{'method': 'UFF', 'basis': None}' molecule_hash='58054eb')>,
 <Result(driver='gradient' model='{'method': 'UFF', 'basis': None}' molecule_hash='f2a154b')>,
 <Result(driver='gradient' model='{'method': 'UFF', 'basis': None}' molecule_hash='3c65f4c')>]

If desired, we can also look at the standard output of the geomeTRIC program:

print(opt.stdout)

9 internal coordinates being used (instead of 9 Cartesians)
Internal coordinate system (atoms numbered from 1):
Distance 1-2
Distance 1-3
Angle 2-1-3
Translation-X 1-3
Translation-Y 1-3
Translation-Z 1-3
Rotation-A 1-3
Rotation-B 1-3
Rotation-C 1-3
<class 'geometric.internal.Distance'> : 2
<class 'geometric.internal.Angle'> : 1
<class 'geometric.internal.TranslationX'> : 1
<class 'geometric.internal.TranslationY'> : 1
<class 'geometric.internal.TranslationZ'> : 1
<class 'geometric.internal.RotationA'> : 1
<class 'geometric.internal.RotationB'> : 1
<class 'geometric.internal.RotationC'> : 1
Step    0 : Gradient = 1.686e-01/1.827e-01 (rms/max) Energy =  0.0180063792
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 1.60000e-01 1.14610e+00 1.14610e+00
Step    1 : Displace = 1.073e-01/1.303e-01 (rms/max) Trust = 1.000e-01 (=) Grad = 1.080e-01/1.234e-01 (rms/max) E (change) =  0.0065007416 (-1.151e-02) Quality = 0.413
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 1.95811e-01 4.90425e-01 1.14610e+00
Step    2 : Displace = 8.518e-02/1.017e-01 (rms/max) Trust = 1.000e-01 (=) Grad = 1.592e-02/2.161e-02 (rms/max) E (change) =  0.0001690891 (-6.332e-03) Quality = 0.307
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 2.02540e-01 4.27767e-01 1.14610e+00
Step    3 : Displace = 1.326e-02/1.817e-02 (rms/max) Trust = 1.000e-01 (=) Grad = 1.744e-03/2.373e-03 (rms/max) E (change) =  0.0000051046 (-1.640e-04) Quality = 0.317
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 1.68278e-01 4.49086e-01 1.14610e+00
Step    4 : Displace = 4.295e-03/4.443e-03 (rms/max) Trust = 1.000e-01 (=) Grad = 1.311e-04/1.413e-04 (rms/max) E (change) =  0.0000000164 (-5.088e-06) Quality = 0.287
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 1.63098e-01 4.53449e-01 1.14610e+00
Step    5 : Displace = 2.091e-04/2.541e-04 (rms/max) Trust = 1.000e-01 (=) Grad = 9.256e-06/1.015e-05 (rms/max) E (change) =  0.0000000001 (-1.632e-08) Quality = 0.294
Converged! =D

Conclusion

These are some of the capabilities QCEngine offers, check out more documentation. If you like the project, consider starring us on GitHub or if you have any questions, join our Slack channel.