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
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,
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
# 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?