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

[1]:

# Imports used for this tutorial
import numpy as np
import os
from MDMC.control import Control
from MDMC.MD import Atom, Dispersion, LennardJones, Simulation, Universe

Supported DL_POLY version 5.0

[2]:

# Change the number of threads depending on the number of physical cores on your computer
# as it was tested for LAMMPS

[3]:

# 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., mass=36.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.

[4]:

Ar_dispersion = Dispersion(universe,
(Ar.atom_type, Ar.atom_type),
cutoff=8.,
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.

[5]:

# MD Engine setup
simulation = Simulation(universe,
engine="lammps",
time_step=10.18893,
temperature=120.,
traj_step=15)

LAMMPS (29 Sep 2021 - Update 3)
LAMMPS output is captured by PyLammps wrapper
LAMMPS (29 Sep 2021 - Update 3)
LAMMPS output is captured by PyLammps wrapper
Total wall time: 0:00:00
Simulation created with lammps engine and settings:
temperature  120.0


[6]:

# Energy Minimization and equilibration
simulation.minimize(n_steps=50)
simulation.run(n_steps=10000, equilibration=True)


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

First we need some data to refine against:

[7]:

# 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',
'type':'SQw',
'weight':1.,
'auto_scale':True,
'resolution':800}]


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.

[8]:

fit_parameters = universe.parameters
fit_parameters['sigma'].constraints = [2.8,3.8]
fit_parameters['epsilon'].constraints = [0.6, 1.4]

control = Control(simulation=simulation,
exp_datasets=exp_datasets,
fit_parameters=fit_parameters,
minimizer_type="GPO",
reset_config=True,
MD_steps=4000,
equilibration_steps=4000)

Control created with:
- Attributes                              -
Minimizer                             GPO
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.

[9]:

# Run the refinement, i.e. refine the FF parameters against the data
control.refine(n_steps=25)
control.plot_results();

Step         FoM Change state  Pred coords     Pred FoM epsilon (#2)   sigma (#3)
0       19.69     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69        1.024         3.36
1       79.19     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69        1.035        3.532
2       322.5     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69        1.075        2.905
3       411.2     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69       0.9356        3.782
4       71.23     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69        1.396         3.47
5         310     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69       0.7189        2.876
6       64.88     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69       0.9139        3.198
7       207.8     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69       0.8076        3.693
8       65.55     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69       0.6141        3.271
9       731.6     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69        1.174        2.812
10       20.04     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69       0.9896        3.336
11       397.6     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69        1.218        2.966
12       204.4     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69        0.747        3.001
13       61.16     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69       0.7752        3.566
14       400.2     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69        1.141        3.713
15       59.48     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69        0.665        3.619
16       299.6     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69        1.354        3.098
17       60.12     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69        1.304        3.223
18       109.7     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69        1.119        3.104
19       22.72     Accepted [1.0243 kJ / mol, 3.36 Ang]        19.69       0.8776        3.441
20       16.78     Accepted [1.2497910435532642 kJ / mol, 3.3191472337787253 Ang]        16.78         1.25        3.319
21       24.17     Accepted [1.2497910435532642 kJ / mol, 3.3191472337787253 Ang]        16.78       0.6007         3.46
22       31.56     Accepted [1.2497910435532642 kJ / mol, 3.3191472337787253 Ang]        16.78        1.393         3.35
23       23.41     Accepted [1.2497910435532642 kJ / mol, 3.3191472337787253 Ang]        16.78       0.7535        3.386
24       22.85     Accepted [1.2497910435532642 kJ / mol, 3.3191472337787253 Ang]        16.78        1.143        3.292

The refinement has finished.

Minimum measured point is:
(1.2497910435532642 kJ / mol, 3.3191472337787253 Ang) with an FoM of 16.783812233483726.

Minimum point predicted is:
(1.2497910435532642 kJ / mol, 3.3191472337787253 Ang) for an FoM of 16.783812233483726.

Automatic Scale Factors
data/Well_s_q_omega_Ar_data.xml  20.210148
Parameter means = [1.13543659 3.34781567], Parameter errors = [0.16434543 0.03308508]

[ ]: