Fitting a Lennard-Jones potential for bulk Cu

In this tutorial you will learn how to create a dataset of structures, calculate their energies, and use the FitLJ package to fit a Lennard-Jones potential.

Creating the database

A necessary requirements for fitting an interatomic potential is to have a database of structures representative of the possible atomic configurations of the material of interest. The energies of these configurations are used as ground truth for the fitting procedure.

Here, we will use ASE to create a data set of different bulk structures of Cu at different volumes. To make the potenital as general and transferable as possible among different configuratiosn, we will consider three different crystalline systems:

• Simple cubic:

  

• Face-centered cubic:

  

• Body-centerd cubic:

  

The volume of the three crystalline structures is defined by the lattice constant a. The ground truth energies will be calculated using the Effective Medium Theory.

Warning

The Effective Medium Theory used in this tutorial is not intended to produce accurate results and it is intended only for pedagogical purposes. For production simulations, one might want to use a more accurate and computationally intensive electronic structure method such as Density Functional Theory.

The script run_emt uses ASE to loop over different values of the lattice constant a for the three different crystalline structures. The calculated structures and associated energies are then dumped into a database. The latter will be used later for the fitting procedure.

import numpy
from ase import Atoms 
from ase.calculators.emt import EMT 
from ase.lattice.cubic import SimpleCubic, FaceCenteredCubic, BodyCenteredCubic
from ase.db import connect

# Define list of lattice constant values
lattice_constants = numpy.linspace(2.0,6.0,10)

# Connect to the database
database = connect('Cu.db')

# Loop over the crystalline lattice systems 
for lattice in [SimpleCubic, FaceCenteredCubic, BodyCenteredCubic]:

    # Loop over the lattice constant values 
    for a in lattice_constants:

        # Create the crystal unit cell    
        atoms = lattice('Cu',latticeconstant=a)

        # Attach a calculator 
        atoms.calc = EMT()

        # Get the potential energy 
        energy = atoms.get_potential_energy()

        # Write the structure to the database
        database.write(atoms)

Fitting the Lennard-Jones potential

Once the database has been created, the FitLJ program can be used to read the database and use it to fit the Lennard Jones potential parameters \(\sigma\) and \(\epsilon\).

The script fit load the database Cu.db previously created and fit the parameters of a Lennard Jones potenital to the energies of the structures contained in the database.

from fitlj.fit_model import FitModel
from fitlj.io import AtomicSimulationEnvironmentDataBase

# Load the database file 
database = AtomicSimulationEnvironmentDataBase(database_filename='Cu.db')

# Connect to the database 
database = database.connect_to_database()

# Initial guess for sigma
sigma_initial_guess = 15.0

# Initial guess for epsilon 
epsilon_initial_guess = 3.5 

# Fit a Lennard Jones model to the configurations present in the database
model = FitModel(database=database, 
                 lennard_jones_parameters=[sigma_initial_guess, epsilon_initial_guess])
sigma, epsilon = model.fit_lennard_jones_model()

# Print the fitted values
print(f'The fitted value of sigma is: {sigma:.5f} eV.')
print(f'The fitted value of epsilon is: {epsilon:.5f} Ang.')

The optimized values of \(\sigma\) and \(\epsilon\) are printed after the optimization has successfully terminated:

>>>python3 fit.py
Optimization terminated successfully.
         Current function value: 730.106876
         Iterations: 22
         Function evaluations: 72
The fitted value of sigma is: 19.05469 eV.
The fitted value of epsilon is: -0.50313 Ang.

Warning

In order for the optimization procedure to be succesfull, reasonable starting guesses for the \(\sigma\) and \(\epsilon\) parameters must be chosen. If the starting guesses are very far from the final, optimized values, the fitting procedure might not be successful.