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 directory samples/minimisation/ you will find example inputs for the minimisation of the water molecule that you have already calculated a single point energy for. In the NWChem case, this is min_water.nwchem.py with equivalents 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.

In addition to the usual ChemShell output, DL-FIND will create further files that can be used to help monitor the course of the optimisation. In particular, the file path.xyz contains a list of XYZ geometries for each cycle that can be read into a visualiser to see how the optimisation has progressed.

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 NO2. First we read in the starting structure as before:

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

Then define the level of theory. NO2 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)

Finally the minimisation is requested in the same way as for water. When it has finished, visualise the optimisation in path.xyz and note the final conformation the molecule adopts.

Now take a look at min_no2p.nwchem.py. Positively charged NO2 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 before. When it is finished visualise the optimisation trajectory in path.xyz. How does the final structure compare with the uncharged case?