.. _lennard-jones-label:
Lennard-Jones fluid
*******************
.. container:: hatnote
The very basics of LAMMPS through a simple example
.. figure:: ../figures/level1/lennard-jones-fluid/lennard-jones-fluid-avatar-dark.webp
:alt: The binary mixture simulated during Tutorial 1. The atoms of type 1 are
represented as small green spheres and the atoms of type 2 as large blue spheres.
:height: 250
:align: right
:class: only-dark
.. figure:: ../figures/level1/lennard-jones-fluid/lennard-jones-fluid-avatar-light.webp
:alt: The binary mixture simulated during Tutorial 1. The atoms of type 1 are
represented as small green spheres and the atoms of type 2 as large blue spheres.
:height: 250
:align: right
:class: only-light
.. container:: abstract
The objective of this tutorial is to perform the simulation of a binary
fluid using LAMMPS.
.. container:: abstract
The system is a Lennard-Jones fluid composed of neutral particles with two
different diameters, contained within a cubic box with periodic boundary conditions
In this tutorial, the temperature of the system is maintained using a Langevin
thermostat :cite:`schneider1978molecular`, and basic quantities are extracted
from the system, including the potential and kinetic energies.
.. container:: abstract
This tutorial illustrates several key ingredients of molecular dynamics
simulations, such as system initialization, energy minimization, integration
of the equations of motion, and trajectory visualization.
.. include:: ../../non-tutorials/needhelp.rst
.. include:: ../../non-tutorials/2Aug2023.rst
My first input
==============
.. container:: justify
To run a simulation using LAMMPS, one needs to write a series of commands
in an input script. For clarity, the input scripts written for this first
tutorial will be divided into five categories which we are going to fill up
one by one.
.. container:: justify
Create a folder, call it *my-first-input/*, and then create a blank
text file in it called *input.lammps*. Copy the following lines
in *input.lammps*, where a line starting with a hash symbol (#)
is a comment ignored by LAMMPS:
.. code-block:: lammps
# PART A - ENERGY MINIMIZATION
# 1) Initialization
# 2) System definition
# 3) Simulation settings
# 4) Visualization
# 5) Run
.. container:: justify
These five categories are not required in every input script, and should not
necessarily be in that exact order. For instance, parts 3 and 4 could be
inverted, or part 4 could be omitted. Note however that LAMMPS reads input
files from top to bottom, therefore the *Initialization* and *System definition*
categories must appear at the top of the input, and the *Run* category at
the bottom.
System initialization
---------------------
.. container:: justify
In the first section of the script, called *Initialization*,
let us indicate to LAMMPS the most basic information
about the simulation, such as:
- the conditions at the boundaries of the box (e.g. periodic or non-periodic),
- the type of atoms (e.g. uncharged single dots or spheres with angular velocities).
Enter the following lines in *input.lammps*:
.. code-block:: lammps
# 1) Initialization
units lj
dimension 3
atom_style atomic
pair_style lj/cut 2.5
boundary p p p
.. container:: justify
The first line, *units lj*, indicates that we want to use the unit system
called *LJ* (Lennard-Jones), in which all quantities are unitless.
.. admonition:: About Lennard-Jones (LJ) units
:class: info
Lennard-Jones (LJ) units are a dimensionless system of units. LJ units are
often used in molecular simulations and theoretical calculations. When using
LJ units:
- energies are expressed in units of :math:`\epsilon`, where :math:`\epsilon`
is the depth of the potential of the LJ interaction,
- distances are expressed in units of :math:`\sigma`, where :math:`\sigma` is
the distance at which the particle-particle potential energy is zero,
- masses are expressed in units of the atomic mass :math:`m`.
All the other quantities are normalized by a combination of :math:`\epsilon`, :math:`\sigma`,
and :math:`m`. For instance, time is expressed in units of :math:`\sqrt{ \epsilon / m \sigma^2}`.
Find details on the |LAMMPS_units|.
.. |LAMMPS_units| raw:: html
LAMMPS website
.. container:: justify
The second line, *dimension 3*, indicates that the simulation
is 3D. The third line, *atom_style atomic*, that the *atomic* style
will be used, therefore each atom is just a dot with a mass.
.. admonition:: About the atom style
:class: info
While we are keeping things as simple as possible in this tutorial,
different *atom_style* will be used in the following tutorials.
Notably, these other atom styles will allow us to create molecules,
i.e. atoms with partial charges and chemical bonds. You can find the complete list
of implemented atom styles from the |atom style page|.
.. |atom style page| raw:: html
atom style page
.. container:: justify
The fourth line, *pair_style lj/cut 2.5*, indicates that atoms
will be interacting through a Lennard-Jones potential with
a cut-off equal to :math:`r_c = 2.5` (unitless)
:cite:`wang2020lennard,fischer2023history`:
.. math::
E_{ij} (r) = 4 \epsilon_{ij} \left[ \left( \dfrac{\sigma_{ij}}{r} \right)^{12}
- \left( \dfrac{\sigma_{ij}}{r} \right)^{6} \right], ~ \text{for} ~ r < r_c,
.. container:: justify
where :math:`r` is the inter-particle distance,
:math:`\epsilon_{ij}` is the depth of potential well that sets the interaction strength, and
:math:`\sigma_{ij}` is the distance parameter or particle effective size.
Here, the indexes *ij* refer to the particle types *i* and *j*.
.. admonition:: About Lennard-Jones potential
:class: info
The Lennard-Jones potential offers a simplified representation that captures
the fundamental aspects of interactions among atoms. It depicts a scenario where two
particles exhibit repulsion at extremely close distances, attraction at moderate
distances, and no interaction at infinite separation. The repulsive part of the
Lennard-Jones potential (i.e. the term :math:`\propto r^{-12}`) is associated
with the Pauli exclusion principle. The attractive part (i.e. the term
in :math:`\propto - r^{-6}`) is linked with the London dispersion forces.
.. container:: justify
The last line, *boundary p p p*, indicates that the periodic boundary
conditions will be used along all three directions of space (the 3 *p* stand
for *x*, *y*, and *z*, respectively).
.. container:: justify
At this point, the *input.lammps* is a LAMMPS input script that does nothing.
You can run it using LAMMPS to verify that the *input* contains
no mistake by typing the following command in the terminal
from the *my-first-input/* folder:
.. code-block:: bw
lmp -in input.lammps
.. container:: justify
Here *lmp* is linked to my compiled LAMMPS version.
Running the previous command should return:
.. code-block:: bw
LAMMPS (2 Aug 2023 - Update 1)
Total wall time: 0:00:00
.. container:: justify
In case there is a mistake in the input script, for example, if
*atom_stile* is written instead of *atom_style*, LAMMPS
gives you an explicit warning:
.. code-block:: bw
LAMMPS (2 Aug 2023 - Update 1)
ERROR: Unknown command: atom_stile atomic (src/input.cpp:232)
Last command: atom_stile atomic
System definition
-----------------
.. container:: justify
Let us fill the *System definition* category of the input script:
.. code-block:: lammps
# 2) System definition
region simulation_box block -20 20 -20 20 -20 20
create_box 2 simulation_box
create_atoms 1 random 1500 341341 simulation_box
create_atoms 2 random 100 127569 simulation_box
.. container:: justify
The first line, *region simulation_box (...)*, creates a region
named *simulation_box* that is a block (i.e. a rectangular cuboid) that
extends from -20 to 20 (no unit) along all 3 directions of space.
.. container:: justify
The second line, *create_box 2 simulation_box*, creates a simulation box based on
the region *simulation_box* with *2* types of atoms.
.. container:: justify
The third line, *create_atoms (...)* creates 1500 atoms of type 1
randomly within the region *simulation_box*. The integer *341341* is a
seed that can be changed in order to create different
initial conditions for the simulation. The fourth line
creates 100 atoms of type 2.
.. container:: justify
If you run LAMMPS, you should see the following information in the terminal:
.. code-block:: bw
(...)
Created orthogonal box = (-20 -20 -20) to (20 20 20)
(...)
Created 1500 atoms
(...)
Created 100 atoms
(...)
.. container:: justify
From what is printed in the terminal, it is clear that
LAMMPS correctly interpreted the commands, and first created
the box with desired dimensions, then 1500 atoms, and then 100
atoms.
Simulation Settings
-------------------
.. container:: justify
Let us fill the *Simulation Settings* category section of
the *input* script:
.. code-block:: lammps
# 3) Simulation settings
mass 1 1
mass 2 1
pair_coeff 1 1 1.0 1.0
pair_coeff 2 2 0.5 3.0
.. container:: justify
The two first commands, *mass (...)*, attribute a mass
equal to 1 (unitless) to both atoms of type 1 and 2.
Alternatively, one could have written
these two commands into one single line: *mass * 1*,
where the star symbol means *all* the atom types of the simulation.
.. container:: justify
The third line, *pair_coeff 1 1 1.0 1.0*, sets the Lennard-Jones
coefficients for the interactions between atoms of type 1,
respectively the energy parameter :math:`\epsilon_{11} = 1.0` and the distance
parameter :math:`\sigma_{11} = 1.0`.
.. container:: justify
Similarly, the last line sets the Lennard-Jones coefficients for
the interactions between atoms of type 2, :math:`\epsilon_{22} = 0.5`,
and :math:`\sigma_{22} = 3.0`.
.. admonition:: About cross parameters
:class: info
By default, LAMMPS calculates the cross coefficients between the different atom types
using geometric average: :math:`\epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}}`,
:math:`\sigma_{ij} = \sqrt{\sigma_{ii} \sigma_{jj}}`. In the present case,
and even without specifying it explicitly, we thus have:
- :math:`\epsilon_{12} = \sqrt{1.0 \times 0.5} = 0.707`, and
- :math:`\sigma_{12} = \sqrt{1.0 \times 3.0} = 1.732`.
When necessary, cross-parameters can be explicitly specified
by adding the following line into the input file: *pair_coeff 1 2 0.707 1.732*.
This can be used for instance to increase the attraction between particles
of type 1 and 2, without affecting the interactions between particles of the same type.
Note that the arithmetic rule, also known as
Lorentz-Berthelot rule :cite:`lorentz1881ueber,berthelot1898melange`, where
:math:`\epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}}`,
:math:`\sigma_{ij} = (\sigma_{ii}+\sigma_{jj})/2`, is more common than the
geometric rule. However, neither the geometric nor the arithmetic rules are
based on rigorous arguments, so here the geometric rule will do just fine.
.. container:: justify
Due to the chosen Lennard-Jones parameters, the two types of particles
are given different effective diameters, as can be seen by plotting
:math:`E_{11} (r)`,
:math:`E_{12} (r)`,
and :math:`E_{22} (r)`.
.. figure:: ../figures/level1/lennard-jones-fluid/lennard-jones.png
:alt: Lennard Jones potential
:class: only-light
:name: fig-lennard-jones
.. figure:: ../figures/level1/lennard-jones-fluid/lennard-jones-dm.png
:alt: Lennard Jones potential
:class: only-dark
.. container:: figurelegend
Figure: The Lennard-Jones potential :math:`E_{ij} (r)`
as a function of the inter-particle distance, where
:math:`i, j = 1 ~ \text{or} ~ 2`. This figure was generated using Python
with Matplotlib Pyplot, and the notebook can be accessed |lennard-jones-pyplot.ipynb|.
The Pyplot parameters used for all figures can be accessed in a |pyplot-perso|.
.. |lennard-jones-pyplot.ipynb| raw:: html
from Github
.. |pyplot-perso| raw:: html
dedicated repository
Energy minimization
-------------------
.. container:: justify
The system is now fully parametrized. Let us fill the two last remaining sections
by adding the following lines into *input.lammps*:
.. code-block:: lammps
# 4) Visualization
thermo 10
thermo_style custom step temp pe ke etotal press
# 5) Run
minimize 1.0e-4 1.0e-6 1000 10000
.. container:: justify
The *thermo* command asks LAMMPS to print
thermodynamic information (e.g. temperature, energy) in the
terminal every given number of steps, here 10 steps.
The *thermo_style custom* requires LAMMPS to print
the system temperature (*temp*), potential energy (*pe*),
kinetic energy (*ke*), total energy (*etotal*),
and pressure (*press*). Finally, the *minimize* command
instructs LAMMPS to perform an energy minimization of the system.
.. admonition:: About energy minimization
:class: info
An energy minimization procedure consists of adjusting the coordinates of
the atoms that are too close to each other until one of the stopping
criteria is reached. By default, LAMMPS uses the conjugate
gradient (CG) algorithm :cite:`hestenes1952methods` (see all the other
implemented methods on the |min_style| page), which runs
until one of the following criteria is reached:
- The change in energy between two iterations is less than 1.0e-4.
- The maximum force between two atoms in the system is lower than 1.0e-6.
- The maximum number of iterations is 1000.
- The maximum number of times the force and the energy have been evaluated is 10000.
.. |min_style| raw:: html
min style
.. container:: justify
Now running the simulation, we can see how the thermodynamic
variables evolve as the simulation progresses:
.. code-block:: bw
Step Temp PotEng KinEng TotEng Press
0 0 78840982 0 78840982 7884122
10 0 169.90532 0 169.90532 17.187291
20 0 -0.22335386 0 -0.22335386 -0.0034892297
30 0 -0.31178296 0 -0.31178296 -0.0027290466
40 0 -0.38135002 0 -0.38135002 -0.0016419218
50 0 -0.42686621 0 -0.42686621 -0.0015219081
60 0 -0.46153953 0 -0.46153953 -0.0010659992
70 0 -0.48581568 0 -0.48581568 -0.0014849169
80 0 -0.51799572 0 -0.51799572 -0.0012995545
(...)
.. container:: justify
These lines give us information about
the progress of the energy minimization. First, at the start
of the simulation (Step 0), the energy in the system is
huge: 78840982 (unitless). This was expected because
the atoms have been created at random positions within the
simulation box and some of them are probably overlapping,
resulting in a large initial energy which is the consequence
of the repulsive part of the Lennard-Jones interaction
potential. As the energy minimization progresses, the energy
rapidly decreases and reaches a negative value, indicating that the atoms have been
displaced at reasonable distances from each other.
.. admonition:: On the temperature during energy minimization
:class: info
As a side note, during energy minimization both temperature and kinetic energy remain equal to
their initial values of 0. This is expected as the conjugate gradient
algorithm only affects the positions of the particles based on the
forces between them, without affecting their velocities.
.. container:: justify
Other useful information has been printed in the terminal, for example, LAMMPS
tells us that the first of the four criteria to be satisfied was the energy:
.. code-block:: bw
Minimization stats:
Stopping criterion = energy tolerance
Molecular dynamics
------------------
.. container:: justify
The system is now ready. Let us continue by completing the input script and
adding commands to perform a molecular dynamics simulation, starting from the
final state of the previous energy minimization step.
.. admonition:: Background Information -- What is molecular dynamics?
:class: info
Molecular dynamics (MD) is based on the numerical solution of the Newtonian
equations of motion for every atom :math:`i`,
.. math::
\sum_{j \ne i} \boldsymbol{F}_{ji} = m_i \times \boldsymbol{a}_i,
where :math:`\sum` is the sum over all the atoms other than :math:`i`,
:math:`\boldsymbol{F}_{ji}` the force between the atom pairs :math:`j-i`,
:math:`m_i` the mass of atom :math:`i`, and :math:`\boldsymbol{a}_i` its acceleration.
The Newtonian equations are solved at every step to predict the
evolution of the positions and velocities of atoms and molecules over
time. Then, the velocity and position of each atom are updated according to the
calculated acceleration, typically using the Verlet algorithm, or similar.
More information can be found in Refs. :cite:`allen2017computer,frenkel2023understanding`.
.. container:: justify
In the same input script, after the *minimization* command, add the following
lines:
.. code-block:: lammps
# PART B - MOLECULAR DYNAMICS
# 4) Visualization
thermo 50
.. container:: justify
Since LAMMPS reads the input from top to bottom, these lines will be
executed after the energy minimization. There is no need to re-initialize
or re-define the system. The *thermo* command is called a second time within
the same input, so the previously entered value of 10 will be replaced by
the value of 50 as soon as *PART B* starts.
.. container:: justify
Then, let us add a second *Run* section:
.. code-block:: lammps
# 5) Run
fix mynve all nve
fix mylgv all langevin 1.0 1.0 0.1 1530917
timestep 0.005
run 10000
.. container:: justify
The *fix nve* is used to update the positions and the velocities of the
atoms in the group *all* at every step. The group *all* is a default group
that contains every atom.
.. container:: justify
The second fix applies a Langevin thermostat to the atoms of the group
*all*, with a desired initial temperature of 1.0 (unitless), and a final
temperature of 1.0 as well :cite:`schneider1978molecular`. A *damping* parameter
of 0.1 is used. The *damping* parameter determines how rapidly the temperature
is relaxed to its desired value. The number *1530917* is a seed, you can
change it to perform statistically independent simulations. Finally, the last
two lines set the value of the *timestep* and the number of steps for the *run*,
respectively, corresponding to a total duration of 50 (unitless).
.. admonition:: What is a fix?
:class: info
In LAMMPS, a *fix* is a command that performs specific tasks during a simulation,
such as imposing constraints, applying forces, or modifying particle properties.
Other LAMMPS-specific terms are defined in the :ref:`glossary-label`.
.. container:: justify
After running the simulation, similar lines should appear in the terminal:
.. code-block:: bw
Step Temp PotEng KinEng TotEng Press
388 0 -0.95476642 0 -0.95476642 -0.000304834
400 0.68476875 -0.90831467 1.0265112 0.11819648 0.023794293
500 0.97168188 -0.56803405 1.4566119 0.88857783 0.02383215
600 1.0364167 -0.44295618 1.5536534 1.1106972 0.027985679
700 1.010934 -0.39601767 1.5154533 1.1194356 0.023064983
800 0.98641731 -0.37866057 1.4787012 1.1000406 0.023131153
900 1.0074571 -0.34951264 1.5102412 1.1607285 0.023520785
(...)
.. container:: justify
The second column shows that the temperature *Temp* starts from 0, but rapidly
reaches the requested value and stabilize itself near :math:`T=1`.
.. container:: justify
From what has been printed in the *log* file, one can plot the potential
energy (:math:`p_\text{e}`) and the kinetic energy (:math:`k_\text{e}`) of
the system over time. The potential energy, :math:`p_\text{e}`, rapidly
decreases during energy minimization. Then, after the molecular dynamics
simulation starts, :math:`p_\text{e}` increases until it reaches a plateau
value of about -0.25. The kinetic energy, :math:`k_\text{e}`, is equal to
zero during energy minimization and then increases during molecular
dynamics until it reaches a plateau value of about 1.5.
.. figure:: ../figures/level1/lennard-jones-fluid/energy.png
:alt: Result tutorial molecular dynamics simulation: Energy plot over time
:class: only-light
.. figure:: ../figures/level1/lennard-jones-fluid/energy-dm.png
:alt: Result tutorial molecular dynamics simulation: Energy plot over time
:class: only-dark
.. container:: figurelegend
Figure: a) Potential energy (:math:`p_\text{e}`) of the binary mixture as a function
of the time :math:`t`. b) Kinetic energy (:math:`k_\text{e}`) as a function of :math:`t`.
Trajectory visualization
------------------------
.. container:: justify
The simulation is running well, but we would like to
visualize the trajectories of the atoms. To do so, we first need
to print the positions of the atoms in a file at a regular interval.
.. container:: justify
Add the following command to the *input.lammps* file, in the *Visualization*
section of *PART B*:
.. code-block:: lammps
dump mydmp all atom 100 dump.lammpstrj
.. container:: justify
Run the *input.lammps* using LAMMPS again. A file named *dump.lammpstrj*
must appear within *my-first-input/*. A *.lammpstrj* file can
be opened using VMD. With Ubuntu/Linux, you can simply execute in the terminal:
.. code-block:: bw
vmd dump.lammpstrj
.. container:: justify
Otherwise, you can open VMD and import the *dump.lammpstrj*
file manually using *File -> New molecule*.
.. container:: justify
By default, you should see a cloud of lines, but you can improve the
representation (see this :ref:`vmd-label` for basic instructions).
.. figure:: ../figures/level1/lennard-jones-fluid/first-input-light.png
:alt: binary fluid simulated by LAMMPS and visualized with VMD
:class: only-light
.. figure:: ../figures/level1/lennard-jones-fluid/first-input-dark.png
:alt: binary fluid simulated by LAMMPS and visualized with VMD
:class: only-dark
.. container:: figurelegend
Figure: View of a slice of the system using VMD, with both
types of atoms represented as spheres. See the corresponding |my_first_input_video|.
.. |my_first_input_video| raw:: html
video
Improving the script
====================
.. container:: justify
Let us improve the input script and perform slightly more advanced operations,
such as imposing a specific initial positions to the atoms, and restarting the simulation
from a previously saved configuration.
Control the initial atom positions
----------------------------------
.. container:: justify
Create a new folder next to *my-first-input/*, and call it *improved-input/*.
Then, create a new input file within *improved-input/* and call it
*input.min.lammps*.
.. container:: justify
Similarly to what has been done previously, copy the following lines
into *input.min.lammps*:
.. code-block:: lammps
# 1) Initialization
units lj
dimension 3
atom_style atomic
pair_style lj/cut 2.5
boundary p p p
.. container:: justify
To create the atoms of types 1 and 2 in two separate regions, let us create
three separate regions: A cubic region for the simulation box and two
additional regions for placing the atoms:
.. code-block:: lammps
# 2) System definition
region simulation_box block -20 20 -20 20 -20 20
create_box 2 simulation_box
region region_cylinder_in cylinder z 0 0 10 INF INF side in
region region_cylinder_out cylinder z 0 0 10 INF INF side out
create_atoms 1 random 1000 341341 region_cylinder_out
create_atoms 2 random 150 127569 region_cylinder_in
.. container:: justify
The *side in* and *side out* keywords
are used to define regions that are respectively inside
and outside of the cylinder of radius 10. Then, copy similar lines
as previously into *input.min.lammps*:
.. code-block:: lammps
# 3) Simulation settings
mass 1 1
mass 2 1
pair_coeff 1 1 1.0 1.0
pair_coeff 2 2 0.5 3.0
# 4) Visualization
thermo 10
thermo_style custom step temp pe ke etotal press
dump mydmp all atom 10 dump.min.lammpstrj
# 5) Run
minimize 1.0e-4 1.0e-6 1000 10000
write_data minimized_coordinate.data
.. container:: justify
The main novelty, compared to the previous input script, is the *write_data*
command. This command is used to print the final state of the simulation in
a file named *minimized_coordinate.data*. Note that the *write_data* command
is placed after the *minimize* command. This *.data* file will be used later
to restart the simulation from the final state of the energy minimization step.
.. container:: justify
Run the *input.min.lammps* script using LAMMPS.
.. container:: justify
As soon as the simulation starts, a new dump file named *dump.min.lammpstrj*
must appear in the folder. This *.lammpstrj* can be used to visualize the
atom's trajectories during minimization using VMD. At the end of the simulation,
a file named *minimized_coordinate.data* is created by LAMMPS.
.. container:: justify
If you open *minimized_coordinate.data* with a text editor, you can see that
it contains all the information necessary to restart the simulation, such as
the number of atoms, the box size, the *masses*, and the *pair_coeffs*:
.. code-block:: lammps
1150 atoms
2 atom types
-20 20 xlo xhi
-20 20 ylo yhi
-20 20 zlo zhi
Masses
1 1
2 1
Pair Coeffs # lj/cut
1 1 1
2 0.5 3
(...)
.. container:: justify
The *minimized_coordinate.data* file also contains the final
positions of the atoms:
.. code-block:: lammps
(...)
Atoms # atomic
970 1 4.4615279184230525 -19.88248310680258 -19.497251754277872 0 0 0
798 1 1.0773937287460968 -17.57843015813612 -19.353475858951473 0 0 0
21 1 -17.542385434367777 -16.647460269156497 -18.93914807895693 0 0 0
108 1 -15.96241088290946 -15.956274144833264 -19.016419910024062 0 0 0
351 1 0.08197850837343444 -16.852380573900156 -19.28249747472579 0 0 0
402 1 -5.270160783673711 -15.592291204068946 -19.6382667867645 0 0 0
(...)
.. container:: justify
The first five columns of the *Atoms* section correspond (from left to right)
to the atom indexes (from 1 to the total number of atoms, 1150), the atom types (1 or 2
here), and the atoms positions :math:`x`, :math:`y`, :math:`z`. The last
three columns are image flags that keep track of which atoms crossed the
periodic boundary.
Restarting from a saved configuration
-------------------------------------
.. container:: justify
Let us create a new input file and start a molecular dynamics simulation
directly from the previously saved configuration. Within *improved-input/*,
create a new file named *input.md.lammps* and copy the same lines as previously:
.. code-block:: lammps
# 1) Initialization
units lj
dimension 3
atom_style atomic
pair_style lj/cut 2.5
boundary p p p
.. container:: justify
Here, instead of creating a new region and adding atoms to it, we can simply
import the previously saved configuration by adding the following command
to input.md.lammps:
.. code-block:: lammps
# 2) System definition
read_data minimized_coordinate.data
.. container:: justify
By visualizing the previously generated *dump.min.lammpstrj*
file, you may have noticed that some atoms have moved from one region to
the other during minimization. To start the simulation from a clean slate, with
only atoms of type 2 within the cylinder and atoms of type 1 outside the
cylinder, let us delete the misplaced atoms by adding the following commands
to *input.md.lammps*:
.. code-block:: lammps
read_data minimized_coordinate.data
region region_cylinder_in cylinder z 0 0 10 INF INF side in
region region_cylinder_out cylinder z 0 0 10 INF INF side out
group group_type_1 type 1
group group_type_2 type 2
group group_region_in region region_cylinder_in
group group_region_out region region_cylinder_out
group group_type_1_in intersect group_type_1 group_region_in
group group_type_2_out intersect group_type_2 group_region_out
delete_atoms group group_type_1_in
delete_atoms group group_type_2_out
.. container:: justify
The two first *region* commands recreate
the previously defined regions, which is necessary since
regions are not saved by the *write_data* command.
.. container:: justify
The first two *group* commands are used to create groups containing
all the atoms of type 1 and all the atoms of type 2, respectively.
The next two *group* commands create atom groups based on their
positions at the beginning of the simulation, i.e. when the commands
are being read by LAMMPS. The last two *group* commands create atom groups
based on the intersection between the previously defined groups.
.. container:: justify
Finally, the two *delete_atoms* commands delete the
atoms of type 1 that are located within the cylinder and the atoms of
type 2 that are located outside the cylinder, respectively.
.. container:: justify
When you run the *input.md.lammps* input using LAMMPS, you
can see in the *log* file how many atoms are in each group,
and how many atoms have been deleted:
.. code-block:: bw
1000 atoms in group group_type_1
150 atoms in group group_type_2
149 atoms in group group_region_in
1001 atoms in group group_region_out
0 atoms in group group_type_1_in
1 atoms in group group_type_2_out
Deleted 0 atoms, new total = 1150
Deleted 1 atoms, new total = 1149
.. container:: justify
Add the following lines into *input.md.lammps*.
Note the absence of *Simulation settings* section,
because the settings are taken from the *.data* file.
.. code-block:: lammps
# 4) Visualization
thermo 1000
dump mydmp all atom 1000 dump.md.lammpstrj
.. container:: justify
Let us extract the number of atoms of each type
inside the cylinder as a function of time, by
adding the following commands to *input.md.lammps*:
.. code-block:: lammps
variable n_type1_in equal count(group_type_1,region_cylinder_in)
variable n_type2_in equal count(group_type_2,region_cylinder_in)
fix myat1 all ave/time 10 200 2000 v_n_type1_in &
file output-population1vstime.dat
fix myat2 all ave/time 10 200 2000 v_n_type2_in &
file output-population2vstime.dat
.. container:: justify
The two *variables* are used to count the number of atoms of a specific
group in the *region_cylinder_in* region.
.. container:: justify
The two *fix ave/time* are calling the previously defined variables and are
printing their values into text files. By using *10 200 2000*, variables are
evaluated every 10 steps, averaged 200 times, and printed in the *.dat* files
every 2000 steps.
.. container:: justify
In addition to counting the atoms in each region, let us also extract the
coordination number per atom between atoms of types 1 and 2. The
coordination number is a measure of the average number of type 2 atoms
in the vicinity of type 1 atoms, serving as a good indicator of
the degree of mixing in a binary mixture. Add the following lines into
*input.md.lammps*:
.. code-block:: lammps
compute coor12 group_type_1 coord/atom cutoff 2.0 group group_type_2
compute sumcoor12 all reduce ave c_coor12
fix myat3 all ave/time 10 200 2000 &
c_sumcoor12 file coordinationnumber12.dat
.. container:: justify
The *compute ave* is used to average the per atom
coordination number that is calculated by the *coord/atom* compute.
This averaging is necessary as *coord/atom* returns an array where each value corresponds
to a certain couple of atoms i-j. Such an array can't be printed by *fix ave/time*.
Finally, let us complete the script by adding the following lines
to *input.md.lammps*:
.. code-block:: lammps
# 5) Run
velocity all create 1.0 4928459 mom yes rot yes dist gaussian
fix mynve all nve
fix mylgv all langevin 1.0 1.0 0.1 1530917 zero yes
timestep 0.005
run 300000
write_data mixed.data
.. container:: justify
There are a few differences from the previous simulation.
First, the *velocity create* command attributes an initial velocity to every
atom. The initial velocity is chosen so that the average initial
temperature is equal to 1 (unitless). The additional
keywords ensure that no linear momentum (*mom yes*) and no angular
momentum (*rot yes*) are given to the system and that the generated
velocities are distributed as a Gaussian. Another improvement
is the *zero yes* keyword in the Langevin thermostat, which
ensures that the total random force applied to the atoms is equal to zero.
.. container:: justify
Run *input.md.lammps* using LAMMPS and visualize the trajectory
using VMD.
.. figure:: ../figures/level1/lennard-jones-fluid/mixing-vmd-light.png
:alt: LAMMPS VMD tutorial molecular dynamics simulation
:class: only-light
.. figure:: ../figures/level1/lennard-jones-fluid/mixing-vmd-dark.png
:alt: LAMMPS VMD tutorial molecular dynamics simulation
:class: only-dark
.. container:: figurelegend
Figure: Evolution of the system during mixing. The three snapshots show
respectively the system at :math:`t=0` (left panel),
:math:`t=75` (middle panel), and :math:`t=1500` (right panel).
.. container:: justify
After running *input.md.lammps* using LAMMPS, you can observe the number
of atoms in each region from the generated data files, as
well as the evolution of the coordination number due to mixing:
.. figure:: ../figures/level1/lennard-jones-fluid/mixing.png
:alt: Result tutorial molecular dynamics simulation: Energy plot over time
:class: only-light
.. figure:: ../figures/level1/lennard-jones-fluid/mixing-dm.png
:alt: Result tutorial molecular dynamics simulation: Energy plot over time
:class: only-dark
.. container:: figurelegend
Figure: Evolution of the number of atoms within the *region_cylinder_in* region
as a function of time (a), and evolution of the coordination number
between atoms of types 1 and 2 (b).
.. include:: ../../non-tutorials/accessfile.rst
Going further with exercises
============================
.. include:: ../../non-tutorials/link-to-solutions.rst
Solve Lost atoms error
----------------------
.. container:: justify
For this exercise, the following input script is provided:
.. code-block:: lammps
units lj
dimension 3
atom_style atomic
pair_style lj/cut 2.5
boundary p p p
region simulation_box block -20 20 -20 20 -20 20
create_box 1 simulation_box
create_atoms 1 random 1000 341841 simulation_box
mass 1 1
pair_coeff 1 1 1.0 1.0
dump mydmp all atom 100 dump.lammpstrj
thermo 100
thermo_style custom step temp pe ke etotal press
fix mynve all nve
fix mylgv all langevin 1.0 1.0 0.1 1530917
timestep 0.005
run 10000
.. container:: justify
As it is, this input returns one of the most common
error that you will encounter using LAMMPS:
.. code-block:: bash
ERROR: Lost atoms: original 1000 current 984
.. container:: justify
The goal of this exercise is to fix the *Lost atoms* error without
using any other command than the ones already present. You can
only play with the values of the parameters and/or replicate every
command as many times as needed.
.. admonition:: Note
:class: info
This script is failing because particles are created randomly in space, some
of them are likely overlapping, and no energy minimization is performed prior
to start the molecular dynamics simulation.
Create a demixed dense phase
----------------------------
.. container:: justify
Starting from one of the *input* created in this tutorial, fine-tune the
parameters such as particle numbers and interaction to create a simulation
with the following properties:
- the density in particles must be high,
- both particles of type 1 and 2 must have the same size,
- particles of type 1 and 2 must demix.
.. figure:: ../figures/level1/lennard-jones-fluid/demixing-light.png
:alt: VMD/LAMMPS exercise molecular dynamics simulation: demixing lennard
jones fluids
:class: only-light
.. figure:: ../figures/level1/lennard-jones-fluid/demixing-dark.png
:alt: VMD/LAMMPS exercise molecular dynamics simulation: demixing lennard
jones fluids
:class: only-dark
.. container:: figurelegend
Figure: Snapshots taken at different times showing the particles of type 1
and type 2 progressively demixing and forming large demixed areas.
.. admonition:: Hint
:class: info
An easy way to create a dense phase is to allow the box dimensions to relax
until the vacuum disappears. You can do that by replacing the *fix nve* with *fix nph*.
From atoms to molecules
-----------------------
.. container:: justify
Add a bond between particles of *type 2* to create dumbbell molecules instead
of single particles.
.. figure:: ../figures/level1/lennard-jones-fluid/dumbell-dark.png
:alt: Dumbbell Lennard-Jones molecules simulated using LAMMPS
:class: only-dark
.. figure:: ../figures/level1/lennard-jones-fluid/dumbell-light.png
:alt: Dumbbell Lennard-Jones molecules simulated using LAMMPS
:class: only-light
.. container:: figurelegend
Figure: Dumbbell molecules made of 2 large spheres mixed with smaller
particles (small spheres). See the corresponding |dumbell_video|.
.. |dumbell_video| raw:: html
video
.. container:: justify
Similarly to the dumbbell molecules, create a small polymer,
i.e. a long chain of particles linked by bonds and angles.
.. figure:: ../figures/level1/lennard-jones-fluid/polymer-dark.png
:alt: Polymer Lennard-Jones molecules simulated using LAMMPS
:class: only-dark
.. figure:: ../figures/level1/lennard-jones-fluid/polymer-light.png
:alt: Polymer Lennard-Jones molecules simulated using LAMMPS
:class: only-light
.. container:: figurelegend
Figure: A single small polymer molecule made of 9 large spheres mixed with
smaller particles. See the corresponding |polymer_video|.
.. |polymer_video| raw:: html
video
.. admonition:: Hints
:class: info
.. container:: justify
Use a *molecule template* to easily insert as many atoms connected
by bonds (i.e. molecules) as you want. A molecule template typically
begins as follows:
.. code-block:: lammps
2 atoms
1 bonds
Coords
1 0.5 0 0
2 -0.5 0 0
(...)
.. container:: justify
A bond section also needs to be added, see this
|molecule_template_lammps| for details on the formatting of a
molecule template.
.. |molecule_template_lammps| raw:: html
page