.. |NO2| replace:: NO\ :sub:`2`\ 

*******************
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 |NO2|. 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. |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)

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