*********************************
Classical forcefield calculations
*********************************

Until now our examples have used quantum mechanical programs for the evaluation
of energies and gradients. From a user point of view first principle ab initio
calculations are simple to set up as they require only a geometry, charge,
spin multiplicity and basis set. 

In this tutorial we look at classical calculations which require a forcefield
to be defined, using the example of the `GULP <http://gulp.curtin.edu.au/gulp/>`_ package. The examples can be found
in ``samples/forcefields_gulp``.

Running GULP from ChemShell
===========================

As with the QM programs we have used in earlier tutorials, GULP is an external 
code and must be obtained separately from its developers. Once you have 
obtained a copy, you can run GULP either as a standalone program or as a 
"linked-in" ChemShell library - see the installation guide for details.


Shell model calculations for ionic materials
============================================

ChemShell recognises structures with periodic boundary conditions and can 
pass these to GULP for tasks such as lattice relaxation. 

In this example, ``mgo.py``, we relax the atomic positions of an MgO 
lattice. First, the unrelaxed structure is read in from disk::

 # Read in MgO crystal data from disk
 mgo = Fragment(coords="mgo.cjson", shellatoms=[0,1,2,3,4,5,6,7])

Note here that ChemShell supports multiple file formats for reading in 
atomic data, and here we are using the `Chemical JSON <http://wiki.openchemistry.org/Chemical_JSON>`_ 
format. This goes 
beyond XYZ to allow us to specify a unit cell and fractional coordinates 
for example (take a look at the ``mgo.cjson`` file to see how these 
are defined).

As we are performing a shell model calculation, shells must be added
to the MgO when it is read in. This is carried out by adding 
a ``shellatoms=`` option to the ``Fragment()`` line. Here, we are 
adding shells to both the Mg and O ions.

Next, the forcefield must be defined. This is written in standard GULP
format, stored in a triple-quoted (multiline) string that is passed 
to GULP in its input file::

 # Define GULP forcefield for MgO
 ff = '''species
 Mg core  1.580
 Mg shel  0.420
 O  core  0.513
 O  shel -2.513
 buckingham
 Mg shel O shel 2457.243 0.2610  0.00 0.0 10.0
 O  shel O shel   25.410 0.6937 32.32 0.0 12.0
 spring
 Mg 349.95
 O   20.53
 '''

There are `libraries of potentials <https://www.ucl.ac.uk/klmc/Potentials/>`_ 
available on the UCL Materials Chemistry web site.

Next, the forcefield calculations is set up in a similar way to QM 
calculations. An additional option ``ff=`` is used to specify the forcefield
definition::

 # Set up forcefield calculation
 my_theory = GULP(frag=mgo, ff=ff)

Finally, an optimisation is requested in the usual way::

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

Note that DL-FIND does not currently support optimisation of cell constants, 
so the energy is minimised under the constraint of a fixed unit cell.


Molecular mechanics for organic systems
=======================================

As well as techniques for materials modelling, GULP also supports molecular 
mechanics (MM) calculations typically used to model organic molecules.

In the example ``butanol.py`` we see how an MM calculation is set up 
for a butanol molecule. As with MgO, the butanol structure is imported 
from a Chemical JSON file, this time in order to specify partial charges as 
well as atomic positions::

 # Read in butanol structure from disk
 butanol = Fragment(coords="butanol.cjson")

Note that the molecular mechanics forcefield we will be using does not 
require shells so we do not add shells in this line. 

Secondly, the forcefield is defined in a standard GULP format as before. 
Importantly, the keyword ``molmec`` is included to specify that it is 
a molecular mechanics forcefield::

 # molmec removes Coulomb terms between bonded atoms
 # and between atoms bonded to a common third atom
 keyword molmec

After this comes a list of bond, angle, torsion and van der Waals terms
characteristic of a molecular mechanics forcefield.

When the GULP calculation is set up, be sure to set the option ``molecule=True``
for a molecular mechanics calculation::

 # Set up forcefield calculations
 my_theory = GULP(frag=butanol, ff=ff, molecule=True)

A geometry optimisation is then requested in the usual way. As this is 
a gas phase molecular system there are no periodic boundary conditions
to worry about.