INVESTIGATE TEMPERATURE AND DENSITY DISTRIBUTION OF ARGON CONFINED IN SILICON WALLS BY LAMMPS

I.     Simulation and Input Script

Similar to building cubic box containing liquid argon, now similarly we use LAMMPS to build a rectangular box which contains both argon and silicon. In this box, silicon is built at both sides to create walls and liquid argon is placed at the middle. The walls are called hot wall and cold wall respectively in the z direction due to the fact that the higher temperature will be applied in left wall and lower temperature will be applied in the right wall to investigate the heat transfer.

 Basically, we can choose the simulation box size depends on our purposes of simulation. However, one important thing is that we should choose the sizes (width, height, and length) such that they are the positive numbers multiples with the lattice constant of the walls. This is easier for us to control the number of atoms in the box. In this example, we use the sizes as follows:

-          x direction:  5 x a

-          y direction:  5 x a

-          z direction:  24 x a in which the length for the wall at each side is 5 x a

where a is the lattice of constant of silicon

Most of researchers in many papers usually place the width at x direction, the height at y direction, and the length as z direction. Therefore, we should follow this as a convention.  The details about the simulation box dimension will be described in the script explanation section.

For the simulation process, in this kind of heat transfer problem, we have two separate steps called NVT ensemble and NVE ensemble, respectively. This is different from the cubic box example which only has the NVT step.

-          NVT ensemble: NVT means constant number of atom number, constant of volume, and constant of temperature. This step is aimed to created a equilibrium system at room temperature (300K) before we do the temperature difference at both sides to investigate the heat transfer. In this step, the Maxwell-Boltzmann velocity distribution is used as the initial velocity for all molecules and Nose-Hoover thermostat is used to maintain the system temperature at 300K.

-          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 wall and one cold wall) with the method of energy conservation. This step is executed right after the NVT is completed.

In NVE step, when the system approach the steady state, we can start to calculate the temperature values of atoms as well as the atom densities. We can divide the rectangular to many smaller regions call slab bins. We call the system having coarse bins if the number of bins is less than 100 bins, and fin bins if the number of bins is greater than 200 bins.

You can refer to the following figures for the simulation box schema and process chart.

1.    LAMMPS Input script example

In addition to the script input below, you'll need Si.sw file for the silicon potential. (put them together with input script)

You can google it "si.sw potential file" or go to the link below to get Si.sw

https://github.com/lammps/lammps/blob/master/potentials/Si.sw

a.      NVT: ----------------------------------------------------------------------------------------------------------------------------------

units          metal

dimension 3

boundary p p f 

atom_style     atomic

lattice  diamond 5.4307

region   sim_domain block 0 27.15 0 27.15 -34.58 99.76 units box

region   lohalf block 0 27.15 0 27.15 -32.58 0 units box

region   hihalf block 0 27.15 0 27.15 65.18 97.76 units box

region   midfluid block 0 27.15 0 27.15 6.3958 59.7758 units box

region   hol1 block 0 27.15 0 27.15 -33 -31 units box

region   hol2 block 0 27.15 0 27.15 -31 -27 units box

region   col1 block 0 27.15 0 27.15 97 99 units box

region   col2 block 0 27.15 0 27.15 93 97 units box 

create_box     2 sim_domain

create_atoms   1 region lohalf

create_atoms   1 region hihalf

lattice        fcc 5.4307

create_atoms   2 region midfluid 

group bksilicon type 1

group bkargon type 2

group ho-wall region lohalf

group co-wall region hihalf

group ho-layer1 region hol1

group co-layer1 region col1

group ho-layer2 region hol2

group co-layer2 region col2

group realhotwall subtract ho-wall ho-layer1

group realcoldwall subtract co-wall co-layer1

group ho-nve subtract ho-wall ho-layer1 ho-layer2

group co-nve subtract co-wall co-layer1 co-layer2

group myall subtract all ho-layer1 co-layer1

mass           1 28.0855

mass           2 39.948002

neighbor 0.1 bin

neigh_modify   delay 10 check no

velocity all create 300.0 3213112 mom yes rot yes dist gaussian

pair_style hybrid sw lj/cut 10.0

pair_coeff * * sw Si.sw Si NULL

pair_coeff 1 2 lj/cut 0.005 2.750

pair_coeff 2 2 lj/cut 0.0103 3.405 

fix 5 myall nvt temp 300.0 300.0 0.01

thermo_style custom step etotal

thermo         10000

dump 1 all atom 10 tmp.dump

dump 2 all xyz 50000 mycawa.xyz

timestep 0.001

restart 2000000 silconargon160_90K_NTC.restart

run            2000000

b.      NVE ----------------------------------------------------------------------------------------------------------------------------------

read_restart silconargon160_90K_NTC.restart.2000000

pair_style hybrid sw lj/cut 10.0

pair_coeff * * sw Si.sw Si NULL

pair_coeff 1 2 lj/cut 0.005 2.750

pair_coeff 2 2 lj/cut 0.0103 3.405

fix 3 myall nve

fix 4 ho-layer2 langevin 160 160 0.01 699483 tally yes

fix 11 co-layer2 langevin 90 90 0.01 890765 tally yes

compute ske bksilicon ke/atom

compute ake bkargon ke/atom

variable temp atom (c_ske+c_ake)/1.5*1.6/1.3806503*10000

compute cc1 all chunk/atom bin/1d z lower 0.005 units reduced

fix 33 all ave/chunk 1 4000000 8000000 cc1 v_temp file all200.profile

compute cc2 bksilicon chunk/atom bin/1d z lower 0.005 units reduced

fix 44 bksilicon ave/chunk 1 4000000 8000000 cc2 v_temp file si200.profile

compute cc3 bkargon chunk/atom bin/1d z lower 0.005 units reduced

fix 55 bkargon ave/chunk 1 4000000 8000000 cc3 v_temp file ar200.profile

compute cc4 all chunk/atom bin/1d z lower 0.005 units reduced

fix 66 all ave/chunk 1 4000000 8000000 cc4 density/number norm sample file denall.txt

compute cc5 bksilicon chunk/atom bin/1d z lower 0.005 units reduced

fix 77 bksilicon ave/chunk 1 4000000 8000000 cc5 density/number norm sample file densi.txt

compute cc6 bkargon chunk/atom bin/1d z lower 0.005 units reduced

fix 88 bkargon ave/chunk 1 4000000 8000000 cc6 density/number norm sample file denar.txt

thermo 10000

dump 4 all xyz 200000 mycawa_NVE.xyz

dump 5 all custom 10000 NVEdensity_*.cfg id type x y z

timestep 0.001

run 6000000

2.    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. Details will be explained later.

I.     Run simulation

The step-by-step guideline for running simulation in work stations or clusters is presented clearly in Chapter III. Only one important notice here is that as mentioned, 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.

II.      Expected Results

After the simulation is finished, based on the output values about temperature, number density, we have to process the values to have temperature distribution and density distribution silicon and argon in simulation box. Followings are the distributions obtained using Excel.

Calculation of Number Density:

Please, read the given link very carefully (https://lammps.sandia.gov/doc/fix_ave_chunk.html)

The density/number value means the number density is computed for each chunk, i.e. number/volume. The output values are in units of 1/volume or density (number_of_atom_in_bin/volume). See the units command doc page for the definition of density for each choice of units, e.g. gram/cm^3. If the chunks defined by the compute chunk/atom command are spatial bins, the volume is the bin volume. Otherwise it is the volume of the entire simulation box.

ρ=(Nσ^3)/V , where, ρ = Number Density, σ = Molecular size, V = Volume

Therefore, to get the number density from LAMMPS output file, number/density column should be multiplied by the cube of molecular size.

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

Now, I think you can navigate your own path over the sea of internet to study MD........ Bon voyage!