# 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 scipy.interpolate import interp2d
from MDMC.control import Control
from MDMC.MD import Atom, Dispersion, LennardJones, Simulation, Universe
```

```
[2]:
```

```
# 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"
```

```
[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.)
# 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 2)
using 4 OpenMP thread(s) per MPI task
LAMMPS output is captured by PyLammps wrapper
LAMMPS (29 Sep 2021 - Update 2)
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
```

```
[6]:
```

```
# Energy Minimization and equilibration
simulation.minimize(n_steps=5000)
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',
'reader':'xml_SQw',
'weight':1.,
'auto_scale':True,
'resolution':None}]
```

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
control = Control(simulation=simulation,
exp_datasets=exp_datasets,
fit_parameters=fit_parameters,
MC_norm=30,
minimizer_type="MMC",
reset_config=False,
MD_steps=570)
```

```
Control created with:
Minimizer MMC
MC norm 30
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=3)
```

```
Step FoM Change state epsilon sigma
0 420.5 Accepted 1.024 3.36
1 427.6 Rejected 1.016 3.36
2 432.3 Accepted 1.025 3.375
3 406.2 Accepted 1.018 3.383
Final Parameters
epsilon sigma
1.017995 3.382975
Automatic Scale Factors
data/Well_s_q_omega_Ar_data.xml 0.219311
```

```
[ ]:
```

```
```