Argon A-to-Z

This tutorial demonstrates a-to-z how to optimise Lennard Jones parameters for liquid argon, and without going into details. For details see other tutorials and wider MDMC documentation.

# Imports used for this tutorial
import numpy as np
import os
from scipy.interpolate import interp2d
from MDMC.control import Control
from MDMC.MD import Atom, Dispersion, LennardJones, Simulation, Universe
Supported DL_POLY version 4.10
# Change the number of threads depending on the number of physical cores on your computer
# as it was tested for LAMMPS
os.environ["OMP_NUM_THREADS"] = "4"
# Build universe with density 0.0176 atoms per AA^-3
density = 0.0176
# This means cubic universe of side:
# 23.0668 A will contain 216 Ar atoms
# 26.911 A will contain 343 Ar atoms
# 30.7553 A will contain 512 Ar atoms
# 38.4441 A will contain 1000 Ar atoms
universe = Universe(dimensions=23.0668)
Ar = Atom('Ar', charge=0.)
# Calculating number of Ar atoms needed to obtain density
n_ar_atoms = int(density * np.product(universe.dimensions))
print(f'Number of argon atoms = {n_ar_atoms}')
universe.fill(Ar, num_struc_units=(n_ar_atoms))
Universe created with:
  Dimensions       [23.07, 23.07, 23.07]
  Force field                       None
  Number of atoms                      0

Number of argon atoms = 216

In the Jupyter cell above, a box of Argon atoms is set up. However, at this point there is no interaction forces between the argon atoms! In the cell below an appropriate (for argon) force-field interaction potential is defined.

Ar_dispersion = Dispersion(universe,
                           (Ar.atom_type, Ar.atom_type),
                           function=LennardJones(epsilon=1.0243, sigma=3.36))

In this case the interaction potential chosen is the humble Lennard Jones (to get info see doc or type help(LennardJones)).

Also, a cutoff value is chosen (see help(Dispersion) for more info). A rule of thumb for Lennard-Jones is to pick cutoff=2.5*sigma. The value for argon is recommended to be between 8 and 12 ang. cutoff is not a force-field parameter and therefore will not be refined. Ideally, and for any system you want to pick at value of the cutoff which is small while not compromising accuracy. For this system picking a value between 8 and 12 ang is found to give near identifical results.

Next (and before starting the refinement), we set up the MD engine and equilibrate the system. Note with MDMC the equilibration only needs to be done once.

# MD Engine setup
simulation = Simulation(universe,
LAMMPS (29 Sep 2021 - Update 3)
  using 4 OpenMP thread(s) per MPI task
LAMMPS output is captured by PyLammps wrapper
LAMMPS (29 Sep 2021 - Update 3)
  using 4 OpenMP thread(s) per MPI task
LAMMPS output is captured by PyLammps wrapper
Total wall time: 0:00:00
Simulation created with lammps engine and settings:
  temperature  120.0

# Energy Minimization and equilibration
simulation.minimize(n_steps=5000), equilibration=True)

OK; time to set up the actual refinement of the force-field parameters.

First we need some data to refine against:

# exp_datasets is a list of dictionaries with one dictionary per experimental
# dataset
# Dataset from: van Well et al. (1985). Physical Review A, 31(5), 3391-3414
# resolution is None as the original author already accounted for instrument resolution
exp_datasets = [{'file_name':'data/Well_s_q_omega_Ar_data.xml',

The number of MD_steps specified must be large enough to allow for successful calculation of all observables. This depends the type of the dataset provided and the value of the traj_step (specified when creating the Simulation). If a value for MD_steps is not provided, then the minimum number needed will be used automatically.

Additionally, some observables will have an upper limit on the number of MD_steps that can be used in calculating their dependent variable(s). In these cases, the number of MD_steps is rounded down to a multiple of this upper limit so that we only run steps that will be useful. For example, if we use 1000 MD_steps in calculation, but a value of 2500 is provided, then we will run 2000 steps and use this to calculate the variable twice, without wasting time performing an additional 500 steps.

fit_parameters = universe.parameters

control = Control(simulation=simulation,
Control created with:
  Minimizer                             MMC
  FoM type               ChiSquaredExpError
  Number of observables                   1
  Number of parameters                    2

And finally start the refinement! Bump up n_steps from 3 when you are ready.

# Run the refinement, i.e. refine the FF parameters against the data
Step         FoM Change state epsilon (#2)   sigma (#3)
   0       505.8     Accepted        1.024         3.36
   1       507.3     Accepted        1.025        3.394
   2       511.8     Accepted        1.026        3.395
   3       533.3     Rejected        1.032        3.412

Final Parameters
 epsilon (#2)  sigma (#3)
     1.025789    3.394878

Automatic Scale Factors
  data/Well_s_q_omega_Ar_data.xml  0.225549
[ ]: