.. _carbon-nanotube-label:
Pulling on a carbon nanotube
****************************
.. container:: hatnote
Stretching a carbon nanotube until it breaks
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/CNT_dark.webp
:alt: carbon nanotube image in vacuum
:height: 250
:align: right
:class: only-dark
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/CNT_light.webp
:alt: carbon nanotube image in vacuum
:height: 250
:align: right
:class: only-light
.. container:: abstract
The objective of this tutorial is to impose the deformation
of a carbon nanotube (CNT) using LAMMPS.
.. container:: abstract
In this tutorial, a small carbon nanotube (CNT) is simulated
within an empty box using LAMMPS. An external
force is imposed on the CNT, and its deformation is measured over time.
.. container:: abstract
To illustrate the difference between classical and reactive force fields,
this tutorial is divided into two parts. Within the first part, a classical
force field is used and the bonds between the atoms of the CNT are
unbreakable. Within the second part, a reactive force field
(named AIREBO :cite:`stuart2000reactive`) is used, allowing for the breaking
of chemical bonds when the CNT undergoes strong deformation.
.. include:: ../../non-tutorials/recommand-lj.rst
.. include:: ../../non-tutorials/needhelp.rst
.. include:: ../../non-tutorials/2Aug2023.rst
Unbreakable bonds
=================
.. container:: justify
With most classical molecular dynamics force fields, the chemical bonds
between the atoms are set at the start of the simulation. Regardless of the
forces applied to the atoms during the simulations, the bonds remain intact.
The bonds between neighbor atoms typically consist of springs with
given equilibrium distances :math:`r_0` and a constant :math:`k_b`:
:math:`U_b = k_b \left( r - r_0 \right)^2`.
Additionally, angular and dihedral constraints are usually applied to maintain
the relative orientations of neighbor atoms.
Create topology with VMD
------------------------
.. container:: justify
The first part of this tutorial is dedicated to creating
the initial topology with VMD. You can skip this part by
downloading directly the CNT topology |download_cnt_molecular_data|,
and continue with the LAMMPS part of the tutorial.
.. |download_cnt_molecular_data| raw:: html
here
.. admonition:: Why use a preprocessing tool?
:class: info
When the system has a complex topology, like is the case of a CNT,
it is better to use an external preprocessing tool to create it as it would be
difficult (yet not impossible) to place the atoms in their correct position
using only LAMMPS commands. Many preprocessing tools exist, see
this |prepross| on the LAMMPS website.
.. |prepross| raw:: html
non-exhaustive list
.. container:: justify
Open VMD, go to Extensions, Modeling, and then Nanotube Builder.
A window named Carbon Nanostructures opens up, allowing us to choose
between generating a sheet or a nanotube, made either of graphene or
of Boron Nitride (BN). For this tutorial, let us generate a carbon nanotube.
Keep all default values, and click on *Generate Nanotube*.
.. container:: justify
At this point, this is not a molecular dynamics simulation,
but a cloud of unconnected dots. In the VMD terminal, set the
box dimensions by typing the following commands in the VMD terminal:
.. code-block:: bw
molinfo top set a 80
molinfo top set b 80
molinfo top set c 80
.. container:: justify
The values of 80 in each direction have been chosen
so that the box is much larger than the carbon nanotube.
.. container:: justify
To generate the initial LAMMPS data file, let us use *Topotools*:
to generate the LAMMPS data file, enter the following command
in the VMD terminal:
.. code-block:: bw
topo writelammpsdata cnt_molecular.data molecular
.. container:: justify
Here *molecular* refers to the LAMMPS *atom_style*, and
*cnt_molecular.data* is the name of the file.
.. admonition:: About TopoTools
:class: info
*Topotools* deduces the location of bonds, angles,
dihedrals, and impropers from the respective positions of the atoms,
and generates a *.data* file that can be read by LAMMPS :cite:`kohlmeyer2017topotools`.
More details about *Topotools* can be found on the
personal page of |Axel_webpage|.
.. |Axel_webpage| raw:: html
Axel Kohlmeyer
.. container:: justify
The parameters of the constraints (bond length,
dihedral coefficients, etc.) will be given later.
A new file named *cnt_molecular.data* has been created, it starts
like that:
.. code-block:: lammps
700 atoms
1035 bonds
2040 angles
4030 dihedrals
670 impropers
1 atom types
1 bond types
1 angle types
1 dihedral types
1 improper types
-40.000000 40.000000 xlo xhi
-40.000000 40.000000 ylo yhi
-12.130411 67.869589 zlo zhi
(...)
.. container:: justify
The *cnt_molecular.data* file contains information
about the positions of the carbon atoms, as well as the
identity of the atoms that are linked by *bonds*, *angles*, *dihedrals*,
and *impropers* constraints.
.. container:: justify
Save the *cnt_molecular.data* file in a folder named *unbreakable-bonds/*.
The LAMMPS input
----------------
.. container:: justify
Create a new text file within *unbreakable-bonds/* and name
it *input.lammps*. Copy the following lines into it:
.. code-block:: lammps
variable T equal 300
units real
atom_style molecular
boundary f f f
pair_style lj/cut 14
bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style harmonic
special_bonds lj 0.0 0.0 0.5
read_data cnt_molecular.data
.. container:: justify
The chosen unit system is *real* (therefore distances are in Ångstrom and time in femtosecond),
the *atom_style* is molecular (therefore atoms are dots that can be bonded with each other),
and the boundary conditions are fixed. The boundary conditions
do not matter here, as the box boundaries were placed far from the CNT.
.. container:: justify
Just like in :ref:`lennard-jones-label`,
the pair style is *lj/cut* (i.e. a Lennard-Jones potential
with a short-range cutoff) with parameter 14, which means that only the atoms
closer than 14 Ångstroms from each other interact through a Lennard-Jones
potential.
.. container:: justify
The *bond_style*, *angle_style*, *dihedral_style*, and *improper_style* commands specify the
different potentials used to restrain the relative positions of the
atoms. For more details about the potentials used here, you can have a look
at the LAMMPS website, see for example
the page of the |OPLS| :cite:`jorgensen1988opls`.
.. |OPLS| raw:: html
OPLS dihedral style
.. container:: justify
The last command, *read_data*, imports the *cnt_molecular.data* file
previously generated with VMD, which contains
information about the box size, atom positions, etc.
.. admonition:: About interaction between neighbors atoms
:class: info
Atoms connected by a bond do not typically interact through
Lennard-Jones interaction. Therefore, atoms that are
bounded must be excluded from the Lennard-Jones potential calculation.
Here, this is done by the *special_bonds* command. The three numbers of the
*special_bonds* command are weighting factors for the
Lennard-Jones interaction between atoms connected by a bond
(respectively directly bounded :math:`C-C`, separated by two bonds :math:`C-C-C`,
and separated by three bonds :math:`C-C-C-C`). For instance, the
first weighting factor, with a value of 0, imposes that two atoms connected
by a bond do not interact through a Lennard-Jones potential (therefore
they only interact through the harmonic potential that bonds the atoms
of the graphene).
.. container:: justify
We need to specify the parameters of both bonded and
non-bonded potentials. Here, the parameters are taken from the OPLS-AA
(Optimised Potentials for Liquid Simulations-All-Atom) force
field :cite:`jorgensenDevelopmentTestingOPLS1996`.
Create a new text file in the *unbreakable-bonds/*
folder and name it *parm.lammps*. Copy the following lines into it:
.. code-block:: lammps
pair_coeff 1 1 0.066 3.4
bond_coeff 1 469 1.4
angle_coeff 1 63 120
dihedral_coeff 1 0 7.25 0 0
improper_coeff 1 5 180
.. container:: justify
The *pair_coeff* command sets the parameters for the non-bonded Lennard-Jones
interaction :math:`\epsilon_{11} = 0.066 \, \text{kcal/mol}`
and :math:`\sigma_{11} = 3.4 \, \text{Å}` for the only type of atom of the
simulation; the carbon atom of type 1.
.. container:: justify
The *bond_coeff* provides the equilibrium distance :math:`r_0= 1.4 \, \text{Å}` as
well as the spring constant :math:`k_b = 469 \, \text{kcal/mol/Å}^2` for the harmonic
potential imposed between two neighboring carbon atoms,
where the potential is :math:`U_b = k_b ( r - r_0)^2`. The
*angle_coeff* gives the equilibrium angle :math:`\theta_0` and
constant for the potential between three neighbor atoms :
:math:`U_\theta = k_\theta ( \theta - \theta_0)^2`. The *dihedral_coeff*
and *improper_coeff* gives the potential for the constraints
between 4 atoms.
.. container:: justify
The file *parm.lammps* is included in the simulation by adding the
following line into the *input.lammps* file:
.. code-block:: lammps
include parm.lammps
Prepare the initial state
-------------------------
.. container:: justify
Before starting the molecular dynamics simulation, let us make sure that we
start from a clean initial state by recentering the CNT at the origin (0, 0, 0).
In addition, let us make sure that the box boundaries are symmetric with
respect to (0, 0, 0), which is not initially the case, as seen in *cnt_molecular.data*:
.. code-block:: lammps
-40.000000 40.000000 xlo xhi
-40.000000 40.000000 ylo yhi
-12.130411 67.869589 zlo zhi
.. container:: justify
Let us recenter the CNT by adding the following lines
to *input.lammps*:
.. code-block:: lammps
group carbon_atoms type 1
variable carbon_xcm equal -1*xcm(carbon_atoms,x)
variable carbon_ycm equal -1*xcm(carbon_atoms,y)
variable carbon_zcm equal -1*xcm(carbon_atoms,z)
displace_atoms carbon_atoms &
move ${carbon_xcm} ${carbon_ycm} ${carbon_zcm}
.. container:: justify
The first command includes all the atoms of type 1
(i.e. all the atoms here) in a group named *carbon_atoms*.
The 3 variables, *carbon_xcm*, *carbon_ycm*, and *carbon_zcm*
are used to measure the current position of the group *carbon_atoms* along
all 3 directions, respectively. Then, the *displace_atoms*
command moves the group *carbon_atoms*, ensuring that its center of mass
is located at the origin (0, 0, 0).
.. container:: justify
Let us also change the box boundaries by adding the following line into *input.lammps*:
.. code-block:: lammps
change_box all x final -40 40 y final -40 40 z final -40 40
.. admonition:: Note
:class: info
Such a cleaner and more symmetrical initial state can simplify
future data analysis, but won't make any difference to
the molecular dynamics.
.. container:: justify
A displacement will be imposed on the edges of the CNT. To do so, let us isolate the
atoms from the two edges and place them into groups named *rtop*
and *rbot*, respectively. Add the following lines into *input.lammps*:
.. code-block:: lammps
variable zmax equal bound(carbon_atoms,zmax)-0.5
variable zmin equal bound(carbon_atoms,zmin)+0.5
region rtop block INF INF INF INF ${zmax} INF
region rbot block INF INF INF INF INF ${zmin}
region rmid block INF INF INF INF ${zmin} ${zmax}
.. container:: justify
The variable :math:`z_\mathrm{max}` corresponds to
the coordinate of the last atoms along :math:`z` minus 0.5
Ångstroms, and :math:`z_\mathrm{min}` to the coordinate of
the first atoms along :math:`z` plus 0.5 Ångstroms. Then, three
regions are defined and correspond respectively to: :math:`z < z_\mathrm{min}`,
(*rbot*, for region bottom)
:math:`z_\mathrm{min} > z > z_\mathrm{max}`
(*rmid*, for region middle), and
:math:`z > z_\mathrm{max}`
(*rtop*, for region top).
.. container:: justify
Finally, let us define 3 groups of atoms
corresponding to the atoms located in each of the three regions,
respectively, by adding to *input.lammps*:
.. code-block:: lammps
group carbon_top region rtop
group carbon_bot region rbot
group carbon_mid region rmid
.. container:: justify
The atoms of the edges as selected within the *carbon_top*
and *carbon_bot* groups can be represented with a different color.
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/colored-edge-undef-dark.png
:alt: CNT in graphene in vacuum image VMD
:class: only-dark
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/colored-edge-undef-light.png
:alt: CNT in graphene in vacuum image VMD
:class: only-light
.. container:: figurelegend
Figure: CNT with the atoms from the *carbon_mid* group in pink,
and the atoms from the *carbon_top* and the *carbon_bot* groups in white.
.. container:: justify
When running a simulation, the number of atoms in each group is printed in
the terminal (and in the *log.lammps* file). Always make sure that the number
of atoms in each group corresponds to what is expected, just like here:
.. code-block:: bash
10 atoms in group carbon_top
10 atoms in group carbon_bot
680 atoms in group carbon_mid
.. container:: justify
Finally, to start from a less ideal state and create a system with some defects,
let us randomly delete a small fraction of the carbon atoms.
To avoid deleting atoms that are too close to the edges,
let us define a new region name *rdel* that
starts :math:`2\,Å` from the CNT edges.
.. code-block:: lammps
variable zmax_del equal ${zmax}-2
variable zmin_del equal ${zmin}+2
region rdel block INF INF INF INF ${zmin_del} ${zmax_del}
group rdel region rdel
delete_atoms random fraction 0.02 no rdel NULL 482793 bond yes
.. container:: justify
The *delete_atoms* command randomly deletes :math:`2\,\%` of the atoms
from the *rdel* group (i.e. about 10 atoms).
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/colored-edge-deleted-dark.png
:alt: CNT in graphene in vacuum image VMD with deleted atoms
:class: only-dark
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/colored-edge-deleted-light.png
:alt: CNT in graphene in vacuum image VMD with deleted atoms
:class: only-light
.. container:: figurelegend
Figure: CNT with 10 randomly deleted atoms. The 10 deleted atoms were chosen randomly
from the central part of the CNT.
The molecular dynamics
----------------------
.. container:: justify
Let us specify the thermalization and the dynamics of the
system. Add the following lines into *input.lammps*:
.. code-block:: lammps
reset_atoms id sort yes
velocity carbon_mid create ${T} 48455 mom yes rot yes
fix mynve all nve
compute Tmid carbon_mid temp
fix myber carbon_mid temp/berendsen ${T} ${T} 100
fix_modify myber temp Tmid
.. container:: justify
Re-setting the atom IDs is necessary before using the *velocity* command,
this is done by the *reset_atoms* command.
.. container:: justify
The *velocity* command gives initial velocities to
the atoms of the middle group *carbon_mid*, ensuring an initial temperature
of 300 K for these atoms with no overall translational momentum, *mom yes*,
nor rotational momentum, *rot yes*.
.. container:: justify
The *fix nve* is applied to all atoms so that all atom positions are
recalculated at every step, and a *Berendsen* thermostat is applied to the atoms
of the group *carbon_mid* only :cite:`berendsen1984molecular`. The *fix_modify myber*
ensures that the *fix Berendsen* uses the temperature of the group *carbon_mid* as an
input, instead of the temperature of the whole system. This is necessary
to make sure that the frozen edges won't bias the temperature. Note that the atoms
of the edges do not need a thermostat because their motion will
be restrained, see below.
.. admonition:: Deal with semi-frozen system
:class: info
Always be careful when part of a system is frozen,
as is the case here. When some atoms are frozen, the total temperature
of the system is effectively lower than the applied temperature
because the frozen atoms have no thermal motion (their temperature
is therefore :math:`0\,\text{K}`).
.. container:: justify
To restrain the motion of the atoms at the edges, let us add the
following commands to *input.lammps*:
.. code-block:: lammps
fix mysf1 carbon_top setforce 0 0 0
fix mysf2 carbon_bot setforce 0 0 0
velocity carbon_top set 0 0 0
velocity carbon_bot set 0 0 0
.. container:: justify
The two *setforce* commands cancel the forces applied on the
atoms of the two edges, respectively. The cancellation of the forces
is done at every step, and along all 3 directions of space, :math:`x`, :math:`y`,
and :math:`z`, due to the use of *0 0 0*. The two *velocity* commands
set the initial velocities along :math:`x`,
:math:`y`, and :math:`z` to 0 for the atoms of *carbon_top*
and *carbon_bot*, respectively.
.. container:: justify
As a consequence of these last four commands, the atoms of the edges will remain
immobile during the simulation (or at least they would if no other command was
applied to them).
.. admonition:: On imposing a constant velocity to a system
:class: info
The *velocity set* commands impose the velocity of a group of atoms at the start of
a run but do not enforce the velocity during the entire simulation.
When *velocity set* is used in combination with *setforce 0 0 0*,
as is the case here, the atoms
won't feel any force during the entire simulation. According to the Newton equation,
no force means no acceleration, meaning that the initial velocity will persist
during the entire simulation, thus producing a constant velocity motion.
Data extraction
---------------
.. container:: justify
Next, in order to measure the strain and stress suffered by the
CNT, let us extract the distance :math:`L` between
the two edges as well as the force applied on the edges.
.. code-block:: lammps
variable L equal xcm(carbon_top,z)-xcm(carbon_bot,z)
fix at1 all ave/time 10 10 100 v_L file output_cnt_length.dat
fix at2 all ave/time 10 10 100 f_mysf1[1] f_mysf2[1] &
file output_edge_force.dat
.. container:: justify
Let us also add a command to print the atom coordinates in a
*lammpstrj* file every 1000 steps.
.. code-block:: lammps
dump mydmp all atom 1000 dump.lammpstrj
.. container:: justify
Finally, let us check the temperature of the non-frozen group over time
by printing it using a *fix ave/time* command:
.. code-block:: lammps
fix at3 all ave/time 10 10 100 c_Tmid &
file output_temperature_middle_group.dat
.. admonition:: About extracting quantity from variable compute or fix
:class: info
Notice that the values of the force on each edge are
extracted from the two *fix setforce* *mysf1* and *mysf2*, simply by
calling them using *f_*, the same way variables are called
using *v_* and computes are called using *c_*.
.. container:: justify
Let us run a small equilibration step to bring the system
to the required temperature before applying any deformation:
.. code-block:: lammps
thermo 100
thermo_modify temp Tmid
timestep 1.0
run 5000
.. container:: justify
With the *thermo_modify* command, we specify to LAMMPS that we
want the temperature :math:`T_\mathrm{mid}` to be printed in
the terminal, not the temperature of the entire system
(because of the frozen edges, the temperature of
the entire system is not relevant).
.. container:: justify
Let us impose a constant velocity deformation on the CNT by combining
the *velocity set* command with previously defined *fix setforce*.
Add the following lines in the *input.lammps* file,
right after the last *run 5000* command:
.. code-block:: lammps
# 2*0.0005 A/fs = 0.001 A/fs = 100 m/s
velocity carbon_top set 0 0 0.0005
velocity carbon_bot set 0 0 -0.0005
run 10000
.. container:: justify
The chosen velocity for the deformation is :math:`100\,\text{m/s}`.
The length :math:`L` of the CNT increase linearly over
time for :math:`t > 5\,\text{ps}`, as expected from the imposed constant velocity.
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/length-unbreakable.png
:alt: length of the CNT with time - lammps molecular dynamics
:class: only-light
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/length-unbreakable-dm.png
:alt: length of the CNT with time - lammps molecular dynamics
:class: only-dark
.. container:: figurelegend
Figure: Evolution of the length :math:`L` of the CNT with time.
The CNT starts deforming at :math:`t = 5\,\text{ps}`.
.. container:: justify
The energy, which can be accessed from the log file, shows a non-linear
increase with time once the deformation starts,
which is expected from the typical dependency of bond energy with
bond distance :math:`U_b = k_b \left( r - r_0 \right)^2`.
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/energy-unbreakable-dm.png
:alt: energy of the CNT with time - lammps molecular dynamics
:class: only-dark
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/energy-unbreakable.png
:alt: energy of the CNT with time - lammps molecular dynamics
:class: only-light
.. container:: figurelegend
Figure: Evolution of the total energy of the system with time.
The CNT starts deforming at :math:`t = 5\,\text{ps}`.
.. container:: justify
As always, is it important to ensure that the simulation
behaves as expected by opening the *dump.lammpstrj* file with VMD.
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/colored-edge-def-dark.png
:alt: CNT in graphene in vacuum image VMD before and after deformation
:class: only-dark
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/colored-edge-def-light.png
:alt: CNT in graphene in vacuum image VMD before and after deformation
:class: only-light
.. container:: figurelegend
Figure: CNT before (top) and after (bottom) deformation. See the corresponding |unbreakable_cnt_video|.
.. |unbreakable_cnt_video| raw:: html
video
Breakable bonds
===============
.. container:: justify
When using a classical force field, as we just did, the bonds between the atoms
are non-breakable. Let us perform a similar simulation and deform a small CNT again,
but this time using a reactive force field that allows for the bonds to break
if the applied deformation is large enough.
Input file initialization
-------------------------
.. container:: justify
Create a second folder named *breakable-bonds/* next to *unbreakable-bonds/*,
and create a new input file in it called *input.lammps*. Type into input.lammps*:
.. code-block:: lammps
# Initialisation
variable T equal 300
units metal
atom_style atomic
boundary p p p
pair_style airebo 2.5 1 1
.. container:: justify
The first difference with the previous part
is the unit system, here *metal* instead of *real*, a choice
that is imposed by the AIREBO force field. A second difference
is the use of the *atom_style atomic* instead of *molecular*,
since no explicit bond information is required with AIREBO.
.. admonition:: About metal units
:class: info
With the *metal* units system of LAMMPS, the time is in pico second,
distances are in Ångstrom, and the energy is in eV.
Adapt the topology file
-----------------------
.. container:: justify
Since *bond*, *angle*, and *dihedral* do not need to be explicitly
set when using AIREBO, some small changes need to be made to the
previously generated *.data* file.
.. container:: justify
Duplicate the previous file *cnt_molecular.data*, name the copy *cnt_atom.data*,
and place it within *breakable-bonds/*. Then, remove all bond, angle, and dihedral
information from *cnt_atom.data*. Also, remove the second column of the
*Atoms* table, so that the *cnt_atom.data* looks like the following:
.. code-block:: lammps
700 atoms
1 atom types
-40.000000 40.000000 xlo xhi
-40.000000 40.000000 ylo yhi
-12.130411 67.869589 zlo zhi
Masses
1 12.010700 # CA
Atoms # atomic
1 1 5.162323 0.464617 8.843235 # CA CNT
2 1 4.852682 1.821242 9.111212 # CA CNT
(...)
.. container:: justify
In addition, remove the *Bonds* table that is placed right after the
*Atoms* table (near line 743), as well as the *Angles*, *Dihedrals*,
and *Impropers* tables. The last lines of the file should look like this:
.. code-block:: lammps
(...)
697 1 4.669892 -2.248901 45.824036 # CA CNT
698 1 5.099893 -0.925494 46.092010 # CA CNT
699 1 5.162323 -0.464617 47.431896 # CA CNT
700 1 5.099893 0.925494 47.699871 # CA CNT
.. container:: justify
Alternatively, you can also download the file I generate
by clicking |download_CNT.data|.
.. |download_CNT.data| raw:: html
here
Use of AIREBO potential
-----------------------
.. container:: justify
Then, let us import the LAMMPS data file, and set the
pair coefficients by adding the following lines into *input.lammps*
.. code-block:: lammps
# System definition
read_data cnt_atom.data
pair_coeff * * CH.airebo C
.. container:: justify
Here, there is one single atom type. We impose this type
to be carbon by using the letter *C*.
.. admonition:: Setting AIREBO pair coefficients
:class: info
In the case of multiple atom types, one has to adapt the *pair_coeff* command.
If there are 2 atom types, and both are carbon, it would read: *pair_coeff * * CH.airebo C C*.
If atoms of type 1 are carbon and atoms of type 2 are hydrogen, then *pair_coeff * * CH.airebo C H*.
.. container:: justify
The *CH.airebo* file can be
downloaded by clicking |download_CH.airebo|,
and must be placed within the *breakable-bonds/* folder.
The rest of the *input.lammps* is very similar to the previous one:
.. |download_CH.airebo| raw:: html
here
.. code-block:: lammps
change_box all x final -40 40 y final -40 40 z final -60 60
group carbon_atoms type 1
variable carbon_xcm equal -1*xcm(carbon_atoms,x)
variable carbon_ycm equal -1*xcm(carbon_atoms,y)
variable carbon_zcm equal -1*xcm(carbon_atoms,z)
displace_atoms carbon_atoms move ${carbon_xcm} ${carbon_ycm} ${carbon_zcm}
variable zmax equal bound(carbon_atoms,zmax)-0.5
variable zmin equal bound(carbon_atoms,zmin)+0.5
region rtop block INF INF INF INF ${zmax} INF
region rbot block INF INF INF INF INF ${zmin}
region rmid block INF INF INF INF ${zmin} ${zmax}
group carbon_top region rtop
group carbon_bot region rbot
group carbon_mid region rmid
variable zmax_del equal ${zmax}-2
variable zmin_del equal ${zmin}+2
region rdel block INF INF INF INF ${zmin_del} ${zmax_del}
group rdel region rdel
delete_atoms random fraction 0.02 no rdel NULL 482793
reset_atoms id sort yes
velocity carbon_mid create ${T} 48455 mom yes rot yes
fix mynve all nve
compute Tmid carbon_mid temp
fix myber carbon_mid temp/berendsen ${T} ${T} 0.1
fix_modify myber temp Tmid
.. container:: justify
Note that a large distance of 120 Ångstroms was used for the box size along
the *z* axis, to allow for larger deformation. In addition, the *change_box* command
was placed before the *displace_atoms* to avoid having the CNT crossing the edge of the box.
Start the simulation
--------------------
.. container:: justify
Here, let us impose a constant velocity deformation using the atoms of one
edge while maintaining the other edge fixed (note that for the unbreakable
CNT, the motion was imposed on the 2 edges).
.. container:: justify
First, as an equilibration step, let us set the velocity to 0
for the atoms of both edges. Let us fully constrain the edges.
Add the following lines into LAMMPS:
.. code-block:: lammps
fix mysf1 carbon_bot setforce 0 0 0
fix mysf2 carbon_top setforce 0 0 0
velocity carbon_bot set 0 0 0
velocity carbon_top set 0 0 0
variable L equal xcm(carbon_top,z)-xcm(carbon_bot,z)
fix at1 all ave/time 10 10 100 v_L file output_cnt_length.dat
fix at2 all ave/time 10 10 100 f_mysf1[1] f_mysf2[1] &
file output_edge_force.dat
dump mydmp all atom 1000 dump.lammpstrj
thermo 100
thermo_modify temp Tmid
timestep 0.0005
run 5000
.. container:: justify
Note the relatively small timestep of :math:`0.0005\,\text{ps}` used. A
reactive force field usually requires a smaller timestep than a classical one.
When running *input.lammps* with LAMMPS, you can see that the
temperature deviates from the target temperature of :math:`300\,\text{K}`
at the start of the equilibration, but that
after a few steps, it reaches the target value:
.. code-block:: bw
Step Temp E_pair E_mol TotEng Press
0 300 -5084.7276 0 -5058.3973 -1515.7017
100 237.49462 -5075.4114 0 -5054.5671 -155.05545
200 238.86589 -5071.9168 0 -5050.9521 -498.15029
300 220.04074 -5067.1113 0 -5047.7989 -1514.8516
400 269.23434 -5069.6565 0 -5046.0264 -174.31158
500 274.92241 -5068.5989 0 -5044.4696 -381.28758
600 261.91841 -5065.985 0 -5042.9971 -1507.5577
700 288.47709 -5067.7301 0 -5042.4111 -312.16669
800 289.85177 -5066.5482 0 -5041.1086 -259.84893
900 279.34891 -5065.0216 0 -5040.5038 -1390.8508
1000 312.27343 -5067.6245 0 -5040.217 -465.74352
(...)
Launch the deformation
----------------------
.. container:: justify
After equilibration, let us set the velocity to 15 m/s and run for
a longer duration than previously. Add the following lines into
*input.lammps*:
.. code-block:: lammps
# 0.15 A/ps = 15 m/s
velocity carbon_top set 0 0 0.15
run 280000
.. container:: justify
The CNT should break around step 250000. If not,
use a longer run.
.. container:: justify
By opening the *lammpstrj* file using VMD, it is possible to observe the
bonds breaking at approximately two-thirds of the simulation. If the bonds
do not break, use a longer run.
.. |video_lammps_cnt| raw:: html
this video
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/deformed-dark.png
:alt: carbon nanotube with broken bonds after simulation with LAMMPS and AIREBO
:class: only-dark
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/deformed-light.png
:alt: carbon nanotube with broken bonds after simulation with LAMMPS and AIREBO
:class: only-light
.. container:: figurelegend
Figure: CNT with broken bonds. See the corresponding |breakable_cnt_video|.
.. |breakable_cnt_video| raw:: html
video
.. admonition:: About bonds in VMD
:class: info
Note that VMD guesses bonds based on the distances
between atoms, and not based on the presence of actual
bonds between atoms in the LAMMPS simulation. Therefore the bonds seen
in VMD when using the *DynamicBonds* representation can be misleading.
.. container:: justify
Looking at the evolution of energy again, one can see that the energy is increasing
with the deformation. When bonds break, the energy relaxes abruptly, as can be seen
near $t=210~\text{ps}$ when plotting the evolution of the total energy with time.
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/energy-breakable-dm.png
:alt: energy of the CNT with time - lammps molecular dynamics
:class: only-dark
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/energy-breakable.png
:alt: energy of the CNT with time - lammps molecular dynamics
:class: only-light
.. container:: figurelegend
Figure: Evolution of the total energy of the system with time.
.. include:: ../../non-tutorials/accessfile.rst
.. container:: justify
There is a follow-up to this CNT tutorial as :ref:`mda-label`,
where a post-mortem analysis is performed using Python.
Going further with exercises
============================
.. include:: ../../non-tutorials/link-to-solutions.rst
Plot the strain-stress curves
-----------------------------
.. container:: justify
Adapt the current scripts and extract the strain-stress curves for
the two breakable and unbreakable CNTs:
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/stress-strain-curve-dark.png
:alt: strain strain curve of the CNTs
:class: only-dark
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/stress-strain-curve-light.png
:alt: strain strain curve of the CNTs
:class: only-light
.. container:: figurelegend
Figure: Strain-stain curves for the two CNTs, breakable and unbreakable.
Solve the flying ice cube artifact
----------------------------------
.. container:: justify
The flying ice cube effect is one of the most famous artifacts of
molecular simulations :cite:`wong2016good`.
Download this seemingly simple |input_flying_cube|, which is a simplified
version of the input from the first part of the tutorial.
Run the input with this |data_flying_cube| file
and this |parm_flying_cube| file.
.. |input_flying_cube| raw:: html
input
.. |data_flying_cube| raw:: html
data
.. |parm_flying_cube| raw:: html
parameter
.. container:: justify
When you run this simulation using LAMMPS, you should see that the temperature is
very close to :math:`300\,\text{K}`, as expected.
.. code-block:: bash
Step Temp E_pair E_mol TotEng Press
0 327.4142 589.20707 1980.6012 3242.2444 60.344754
1000 300.00184 588.90015 1980.9013 3185.9386 51.695282
(...)
.. container:: justify
However, if you look at the system using VMD, the atoms are not moving.
.. container:: justify
Can you identify the origin of the issue, and fix the input?
Insert gas in the carbon nanotube
---------------------------------
.. container:: justify
Modify the input from the unbreakable CNT, and add atoms of argon
within the CNT.
.. container:: justify
Use the following *pair_coeff* for the argon,
and a mass of *39.948*:
.. code-block:: lammps
pair_coeff 2 2 0.232 3.3952
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/CNT-gas-dark.png
:alt: CNT with Argon modeled in LAMMPS
:class: only-dark
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/CNT-gas-light.png
:alt: CNT with Argon modeled in LAMMPS
:class: only-light
.. container:: figurelegend
Figure: Argon atoms in a CNT. See the corresponding |gas_cnt_video|.
.. |gas_cnt_video| raw:: html
video
Make a membrane of CNTs
-----------------------
.. container:: justify
Replicate the CNT along the *x*
and *y* direction, and equilibrate the system to
create an infinite membrane made of multiple CNTs.
.. container:: justify
Apply a shear deformation along *xy*.
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/membrane-dark.png
:alt: deformed membrane of CNTs
:class: only-dark
.. figure:: ../figures/level1/breaking-a-carbon-nanotube/membrane-light.png
:alt: deformed membrane of CNTs
:class: only-light
.. container:: figurelegend
Figure: Multiple carbon nanotubes forming a membrane.
.. admonition:: Hint
:class: info
The box must be converted to triclinic to support deformation
along *xy*.