Optimizing molecular structure

Optimizing molecular structure#

The initial structure can be provided manually, or as constructed using SMILES strings (see previous section). Here we will do the latter to illustrate a full structure optimization workflow, with methanol as the example.

To minimize computational cost, we perform the final quantum chemical structure optimization at a minimal level of theory (HF with the STO-3G basis set). For practical calculations, DFT or MP2 as well as a larger basis set should be used.

import numpy as np
import py3Dmol as p3d
import veloxchem as vlx

Initial structure#

We obtain the initial (force field optimized) structure from Open Babel, using the SMILES string (CO).

methanol_ff = vlx.Molecule.read_smiles("CO")

xTB optimization#

Next, set up the xTB driver and perform the structure optimization.

xtb_drv = vlx.XtbDriver()

xtb_opt_drv = vlx.OptimizationDriver(xtb_drv)
xtb_opt_drv.ostream.mute()

xtb_opt_results = xtb_opt_drv.compute(methanol_ff)

methanol_xtb = vlx.Molecule.read_xyz_string(xtb_opt_results["final_geometry"])

HF optimization#

Set up the SCF and optimization drivers and performing the final optimization, using the xTB results as the initial structure.

Comparison#

Visualize the results using py3Dmol:

viewer = p3d.view(viewergrid=(1, 3), width=500, height=200, linked=False)

viewer.addModel(methanol_ff.get_xyz_string(), "xyz", viewer=(0, 0))
viewer.addModel(methanol_xtb.get_xyz_string(), "xyz", viewer=(0, 1))
viewer.addModel(methanol_hf["final_geometry"], "xyz", viewer=(0, 2))

viewer.setViewStyle({"style": "outline", "width": 0.05})
viewer.setStyle({"stick": {}, "sphere": {"scale": 0.25}})

viewer.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

These structures are all very similar, featuring, e.g, the staggered H-O-C-H dihedral. To better see any differences in the structures, the distance matrices can be calculated using functionalities from the molecule object:

#print("From force field (MMFF94):")
#dm_ff = methanol_ff.get_distance_matrix_in_angstrom()
#print(np.around(dm_ff, 3))
#print()

#print("From xTB:")
#dm_xtb = methanol_xtb.get_distance_matrix_in_angstrom()
#print(np.around(dm_xtb, 3))
#print()

#print("From HF:")
#dm_hf = methanol_hf.get_distance_matrix_in_angstrom()
#print(np.around(dm_hf, 3))

This gives the the distances between all pairs of atoms, as expressed in Å.

As long as the atom order is consistent (which it here is), we can print the differences in distance matrices, or maximal absolute difference:

#print("Difference between HF and xTB:")
#print(np.around(dm_hf - dm_xtb, 3))
#print()
#print("Maximal absolute difference between HF and xTB:")
#print(f"{np.max(np.abs(dm_hf - dm_xtb)):.4f} Å")

The largest difference is thus seen to be 0.04 Å, between atoms 2 and 5. These atoms are the oxygen and one hydrogen from CH\(_4\), which has a total distance of 2.05 Å. As any pair-wise differences will propagate along the molecule, it may be more relevant to look at relative differences:

# replace the diagonal elements (zeros) with a large number to avoid nan
#dm_tmp = dm_hf + 1e10 * np.diag(np.ones(len(dm_hf)))

#print("Relative difference between HF and xTB:")
#print(np.around((dm_hf - dm_xtb) / dm_tmp, 3))
#print()
#print("Maximal relative difference between HF and xTB:")
#print(f"{np.max(np.abs((dm_hf - dm_xtb)/dm_tmp)):.4f}")

Now we see that the largest relative difference in atom-pair distance is between atoms 2 and 6, which correspond to the oxygen and the hydrogen directly connected to the oxygen. With this, we have thus found the largest relative deviation in a bond length, which is in most cases more important than in distances between atoms at different sites.

Contents