Energy minimisation

ChemShell contains a wide range of geometry optimisation algorithms through its close integration with the DL-FIND geometry optimisation library. The simplest optimisation task that can be carried out is to minimise the energy of a system on its potential energy surface, which we demonstrate in this tutorial.

Minimisation of water

In the previous tutorial we calculated the single-point energy of a water molecule. We will now look at minimisation of the same molecule using the example inputs in the directory samples/minimisation/. For the NWChem case, the example is min_water.nwchem.py, with equivalents also available for GAMESS-UK and ORCA.

Take a look at the input. You will see that most of the commands are exactly the same as in the single-point case, with the same geometry being read in and the same theory definition. The only difference is the lines requesting the task to be carried out:

# Request an optimisation using DL-FIND
opt = Opt(theory=my_theory)
opt.run()

In this case, we have requested a DL-FIND geometry optimisation by defining an Opt() object. As we haven’t specified any options in this object besides the name of the theory object, it will run an energy minisation by default. As before, calling the run() method will execute the calculation.

When the calculation is run, DL-FIND will print information on the energy, step taken and gradient changes at each cycle, and will continue until convergence on all these measures is reached. When the calculation has finished, compare the energy output at the end with the energy reported in the first cycle to check that it has indeed been lowered.

After the optimisation has completed, the save() method can be used to save the optimised geometry to disk for further analysis:

# Save the optimised geometry to disk
my_mol.save("water_opt.xyz")

Load water.xyz and water_opt.xyz into your visualiser. Note that the starting geometry was close to the minimum in order to make sure that the example ran quickly, but you should still see a small difference in the atomic positions between the starting and optimised geometries.

Minimisation of nitrogen dioxide

In the second example, we look at how charged and open shell calculations can be performed and how these affect the resulting optimised structure for the case of nitrogen dioxide.

First let’s look at the file min_no2.nwchem.py which runs an optimisation on uncharged NO₂. First we read in the starting structure as before:

# Read in NO2 geometry from disk
my_mol = Fragment(coords="no2.xyz")

Then we define the level of theory. NO₂ is an open shell doublet system, which we select with the options scftype="uhf", which specifies an unrestricted HF or DFT calculation, and mult=2 to set the doublet spin multiplicity:

# Set up the QM calculation
my_theory = NWChem(frag=my_mol, method="dft", functional="blyp", basis="3-21g",
                  scftype="uhf", mult=2)

The minimisation is then requested in the same way as for water, except this time we will introduce further DL-FIND options:

# Request an optimisation using DL-FIND
opt = Opt(theory=my_theory, tolerance=0.0003, save_path=True)
opt.run()

Here, we show how to set a tighter threshold for convergence of the optimisation, using the tolerance option, and include the option save_path which instructs ChemShell/DL-FIND to save further input and output files that can be used to help monitor the course of the optimisation. These will be saved in a folder called _dl_find. Of particular interest is the file path.xyz, which contains a series of XYZ geometries corresponding to the full optimisation path. When the calculation is finished, open _dl_find/path.xyz in your visualiser and take a look at how the structure changed as the optimisation progressed. Note the final conformation that the molecule adopts, which will be the same as the geometry saved to no2_opt.xyz.

Now let’s take a look at min_no2p.nwchem.py. Positively charged NO₂ is a closed shell system, so we use the default restricted SCF calculation. The charge is specified by the option charge=1:

# Set up the QM calculation
my_theory = NWChem(frag=my_mol, method="dft", functional="blyp", basis="3-21g",
                   charge=1)

The optimisation is then run in the same way as in the uncharged system. When it is finished read in the two optimised geometries no2_opt.xyz and no2p_opt.xyz into your visualiser. How does the positively charged structure compare with the uncharged case?