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 package. The examples can be found
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
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
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
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 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
# 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
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.