8. LAMMPS-Molecular Dynamics Tutorial

--- Quick Link ---

This tutorial is for beginners who never used LAMMPS ( https://lammps.sandia.gov/ ) before.

A. Investigation of solidification/crystallization of Argon in nano-scale by using Molecular Dynamics simulation with LAMMPS

-----------------------------------------

CONTENT

I. Modeling and simulation script

1.1. Modeling

1.2. Simulation script

1.3. Explaination of script’s function

II. Run simulation

III. Show results

3.1. Data files

3.2. Virtual Molecular Dynamics (VMD)

3.3. Simulation results

------------------------------------------

I. Modeling and simulation script

1.1 Modeling

In this case, we build the virtual box with defined size around ten nanometer, which is then filled with Argon atoms. After setting the chosen bulk temperature and density, we investigate the local temperature and density distribution of Argon in the box at equilibrium.

Fig.1. Scheme of our simulation process

To observing solidification of Ar as well as the crystallization, we produce the solidification progress:

Fig.2. Solidification progress

1.1 Simulation script

First, we create an input.txt file which contains the script series below:

# ---------- Initialize Simulation ---------------------

units metal

dimension 3

boundary p p p

atom_style atomic

lattice fcc 3.65

#-------------------Define Region ---------------------------------------

region box1 block 0 57.4 0 57.4 0 57.4 units box

region box2 block 0 36.5 0 36.5 0 36.5 units box

#-------------------Create Atoms ---------------------------------------

create_box 1 box1

create_atoms 1 region box2

# ---------- Define Interatomic Potential ---------------------

mass 1 39.948002

neighbor 0.1 bin

neigh_modify delay 10 check no

velocity all create 150.0 3213112 mom yes rot yes dist gaussian

pair_style hybrid lj/cut 10.0

pair_coeff 1 1 lj/cut 0.0103 3.405

fix 1 all nvt temp 150.0 10.0 0.01

# ---------- Define Settings ---------------------

compute 2 all chunk/atom bin/1d z lower 0.05 units reduced

# ---------- Run Minimization ---------------------

fix 33 all ave/chunk 1 2000000 2000000 2 v_temp file temp20bin.profile

thermo 100000

dump 2 all xyz 200 argon_150_10_model.xyz

dump 3 all custom 10000 Argon_HT_*.cfg id type x y z

timestep 0.001

restart 6000000 abcdetemp100.restart

run 2000000

1.1 Explanation of scripts’ function

Explanations for the scripts in input.txt file respectively are shown as in the table below:

*The table utilized some information from LAMMPS website http://lammps.sandia.gov

** There are couple ways to generate the bulk density of Argon in box, e.g. change group volume, or change lattice space distant (just only use this way when the objective materials is in amorphous phase that have no relationship with initial set up lattice space constant)

Actually, every single LAMMPS command can be used in various different circumstances and purposes. In order to have more understanding about these codes in or out of this study, we can access the LAMMPS-command section on LAMMPS homepage http://lammps.sandia.gov

II. Run simulation

The input.txt file is used for running MD simulation with LAMMPS software. For installing LAMMPS, we access to LAMMPS download homepage:

http://lammps.sandia.gov/download.html

In our case, we run simulation on a work station with Microsoft Window platform, then we choose window package:

http://rpm.lammps.org/windows.html , download the latest version and install to the work station or your CPU.

(above link has expired. use below.)

https://lammps.sandia.gov/doc/Install_windows.html

Fig.3. LAMMP download website

Fig.4. Installing LAMMPS

When the installing complete, we install MPI based parallelism downloaded from the above link.

After saving the input script, to run LAMMPS simulations open a terminal window as outlined above and change to the directory with your input script and then run the command:

mpiexec -localonly N lmp_mpi.exe < inputfile

N is number of your CPU processor prepared to use

Inputfile is the input file name; in our example case is input.txt

The parallel executable can also run on a single processor by typing something like:

lmp_mpi -in inputfile

And the parallel executable also includes OpenMP multi-threading, which can be combined with MPI using something like:

mpiexec -localonly 2 lmp_mpi -in inputfile -pk omp 2 -sf omp

For running the non-MPI executable, follow these steps:

Get a command prompt by going to Start->Run... , then typing “cmd”.
Move to the directory where you have your input script
At the command prompt, type “lmp_serial -in inputfile”

The serial executable includes support for multi-threading parallelization from the styles in the USER-OMP packages.

To run with, e.g. 4 threads, type “lmp_serial -in inputfile -pk omp 4 -sf omp”

III. Show results

3.1 Data files

As a completed results, we will have .profile files which contains information about N (numbers of atom), local temperature at each bins that we can calculate the density distribution also. Log file can help you to cover again the running process. Microsoft Excel is the consistent choice for opening and doing analysis these files mentioned above. And, .xyz is important file, as the Visual Molecular Dynamics (VMD) input file, which will show us the transient movement of the argon system visually.

Fig.5. Initial input script file before running (200.txt is an example of input.txt)

Fig.6. Data files in running and after completed (the important result file is .xyz file)

3.2 Visual Molecular Dynamics

To open the data results as in .xyz file, we use Visual Molecular Dynamics (VMD) software. Download VMD install package from the link below:

https://www.ks.uiuc.edu/Development/Download/download.cgi?PackageName=VMD

if the link above has expired, google "VMD visual molecular dynamics"

Register your VMD account and login for accessing download. Open the latest updated install file by administration right, follow the installing structure.

After finish the installing, let’s step by step use VMD to show own results:

Open VMD

VMD Main > File > New Molecule > Browse your .xyz file

>Determine file type> automatically

>Speed> Low

>Display> orthographic

>Graphic> Representations> Drawing Method>VDW

>Coloring Method >(upon you) e.g Res Type

>Graphic>Color>Display> Background>White color

>File > New Molecule > Browse your .xyz file > Load

Fig.7. Molecule browser sub-window of VMD Main-window

Fig.8. VMD Main window

Then, we will see the result.

3.3 Simulation result (Ar at vapor state (150K) and solid state (10K))

Fig.9. Results of Ar solidification process

As an analysis, vapor state of gas phase of Ar has shown the random dissipation of group of Ar atom. After cooling down to 10^oK, the atoms link tightly together to show off the ordering structure of FCC.

The video is also presented on YouTube:

3D: https://www.youtube.com/watch?v=MUfjPnns9mw&feature=youtu.be

2D: https://www.youtube.com/watch?v=gvDPdCJ05no&feature=youtu.be

B. Investigate Temperature Distribution of Argon in LAMMPS

1.1 Modeling

We have to build a rectangular box which contains argon. In this box, the hot layer (130K) and cold layer (90K) is built respectively in the z direction at the outmost position of the box as shown in the above figure. So, there will be heat transfer across the simulation domain.

For the simulation process, we have to done the simulation in two separate input file or steps. These steps are NVT ensemble and NVE ensemble.

· NVT ensemble: NVT means constant number of atom number, constant of volume, and constant of temperature. This step help to create an equilibrium state at 100K temperature.

· NVE ensemble: NVE means constant number of atom number, constant of volume, and constant of energy. In this step, different temperature will be applied to both sides (one hot and one cold layer).

additional comments to think about)

First step is NVT for the whole system to reach the Equilibrium State. The second step is not not NVE system and it could be easily misunderstood as NVE since no thermostat will be applied to the middle fluid. However, both hot and cold side will supply and drain heat energy by thermostat. Therefore, the system is more exact to be described as the Steady State in macroscopic view, and combination of different NVT bins (since temperature of each bin is different) for microscopic view.

1.2 Simulation script

First, we create an input.txt file which contains the script series below:

I. NVT (for initialization)

input.txt

# ---------- Initialize Simulation ----------------------------------------------

units metal

dimension 3

boundary p p p

atom_style atomic

lattice fcc 4.8

#-------------------Define Region -----------------------------------------------

region box block 0 36.5 0 36.5 0 73 units box

region mid block 2 33.5 5 33.5 3 69 units box

region col1 block 0 36.5 0 36.5 69 72 units box

region hol1 block 0 36.5 0 36.5 0 3 units box

#-------------------Create Atoms ----------------------------------------------

create_box 1 box

create_atoms 1 region mid

#------------------Divide Groups------------------------------------------

group argon type 1

group ho-layer1 region hol1

group co-layer1 region col1

group myall subtract all ho-layer1 co-layer1

# ---------- Define Interatomic Potential ------------------------------------

mass 1 39.948002

neighbor 0.1 bin

neigh_modify delay 10 check no

velocity all create 100.0 3213112 mom yes rot yes dist gaussian

pair_style lj/cut 10.0

pair_coeff * * 0.0103 3.405

compute ar_temp all temp

# ---------- Run Minimization ---------------------

fix 5 all nvt temp 100.0 100.0 100

thermo_style custom step c_ar_temp etotal

thermo 10000

timestep 0.004

dump d2 all xyz 50000 argon_130_90K_NVT.xyz

restart 50000 argon_130_90K_NVT.restart

run 500000

II. NVE (Thermostat on both ends for heat transfer.)

*make sure you have restart file from first step. filename "argon_130_90K_NVT.restart.500000"

additional comments to think about)

input.txt

# ---------- Initialize Simulation ----------------------------------------------

read_restart argon_130_90K_NVT.restart.500000

#-------------------Define Region -----------------------------------------------

region box block 0 36.5 0 36.5 0 73 units box

region mid block 2 33.5 5 33.5 3 69 units box

region col1 block 0 36.5 0 36.5 69 72 units box

region hol1 block 0 36.5 0 36.5 0 3 units box

#------------------Divide Groups--------------------------------------------------

group ho_layer1 dynamic all region hol1 every 1

group co_layer1 dynamic all region col1 every 1

group myall dynamic all region mid every 1

# ---------- Define Interatomic Potential -------------------------------------

mass 1 39.948002

neighbor 0.1 bin

neigh_modify delay 10 check no

pair_style lj/cut 10.0

pair_coeff * * 0.0103 3.405

# ---------- Define Settings ------------------------------------------------------

compute ke all ke/atom

variable temp atom c_ke/1.5*1.6/1.3806503*10000

compute 2 all chunk/atom bin/1d z lower 0.05 units reduced

# ---------- Run Minimization -------------------------------------------------

fix 33 all ave/chunk 1 200000 1000000 2 v_temp file all200.profile

fix 3 all nve

fix 4 ho_layer1 langevin 130 130 0.01 699483 tally yes

fix 11 co_layer1 langevin 90 90 0.01 890765 tally yes

thermo 10000

dump d3 all xyz 10000 argon_130_90K_NVE.xyz

timestep 0.004

run 500000

1. Script explanation

a. NVT

b. NVE

After the NVT step is finished, we need to run NVE step using the restart file created from NVT step.

II. Run simulation

We have to run the NVT step first to get the restart file. Then the NVE is run separately the input is the restart file generated from NVT step.

III.Result

As a completed results, we will have .profile file which contains information about N (numbers of atom), local temperature at each bins that we can calculate the temperature distribution also. Log file can help you to cover again the running process. Microsoft Excel is the consistent choice for opening and doing analysis these files mentioned above.

Temperature distribution of liquid argon

Figure: Temperature distribution of liquid argon between 130K and 90K

- End.