Pressure equilibration

LAMMPS LAMMPS Mailing List Mirror
1

Dear Lammps users,

I am trying to equilibrate my system at 300k and 0 bar, I suffered from a high pressure at start(~40000bar) so I used nve/lim and made box relaxation and eventually the box pressure started to drop massively, I then tried to equilibrate for a very long time and using the recommended Tdamp and Pdamp ( my time step is 0.001, Tdamp=0.1 and Pdamp=1) but the pressure kept fluctuating from -10000 bar to 10000 bar. Can you advice please what may be the reason for that?

Thanks in advance

Regards

units metal
atom_style charge
boundary p p p

lattice custom 1 a1 6.607 0 0 a2 0 6.607 0 a3 0 0 5.982 basis 0.0 0.75 0.125 basis 0.0 0.25 0.875 basis 0.5 0.75 0.375 basis 0.5 0.25 0.625 basis 0.0 0.75 0.625 basis 0.0 0.25 0.375 basis 0.5 0.75 0.875 basis 0.5 0.25 0.125 basis 0.0 0.0658 0.1955 basis 0.0 0.9342 0.8045 basis 0.5 0.0658 0.3045 basis 0.5 0.9342 0.6955 basis 0.0 0.5658 0.8045 basis 0.0 0.4342 0.1955 basis 0.5 0.5658 0.6955 basis 0.5 0.4342 0.3045 basis 0.6842 0.25 0.9455 basis 0.8158 0.75 0.4455 basis 0.8158 0.25 0.5545 basis 0.6842 0.75 0.0545 basis 0.1842 0.25 0.5545 basis 0.3158 0.75 0.0545 basis 0.3158 0.25 0.9455 basis 0.1842 0.75 0.4455

region R1 block -1 1 -1 1 -1 1

create_box 3 R1

create_atoms 3 region R1 basis 1 1 basis 2 1 basis 3 1 basis 4 1 basis 5 2 basis 6 2 basis 7 2 basis 8 2 basis 9 3 basis 10 3 basis 11 3 basis 12 3 basis 13 3 basis 14 3 basis 15 3 basis 16 3 basis 17 3 basis 18 3 basis 19 3 #basis 20 3 basis 21 3 basis 22 3 basis 23 3 basis 24 3

mass 1 91.224
mass 2 28
mass 3 16
group 1 type 1 # zircon
group 2 type 2 # silicon
group 3 type 3 # oxygen

set group 1 type 1 charge 3.8
set group 2 type 2 charge 2.0
set group 3 type 3 charge -1.45

velocity all create 0 7928459 dist gaussian

dump atom all custom 1000 dump.coordinates id type x y z

System interatomic potential

kspace_style ewald 5e-3

pair_style born/coul/long 10

pair_coeff 1 3 1967.0 0.305004 0.0 0.0 0.0
pair_coeff 2 3 1277.0 0.227225 0.0 0.0 0.0
pair_coeff 3 3 1755.0 0.306820 0.0 0.0 0.0
pair_coeff 1 2 0.0 0.306820 0.0 0.0 0.0
pair_coeff 2 2 0.0 0.306820 0.0 0.0 0.0
pair_coeff 1 1 0.0 0.306820 0.0 0.0 0.0

timestep 0.001
thermo 100

fix 1 all nve
run 10000
unfix 1

fix 1 all nve/limit 1
run 10000
unfix 1

fix 2 all nve/limit 0.8
run 10000

unfix 2

fix 3 all nve/limit 0.6
run 10000
unfix 3

fix 4 all nve/limit 0.4
run 10000
unfix 4

fix 1 all nve/limit 0.3
run 10000
unfix 1

fix 1 all nve/limit 0.2
run 10000
unfix 1

fix 1 all nve/limit 0.1
run 10000
unfix 1

fix 1 all nve/limit 0.07
run 10000
unfix 1

fix 1 all nve/limit 0.06
run 10000
unfix 1

fix 1 all nve/limit 0.05
run 10000
unfix 1

fix 1 all nve/limit 0.04
run 10000
unfix 1

fix 1 all nve/limit 0.03
run 10000
unfix 1

fix 1 all nve/limit 0.02
run 10000
unfix 1

fix 1 all nve/limit 0.01
run 10000
unfix 1

#BOX RELAXATION

fix 1 all box/relax iso 0.0

#minimization

min_style cg

minimize 1e-25 1e-25 5000 100000
run 10000

unfix 1

#Equilibration
fix 1 all npt temp 300.0 300.0 0.01 iso 0.0 0.0 0.1
run 5000000

2

Dear Lammps users,

I am trying to equilibrate my system at 300k and 0 bar, I suffered from a
high pressure at start(~40000bar) so I used nve/lim and made box relaxation
and eventually the box pressure started to drop massively, I then tried to
equilibrate for a very long time and using the recommended Tdamp and Pdamp
( my time step is 0.001, Tdamp=0.1 and Pdamp=1) but the pressure kept
fluctuating from -10000 bar to 10000 bar. Can you advice please what may be
the reason for that?

​please follow the advice given in the mailing list guidelines and *first*
research the mailing list archives. there are a gazillion of discussion
about pressure fluctuations.

...and then you need to simplify your input. a lot of it doesn't make much
sense. e.g. you go through this elaborate sequence of ultra short runs with
nve/limit, where you shorten the limit with each set. the opposite order
would make much more sense. but also, you do a minimization *after* that
sequence. common sense dictates that this should be reversed, i.e. what is
the point of going through an elaborate scheme to avoid having atoms move
too fast while taking a system from high to low potential energy, when the
very purpose of a minimization is the same thing, and without all those
complications. ​if you start from a very bad initial configuration, you may
need to do a second minimization if running the MD picks up too much energy
quickly. but overall, it would probably be much more effective to start
with a more reasonable initial configuration or a more suitable force
field. if you have to remove much potential energy and have an unstable
system, then the problem is usually not how your get it equilibrated but
rather why the structure is unstable, which may be caused by a bad geometry
or bad force field settings.

axel.