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