Orientated quartz with BKS potential

Hi, all

I am using LAMMPS to simulate an orientated alpha-quartz with BKS potential, but I am facing a problem that the constructed structure built by LAMMPS is bombed at the beginning of simulation.

I am very confused that when the orientation is normal, a-axis along x direction, and c-axis along z direction, the results are consistent with other papers’ results. But when I try to do some different orientations, the different situation comes, and I am very grateful if someone could offer some information.

Many thanks.

My script is as follows:

------------------------ INITIALIZATION ----------------------------

units metal
dimension 3
boundary p p p
atom_style charge

----------------------- ATOM DEFINITION ----------------------------

lattice custom5.45 &
a1 0.9095 0.0000 0.0000 &
a2 -0.4547 0.7876 0.0000 &
a3 0.0000 0.0000 1.0000 &
basis 0.4697 0.0000 0.0000 basis 0.0000 0.4697 0.6667 basis 0.5303 0.5303 0.3333 basis 0.4135 0.2669 0.1191 basis 0.2669 0.4135 0.5475 basis 0.7331 0.1466 0.7857 basis 0.5865 0.8534 0.2142 basis 0.8534 0.5865 0.4524 basis 0.1466 0.7331 0.8809 &
orient x 4 5 1 &
orient y 20 -17 5 &
orient z 1 0 -4 &
origin 0.0 0.15 0.31region whole block 0 20 0 20 0 20 units lattice
create_box 2 whole
create_atoms 1 box &
basis 1 1 &
basis 2 1 &
basis 3 1 &
basis 4 2 &
basis 5 2 &
basis 6 2 &
basis 7 2 &
basis 8 2 &
basis 9 2
mass 1 28.0855
mass 2 15.9994
group siliconatoms type 1
group oxygenatoms type 2
set group siliconatoms charge 2.4
set group oxygenatoms charge -1.2

------------------------- SETTINGS ---------------------------------

pair_style buck/coul/long 10.0 10.0
pair_coeff 1 1 0.00000000 0.1000000000 0.000000
pair_coeff 1 2 18003.7572 0.2052048149 133.5381
pair_coeff 2 2 1388.77300 0.3623188406 175.0000
timestep 0.0005

kspace_style ewald/disp 1.0e-4

neighbor 4 bin
neigh_modify delay 0 every 1 check yes page 1000000 one 100000
thermo 1
thermo_style custom step temp vol lx ly lz press pxx pyy pzz pe v_density

velocity all create 300 4928459 dist gaussian
fix 1 all npt temp 300 300 0.05 iso 0 0 0.5 drag 1

run 2000

The error is as follows:

On timestep 1 in your output the temperature went from
300 to 76,000 - obviously the dynamics are bad.

You probably created atoms that are nearly overlapping
across a periodic boundary. Read the lattice command
doc page carefully, experiment with a small system,
visualize the results, dump the coords and verify
that you are getting the system you expect.

Steve