Thin Slab Model Got Twisted During Relaxation

LAMMPS LAMMPS General Discussion
1

Dear LAMMPS Users,

I was trying to construct a thin slab SiC model to calculate its potential energy at 300 K. PBC is applied to the system, and the model is smaller than the box. After energy minimization, I used NVT and NPT ensembles to the system, but during these two relaxations, the shape of the model got changed. It actually got twisted a little bit. I would like to ask if the potential energy is still the correct value. Can I control the shape of the model by immobilizing the atoms in the corners of the slab? I’m using Sep 21 stable version LAMMPS. My code is attached here. Thank you very much!

Simulation of SiC

Initialization

echo both
units metal
dimension 3
boundary p p p
processors 4 4 4
atom_style atomic

Definition of the system box and creation of the atoms

lattice custom 4.3596 &
a1 1.0 0.0 0.0 &
a2 0.0 1.0 0.0 &
a3 0.0 0.0 1.0 &
basis 0.0 0.0 0.0 &
basis 0.5 0.0 0.5 &
basis 0.0 0.5 0.5 &
basis 0.5 0.5 0.0 &
basis 0.25 0.25 0.25 &
basis 0.75 0.25 0.75 &
basis 0.25 0.75 0.75 &
basis 0.75 0.75 0.25

region box block 0.0 360 0.0 200 0.0 200 units box
create_box 2 box

region sic block 60 300 25 175 100 114 units box

create_atoms 2 region sic &
basis 1 1 &
basis 2 1 &
basis 3 1 &
basis 4 1 &
basis 5 2 &
basis 6 2 &
basis 7 2 &
basis 8 2

group SiAtom type 1
group CAtom type 2

neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes

mass 1 28.08
mass 2 12.01

write_data sic.xyz

Potentials settings

pair_style tersoff
pair_coeff * * SiC.tersoff Si C

Energy minimization

dump 1 all atom 500 minimize.atom
thermo_style custom step temp pe ke etotal press vol
thermo 50
min_style cg
minimize 1.0e-8 1.0e-8 10000 20000
undump 1

timestep 0.001
thermo 100
thermo_style custom step temp pe ke etotal press vol

variable steps equal step
variable E equal pe

velocity all create 300 123456 mom yes rot yes dist gaussian

dump 2 all custom 1000 relax1.nvt id type mass x y z
fix 1 all nvt temp 300 300 0.1
run 10000
undump 2
unfix 1

dump 3 all custom 1000 relax2.npt id type mass x y z
fix 2 all npt temp 300 300 0.1 iso 0 0 1
fix 3 all print 1000 “${steps} $E” file energy.txt screen no
run 70000
undump 3
unfix 2
unfix 3

2

There is no such thing. When your system is a finite temperature and in equilibrium then your potential energy will be a distribution, aka as an ensemble. Please pick up a statistical mechanics text book to learn more about ensembles.

No you didn’t. Apart from the incorrect grammar you are using the term “ensemble” incorrectly (like too many people here). You were using fix nvt and fix npt with the former applying velocity verlet time integration and a Nose-Hoover thermostat and the latter in addition a Nose-Hoover barostat. Neither will represent the NVT or NPT ensemble because your system is not a bulk system.

Correct with respect to what. Of course LAMMPS computes the correct potential energy for the model and geometry you input. If that is not the intended geometry or model is not for LAMMPS to decide.

How do you know that undulations are not the expected behavior for your geometry?

3

Thanks for your reply, sir. I am going to learn about it and revise my simulation based on what you mentioned. Thank you.