Optimizing A Molecule’s Geometry with Psi4

Most often when we need to prepare a molecule for some perhaps docking, or some molecular dynamics simulation we would like to have the ligand in a relaxed state where the molecule would have the correct geometry. Most often I work with writing my own internal coordinates for molecules that serve as a good initial guess from SMILES and then I use Psi4, a quantum mechanical software that will optimize the geometry of the ligand.

So let’s get into it. First you can conda install psi4:

conda install -c psi4 psi4 

Next we need a ligand where we are going to optimize the geometry. Let’s look at methanol in Z-Matrix Form where the

• col 1 is the atom 1
• col2 is the bond connecting atom
• col3 is the bond length between two atoms
• col4 is the angle connecting atom
• col5 is the angle between the three atoms
• col6 is the dihedral connecting atom
• col7 is the dihedral angle between the four atoms

H11
O11  H11  0.9600
C11  O11  1.4000  H11  108.0000
H12  C11  1.1000  O11  112.0000  H11    0.0000
H13  C11  1.1000  O11  112.0000  H11  120.0000
H14  C11  1.1000  O11  111.8699  H11 -120.0000

Now if we take look at the geometry and output to some xyz coordinate:

import psi4

psi4.set_memory('1000mb')
psi4.core.set_num_threads(1)

zmatrix = '''\
H11
O11  H11  0.9600
C11  O11  1.4000  H11  108.0000
H12  C11  1.1000  O11  112.0000  H11  -60.0000
H13  C11  1.1000  O11  112.0000  H11   60.0000
H14  C11  1.1000  O11  111.8699  H11 -180.0000
'''

universe = psi4.geometry(zmatrix)
universe.update_geometry()
universe.print_out()

universe.save_xyz_file('test.xyz', False)

We get this:

Now for small molecules, that’s pretty good and easy to write the Z-matrix. However, that geometry I have written it’s only a somewhat good guess.

What we can do is use the template as a good initial guess and allow the QM software to relax the geometry. So let’s set up some options first with psi4. We want to setup what options we would like to pass into psi4:

• scf type: density fitting (these are approximation methods in calculating the potential energy of a system).
• g_covergence: This is geometric convergence and there is a force placed on the molecule which restricts it’s movement. If we place it gau tight it has a certain value that limits the wiggle of the molecule
• freeze_core: Freeze core is the whether you take into account the inner shell electrons when calculating the energy or only the outer valence electrons and their respective orbital.

import psi4

psi4.core.set_num_threads(1)
psi4.set_memory('1000mb')

psi4.set_options({
     'scf_type': 'df',
     'g_convergence': 'gau_tight',
     'freeze_core': 'true',
})

So now, let’s setup our optimization code:

zmatrix = '''\
H11
O11  H11  0.9600
C11  O11  1.4000  H11  108.0000
H12  C11  1.1000  O11  112.0000  H11  -60.0000
H13  C11  1.1000  O11  112.0000  H11   60.0000
H14  C11  1.1000  O11  111.8699  H11 -180.0000
0 1
'''

universe = psi4.geometry(zmatrix)
universe.update_geometry()

psi4.optimize(
    'hf/6-31G*',
    molecule=universe
)

universe.update_geometry()
universe.print_in_input_format()

So what’s going on here is the Z-matrix is getting inputed into the psi4 package and a Hartree-Fock (HF) level of theory with a basis set of (6–31G*) is passed into the optimizer code. We choose hf because we don’t the extended virtual orbitals give if we ran the optimization that included electron correlation like mp2 . It’s excessive to run your optimization with something like “mp2/aug-cc-pvdz”. It’s okay to use something cheap like HF.

The last line is to print in the format we injected the coordinates with and output looks like this:

    H
    O             1    0.946285
    C             2    1.399638      1  109.452715
    H             3    1.087420      2  112.035807      1  -61.226978
    H             3    1.087420      2  112.035807      1   61.226978
    H             3    1.081049      2  107.168437      1 -180.000000

Unfortunately, we have to map it back out to the numbers again but here we have the real coordinates of what our proposed methanol molecule would look like.

Alright that’s it for now.