Energy not constant for amorphous SiO2 film in NVE simulation

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I am simulating amorphous SiO2 film in NVE. But The energy increase steadily. The SiO2 film is 606025 (xyz) angstroms. It is periodic in x and y direction and free in z direction. I used tersoff potential, and the time step is 0.25fs. The amorphous SiO2 is prepared by cooling it down from 6000K with periodic direction in three directions. Then I increased the z dimension slowly and got a film. I first equilibrated the film at 300K by rescaling the velocity for 2.5ns, and equilibrated it in NVE for another 2.5ns. But the total energy in NVE simulation increased steadily by 15eV for 2.5 ns. The temperature also increased. Did anyone have similar experience? Thank you.

The script to prepare the film is as following

units metal
dimension 3
boundary p p p
atom_style atomic

processors 4 4 2

read_data SiO2_amorphous.dat

group oxygen type 1
group silicon type 2
pair_style tersoff
pair_coeff * * SiOC.tersoff O Si

timestep 0.00025

velocity all create 300.0 5418580 mom yes rot yes dist gaussian units box

reset_timestep 0

#equilibrium dimensions calculated using NPT
displace_box all x final 0.666563 61.0694 y final -0.333437 60.0694 z final -14.5882 13.514013 units box

thermo_style custom step atoms temp press pe ke etotal vol
thermo 10000

dump 1 all xyz 1000000 dump.coord.SiO2_mini.*
restart 1000000 restart.sio2.mini.*
min_modify line quadratic
minimize 0.0 1.0e-9 10000 10000

fix NVT all nvt temp 300.0 300.0 0.025
run 1000000
thermo_style custom step atoms temp press pe ke etotal zlo zhi vol

#open up z direction periodic boundary
fix deformZ all deform 500000 z delta -10.0 10.0 units box
run 20000000

Regards,

Liang Chen

2

A 0.25 fmsec timestep is normally overkill for Tersoff.
How much did the Temp increase? Maybe Aidan
has a comment on this.

Steve

3

Thank you, Steve. The increase of the temperature is corresponding to the total energy increase, and it is about 10K for my system.

Liang Chen

4

Energy drift is affected by many things, including the presence of stiff
regions in the potential. It is possible that the parameters you are using
(SiOC.tersoff is not a LAMMPS file, so I can't say what it is) are
unusually stiff. You should also verify that the parameters are correct,
meaning that they accurately reproduce results published by the authors of
that potential. You should also examine the dependence of energy
conservation on timestep size for a perfect crystal.

5

Thank you for your suggestion, Aidan. The tersoff potential SiOC.tersoff is same as the Lammps file SiO.tersoff. I copied the Si and O parameters from the file (Mumetoh 2007).

I also tested the perfect crystal. The energy drift was much smaller. It is about 0.5eV for 2.5ns.
If I first equilibrated the SiO2 crystal using NPT, and then there was no energy drift in the following NVE simulation. But this did not work for the amorphous SiO2 film. Thanks a lot.

Liang Chen

6

A 10 degree drift over 2.5 nanosec (10 million timesteps)
is not a lot. All MD drifts over long timescales. It's just
a 2nd-order integrator.

Steve