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
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., 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)
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
using multi-threaded neighbor list subroutines
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',
'reader':'xml_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]
[ ]: