How to fit a potential using GULP and LAMMPS

How to use GULP

1. Download and make GULP (instructions are included in a README with GULP). Include any packages that you require for the specific potentials you are using.

2. GULP can be run by executing the gulp program and giving it an input file, you will also want to pipe this out to an output file:

path/gulp-x.x/Src/gulp Gulp_Input > Gulp_Output

GULP comes with many different example input files, you will want to run GULP from the folder with the input file.

3. GULP structural data is given in the input file in the following format:

# Start of Structure
vectors
3.1023 0.0000 0.0000
1.4990 2.7214 0.0000
1.5040 0.8336 5.0333
fractional
Fe core 0.0000 0.0000 0.0000
Fe core 0.5000 0.5000 0.5000
O core 0.7500 0.7500 0.7500
O core 0.2500 0.2500 0.2500
space
P 1
observable
energy ev
-15.5000 100.0
end

Fractional can be changed for cartesian and both cores and shells can be defined. Atomic charges can be added to the end of each atom definition:

Fe core 0.0000 0.0000 0.0000 4.0

and many different observables can be defined, in this case the binding energy is given with a weighting.

4. At the top of the GULP input file you need to define what you want GULP to do. For a potential fit:

fit abs conv
maxcyc 5000
title
funny/creative title goes here

In this example, "fit" means we are fitting something (does not need to be a potential) and "maxcyc" defines the maximum number of fitting cycles are done.

5. At the bottom of the input file we define the potentials with their initial parameters and which of these are being fitted:

buck
Fe core O core 1388.0000 2.7600 175.0000 0.00 6.00 1 1 0

This is an example of a Buckingham potential between Fe core and O core. The first three numbers are the three parameters of the Buckingham potential, the fourth and fifth are the minimum and maximum radius over which this potential acts and the last three numbers define which of the parameters are allowed to vary, 1 means it is varied, 0 means it is fixed.

6. Putting this all together we can get an input file:

fit abs conv
maxcyc 5000
title
funny/creative title goes here

# Start of Structure
vectors
3.1023 0.0000 0.0000
1.4990 2.7214 0.0000
1.5040 0.8336 5.0333
fractional
Fe core 0.0000 0.0000 0.0000
Fe core 0.5000 0.5000 0.5000
O core 0.7500 0.7500 0.7500
O core 0.2500 0.2500 0.2500
space
P 1
observable
energy ev
-15.5000 100.0
end

buck
Fe core O core 1388.0000 2.7600 175.0000 0.00 6.00 1 1 0

Which will attempt to fit a Buckingham potential for Fe core and O core to the FeO structure and energy provided. This WILL NOT change the structure itself, only the parameters in the potential that are allowed to vary. When doing these fits, you do not want to allow all parameters to change as we have not provided any forces, only atomic positions. Start with a potential that is already close to what you expect the final potential to be and allow 1 or 2 of the parameters to vary. Keep cycling with parameters can vary between fitting runs, this may need to be done many times depending on how many parameters you have to fit.

Converting potentials from GULP to LAMMPS

Not all potentials that are found in GULP are found in LAMMPS and vice versa, be careful when selecting which potential to use. The same potential are also not guaranteed to have the same parameters in both GULP and LAMMPS, for example the 2-body and 3-body Stillinger-Weber potential in GULP is given by:

$E_{ij}=A_{ij}\cdot \exp \left({\frac {\rho _{ij}}{r_{ij}-r_{max}}}\right)\cdot \left({\frac {B_{ij}}{r_{ij}^{4}}}-1\right)$

$E_{ijk}=K_{ijk}\cdot \exp \left({\frac {\rho _{ij}}{r_{ij}-r_{max}}}\right)\cdot \exp \left({\frac {\rho _{ik}}{r_{ik}-r_{max}}}\right)(\cos(\theta _{jik})-\cos(\theta _{0}))^{2}$

Whereas in LAMMPS this same potential is given by:

$E_{ij}=A_{ij}\epsilon _{ij}\cdot \exp \left({\frac {\sigma _{ij}}{r_{ij}-a_{ij}\sigma _{ij}}}\right)\cdot \left(B_{ij}({\frac {\sigma _{ij}}{r_{ij}}})^{p_{ij}}-({\frac {\sigma _{ij}}{r_{ij}}})^{q_{ij}}\right)$

$E_{ijk}=\lambda _{ijk}\epsilon _{ijk}\cdot \exp \left({\frac {\gamma _{ij}\sigma _{ij}}{r_{ij}-a_{ij}\sigma _{ij}}}\right)\cdot \exp \left({\frac {\gamma _{ik}\sigma _{ik}}{r_{ik}-a_{ik}\sigma _{ik}}}\right)(\cos(\theta _{jik})-\cos(\theta _{0}))^{2}$

In GULP, parameters that appear in both the two-body and three-body potentials don't necessarily need to be the same, whereas in LAMMPS they are due to the way the Stillinger-Weber potential is written in an input for LAMMPS. Within the input file for LAMMPS, a Buckingham potential given above for use in GULP is written as:

pair_style buck 6.00
pair_coeff 1 2 1388.00 2.76 175.00

where the pair_style line defines the type of potential being used and its cutoff distance. In the pair_coeff line the two atoms this acts between is given (if 1 is defined as Fe and 2 is O) along with the three parameters for this potenital.