Day 8 Force-fields, potentials and optimisation methods¶

Introduction¶

As you heard in the lectures, a vital component of most simulations is the specification of a force-field or potential. In this exercise you will investigate the behaviour of two potentials for alumina, a rigid ion pair potential model with fitted ionic charges and a shell model pair potential with full ionic charges. Neither model is completely satisfactory as a model for alumina. By the end of the exercise, you should be able to answer two questions.

1. What kinds of investigation would you use each model for?

2. What further tests would you think it necessary to do for a specific problem?

On the way, we shall also investigate the use of a number of optimisation methods. We shall use the GULP code developed by Julian Gale and co-workers. This is a general-purpose code that, amongst other things, performs optimisations on lattice structures, surfaces and point defects. Full details of the program’s capabilities can be found in the GULP Manual.

Setup (optional)¶

This step is needed only if you run on google colab. Instructions may work on other systems but are not tested.

[ ]:

! pip install data_tutorials ase weas_widget
! apt install gfortran

Now bring the data needed for the exercises.

[ ]:

from data_tutorials.data import get_data
get_data(
    url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_8/optimisation/",
    filename=["bixbyite1.gin","bixbyite2.gin","corundum1.gin","corundum2.gin","rhodium1.gin","rhodium2.gin"],
    folder=".",
)

get_data(
    url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/.raw/",
    filename=["gulp"],
    folder="/usr/local/bin",
)
! chmod 755 /usr/local/bin/gulp

Some cells were added to allow you execute the commands needed for the tutorial. feel free to add more as needed or use the console to carry on the tasks

The pair potential \(V_{ij}(r)\) takes the form

where the first term is the electrostatics term and the short-range Buckingham term has three fitted parameters \(A_{ij}\), \(\rho_{ij}\) and \(C_{ij}\). The rigid ion model used here has a further parameter \(Z_i\). In the shell model used, this is fixed at +3 (Al). The shell model has two further parameters that describe the polarisability of the ions. The model comprises a massless shell (supposed to represent the outer electrons) of charge \(Y_i\) linked to a massive core (representing the nucleus and inner electrons) with charge \(Z_i - Y_i\) by a spring constant \(k_i\). The polarisability of the free ion in vacuum is then given by \(\alpha_i=\frac{Y_i^2e^2}{4\pi\varepsilon_0 k_i}\). All these parameters can be fixed either by using ab initio calculations as a training set or by fitting to experimental data. Details of these force-fields can be found in ref [Bush1994] (shell model) and ref [Gale1992] (rigid ion model).

The input file¶

This is a complex input file since we shall perform a large number of calculations by commenting and un-commenting various commands. Remember that in GULP input files, anything that follows a # on a line is commented out. You are given six input files in all: three structures – corundum (trigonal), bixbyite (cubic) and rhodium oxide (orthorhombic) and two potentials for each structure. An example input file is shown below with added comments for each line.

opti conp full # keywords to optimise structures at constant pressure
#keyword force_minimisation # Minimise to zero force not minimum energy
#keyword conj          # use conjugate gradients minimisation
#keyword dfp           # use Davidon-Fletcher-Powell variable metric method
#keyword rfo           # use rational functional optimisation
#switch bfgs cycle 10  # switch to bgfs method after 10 cycles
#keyword nosym         # Perform minimisation using no symmetry elements

#keyword linmin    # print out details of line minimisations
#terse potentials  # suppress details of potential
#terse derivatives # suppress details of derivatives

maxcyc 100         # maximum number of minimisation cycles
stepmx opt 1.0   # maximum size of minimisation step
update 10        # complete recalculation of Hessian after n cyclea

#keyword phonon      # calculate phonon spectrum using details below
#keyword free_energy # Minimise on the free energy
#keyword zsisa       # Approximate minimisation on free energy

keyword property  # Calculate bulk lattice properties
temperature 298   # Set temperature (K) (for free energy minimisation)
pressure    0     # Set pressure (GPa) and minimise on the enthalpy

title  # Lines after this keyword give the title
alumina corundum structure file
Potential is Bush (1994)
end    # End of title section
cell   # Define unit cell of lattice on the next line
4.7602   4.7602  12.9933  90.000000  90.000000 120.0
frac   # Define the species in the unit cell
Al core 0.000000   0.000000   0.352160
Al shel 0.000000   0.000000   0.352160
O  core 0.306240   0.000000   0.250000
O  shel 0.306240   0.000000   0.250000
Space  # Space group of lattice (standard ordering)
167

supercell 1 1 1  # Define a supercell of size l by n by m unit cells
#shrink 1 1 1    # Construct grid to calculate lattice vibration entropy
#dispersion 3 20 # Calc dispersion curve (3 directions, 20 points each)
#0.0 0.0 0.0 to 0.5 0.0 0.0
#0.0 0.0 0.0 to 0.5 0.5 0.0
#0.0 0.0 0.0 to 0.5 0.5 0.5

species # Define the species in the unit cell (label, charge)
Al core  0.043
Al shel  2.957
O  core  0.513
O  shel -2.513
buckingham # Pair potential type and cutoff (see lectures for details)
Al shel O shel  2409.505 0.2649  0.00 0.0 10.0
O  shel O shel    25.410 0.6937 32.32 0.0 12.0
spring # spring constants for shell model
Al 403.98
O   20.53

#cutp 12.0 mdf 2.0  # Taper function for potential cutoff

output xyz corundum1    # Outputs coordinate file for viewing
#output phonon corundum1 # Outputs dispersion data and dos for plotting

The exercises¶

The Gulp code can be run using the simple command: gulp <infile >outfile at the command line or incorporating it into a suitable jobscript. I strongly suggest that you retain the original input files and copy them each time you want to use them. All input files are in ~/WORKSHOP/Day_5/optimisation

1. Run the six input files (corundum1.gin, corundum2.gin, bixbyite1.gin, bixbyite2.gin, rhodium1.gin, rhodium2.gin) and compare the energies per formula unit for the three structures using the two potentials. Ab initio results [Sarsam2013] predict that the corundum phase is the most stable, the bixbyite phase is 0.2eV/formula unit above that and the \(Rh_2O_3\) orthorhombic phase even higher in energy. The experimental hexagonal lattice parameters for corundum are \(a_0=4.7602A\), \(c_0=12.9933\) and the fractional z coordinates are 0.3522 (Al) and 0.3062 (O). Also compare the calculated elastic and dielectric constants with the experimental values given below (in GPa for elastic constants)

[1]:

%%bash

gulp < corundum1.gin > corundum1.out
gulp < corundum2.gin > corundum2.out
gulp < bixbyite1.gin > bixbyite1.out
gulp < bixbyite2.gin > bixbyite2.out
gulp < rhodium1.gin > rhodium1.out
gulp < rhodium2.gin > rhodium2.out

you can visualise the structures from each of the calculations.

[ ]:

from ase.io import read
from weas_widget import WeasWidget

c1 = read("corundum1.cif")
c2 = read("corundum2.cif")

v=WeasWidget()
v.from_ase([c1,c2])
v.avr.model_style = 1
v

inspect outputs, check each as created above

What are your conclusions about the usefulness of these models in modelling ground-state properties? How about point defects, say a vacancy or doping the lattice with a rare earth ion? 2. There are a number of optimisation schemes available and all of them can be fine-tuned. You will have noticed that the different structures required different numbers of cycles to minimise under the bfgs scheme. Why do you think this is? Consider one of the files and investigate the effects of changing

• the maximum step length of the minimisation (uncomment the stepmx parameter and change the number)

• the optimisation scheme (uncomment dfp, conj, rfo in turn). You will get the most variable behaviour with rhodium1.gin and rhodium2.gin.

[ ]:

%%bash

It can be useful to change the optimisation scheme part-way through the calculation. Try this using the switch command.

1. The various structures cannot transform into each other because they are constrained by symmetry. You can investigate this by uncommenting the keyword nosym instruction and also ensuring that the repeating units are all the same size by uncommenting the supercell instruction and altering the integers to make the number of formula units in the simulation cell the same in all the cases. What is the effect of this? If the structures cannot be made to transform into each other by this kind of activity that suggests that each structure occupies a local minimum. Stronger evidence for this can be found in the next exercise.

[ ]:

%%bash

2. The calculation of the phonon spectrum (lattice dynamical behaviour) is often a sensitive test of the performance of a potential. Take the corundum1.gin and corundum2.gin files and uncomment the keyword phonon line and the dispersion line together with the three following lines and calculate the phonon spectrum for the two potentials. Can you see any negative values for the frequencies (apart from the translational modes)? What does this tell you about the stability of the minima in which the structures sit? Are there any significant differences between the rigid ion and the shell model calculations?

[ ]:

%%bash

3. So far we have looked at the behaviour of the forcefield at standard temperature and pressure. The next two exercises will look at the temperature-dependent behaviour (using the quasi-harmonic approximation). We first consider the temperature behaviour. This is done using the keyword free_energy command which optimises the free energy of the lattice not the configurational energy. To do this we require a reasonable representation of the phonon density of states to calculate the vibrational entropy. We use the shrink command to do this by constructing a grid of reciprocal lattice points. First ensure that the free energy is converged with respect to the grid by increasing the values of the three integers until the optimised free energy no longer changes significantly. Then perform calculations for a range of temperatures and plot the free energy and the lattice volume against temperature. The experimental data is [Fiquet1999]

• Does it make a difference if you use the rigid ion potential in corundum2.gin? What happens with the other lattices? Does using the free energy rather than the configurational energy (as you did before) have any effect on the relative stability of the crystal structures predicted by the forcefield.

[ ]:

%%bash

1. Investigate the effect of pressure on the structures using the pressure command and a range of pressure values (given in GPa). This now causes the program to optimise with respect to the enthalpy. Is there any pressure where the corundum structure becomes more stable than bixbyite?

2. In the light of all this, what kind of simulation would you be prepared to use either of the forcefields for?

3. Do you have any thoughts on why both forcefields fail to predict that corundum is the stable structure for alumina at room temperature and pressure?

[ ]:

%%bash

[ ]: