Molecular dynamics simulation of Argon gas

Input script

#-------------------------------------------------------------------------------------------------------------------------

units metal #What are the units you will use to specify various things in the input file?
boundary p p p #Specify periodic boundary condition are needed in all three faces
atom_style atomic #What style of atoms is to be used in the simulation
log logfile.txt #Write the log file to this text file. All thermodynamic information applicable to the entire system
#================================================
region forbox block 0 45.8 0 45.8 0 45.8 units box
#Refers to an abstract geometric region of space. units box refers to the fact that the size of the box is specified in the units as given in the units command. The name "forbox" refers to the region ID so that you can refer to it somewhere else in this input script.
create_box 1 forbox #Create the box

#=======================================================================
lattice fcc 4.58 # Since we have given fcc as lattice type no need to mention basis for this
#=======================================================================
create_atoms 1 region forbox  basis 1 1 basis 2 1 basis 3 1 basis 4 1 units box
mass 1 39.948 #Mass of atom type 1 is 39.48 [mass units grams/mole]
pair_style lj/cut 10
# k_B = 8.6173303e-5 eV/K #How are atoms interacting. Provide the name of the potential and the cooresponding cut-off distance
pair_coeff 1 1 0.01006418 3.3952 #The coefficient of the lj potential for the interactions of atom type 1 with 1
#=======================================================================
group ar type 1 #Group all the argon types (argon type is of type 1). All atoms of type 1 are in group with the name 'ar'
#=======================================================================
dump dump_1 all custom 1 dump_initial_config.dump id type x y z ix iy iz vx vy vz
#=========================================================
run 1
undump dump_1 # Stop dumping to this file
#=======================================================
minimize 1e-25 1e-19 10000 10000 #<Minimize the energy using a conjugate gradient step.
print "Finished Minimizing"
#minimize etol ftol maxiter for minimizer maxiter for force_energy valuation
variable ener equal pe
#=========================================================
timestep  0.001
velocity all create 300 102939 dist gaussian mom yes rot yes
# Set the velocities of all the atoms so that the temperature of the system
# is 300K. Make the distribution Gaussian.
fix 1 all nve #this is equlibration process.
dump dump_1 all custom 100 trj.dump id type x y z ix iy iz vx vy vz #<Dump all the atoms to the file dump.min>
variable poten equal pe-${ener}
thermo_style  custom step time temp pe ke etotal press vol v_poten
#What to print in the logfile.txt?
#=========================================================
thermo 100 #How frequently to print the thermodynamic information#
run 10000 # run with active settings as many runs as required. timestep*No. of. steps =10ps
undump dump_1 # Stop dumping information to the dump file.
unfix 1
# Unfix the NVE. Additional lines if any will assume that this fix is off.
#End
#------------------------------------------------------------------------------------------------------------------------

Visit LAMMPS Manual for a detailed explanation of various commands.

Visualisation

Explanation 

1. This is a typical example of molecular dynamics simulation. When you observe the atoms are vibrating. The reason for this is, I have performed an equilibration process, which I haven't done in the previous blog.

2. I have plotted this temperature graph from the log file generated (log.txt). If you notice, I have given the initial temperature as 300K before equilibration process. But as you can see, the temperature value drops to 150K and gets constant. There are mathematical reason behind this. If you are interested, contact me so that I can share the calculations. But as for now what you can remember is whatever the temperature you give, after equilibration the system will tend to maintain half the given temperature.

3. This temperature plot can also be used to identify whether or not your system has reached equilibrium. As you can see the plot, there is a minor oscillation of temperature around 150K. If you run the simulation for a greater time step (increase the run command argument value) the oscillation will further reduce, and thus the system will tend to reach equilibrium. Thus, if the temperature becomes constant at a value half the initially set value, with minimum oscillation, then we can say that the system is in equilibrium. If you find more oscillation in temperature, then run your simulation for an increased number of steps so as to reach equilibrium.