SiO2 - Core Shell Model - 3 Body terms

Dear lammps Users,

-I am trying to do an energy minimization simulation of SiO2. The model i am using contains a core shell for the Oxygen atoms and a 3 body term (O - Si - O).
-I am using energy minimization features in gulp (initial and final energies) and lattice constants from the literature as references[1].
-I am using constant pressure relaxation.
-I have tracked the problem in the three body documentation on the website and i understood that the code actually reads everything [Si (central atoms) O O (outer atoms) 1.04 (k including the half) 109 (theta) 1.8 (cutoff)] if i have similar outer atoms. I also understood that the units of these arguments depends on the units defined in the potential file.
-I have tried every combination of different implementation using both software in order to figure out where is the problem.

  • I have started by the coordinates from the same input file in the crystallographic data for both software. Also its a triclinic structure with XY tilt only.

  • However only when i implement all the features of the simulation i find a small drift in the total energy which causes a significant change in the final lattice constants

  • I removed the shells and kept only the two body interactions (same initial and final energies and same final lattice constants )

  • I removed the shells implemented 2 and 3 body terms (same initial and final energies and same final lattice constants )

  • I kept the shells and only implemented 2 body terms (a small change in final energies and different final lattice constants when compared to the gulp exact implementation )

  • I implemented everything (a small change in final energies and different final lattice constants when compared to the gulp exact implementation and literature)

  • Is there a way or is it even necessary, in the three body implementation, to implement the cutoff distance between the outer atoms? In this case O and O?
    -I concluded that there is a problem with the core shell model implementation but i cant recognize what is the problem and i have tried to change everything and monitor what is affected but couldn’t recognize what is the problem.

data.minl.nff.SiO2

SiO2

15 atoms
6 bonds
0 angles
0 dihedrals
0 impropers

3 atom types
1 bond types
0 angle types
0 dihedral types
0 improper types
0 4.867597 xlo xhi
0 4.21546 ylo yhi
0 5.375577 zlo zhi
-2.4337985 0 0 xy xz yz

Masses

1 16.000 #Oc - Oxygen core
2 0.0125 #Os - Oxygen Shell
3 28.0855 #Si

Atoms

1 1 1 0.848 0.119085908 -0.957542974 -1.508393288
2 1 2 -2.848 0.119085908 -0.957542974 -1.508393288
3 1 1 0.848 1.986613631 -1.525829810 0.283455610
4 1 2 -2.848 1.986613631 -1.525829810 0.283455610
5 1 1 0.848 1.545000434 0.375640076 2.075304031
6 1 2 -2.848 1.545000434 0.375640076 2.075304031
7 1 1 0.848 -2.314714193 0.957542818 -0.283455372
8 1 2 -2.848 -2.314714193 0.957542818 -0.283455372
9 1 1 0.848 -0.447186351 1.525829714 -2.075304091
10 1 2 -2.848 -0.447186351 1.525829714 -2.075304091
11 1 1 0.848 -0.888799429 -0.375640112 1.508393288
12 1 2 -2.848 -0.888799429 -0.375640112 1.508393288
13 2 3 4.0 1.055879712 -2.107732612 -0.895951152
14 2 3 4.0 -2.353289843 -0.139447767 0.895951271
15 2 3 4.0 0.080510020 0.139447373 -2.68779993115

Bonds

1 1 1 2
2 1 3 4
3 1 5 6
4 1 7 8
5 1 9 10
6 1 11 12

CS-Info

1 1
2 1
3 2
4 2
5 3
6 3
7 4
8 4
9 5
10 5
11 6
12 6
13 0
14 0
15 0

Potential File

#INITIALIZATION
units metal
atom_style full

#Atom definition
fix csinfo all property/atom i_CSID
read_data data.minl.nff.SiO2 fix csinfo NULL CS-Info
group cores type 1 3
group shells type 2

#FORCE FIELDS
kspace_style ewald 1.0e-6
pair_style hybrid/overlay nb3b/harmonic buck/coul/long/cs 12.0 9.0
pair_coeff * * nb3b/harmonic SiO2.nb3b.harmonic NULL Os Si
pair_coeff * * buck/coul/long/cs 0.00 1.0 0.00
pair_coeff 2 3 buck/coul/long/cs 1283.907 0.32052 10.6620 #Os-Si
pair_coeff 2 2 buck/coul/long/cs 22764.0 0.149 27.8800 #Os-Os
bond_style harmonic
bond_coeff 1 74.92038 0.0

fix 3 all box/relax tri 1013.25 vmax 1.0e-3
thermo_style custom step temp press pe xlo xhi ylo yhi zlo zhi xy xz yz
thermo 200
minimize 1.0e-6 0.001 100000 10000

SiO2.nb3b,harmonic

ielem jelem kelem k theta_0 cutoff

Si Os Os 1.04862 109.47 1.8
Si Os Si 0.0 0.0 0.0
Si Si Os 0.0 0.0 0.0
Si Si Si 0.0 0.0 0.0
Os Os Os 0.0 0.0 0.0
Os Os Si 0.0 0.0 0.0
Os Si Os 0.0 0.0 0.0
Os Si Si 0.0 0.0 0.0

Thank You very much in advance

[1] http://rruff.info/uploads/ZK198_177.pdf

your conclusion is wrong. are you certain, that your "three-body-term"
needs to be implemented in the rather unusual way with nb3b instead of
the much more common angle potential?

to explain: since you have a system with bond and you didn't change
the special_bonds settings, pairwise interactions from 1-2, 1-3, and
1-4 neighbors are excluded from the non-bonded pairwise neighborlists,
and thus the neighborlist for nb3b will be incomplete. in fact, nb3b
was not conceived to model such a situation, but as a standalone
manybody potential.

with an inconsistent model like yours, it is no surprise, that you get
inconsistent results.

axel.

also, there seem to be far too few bonds for a crystal structure.
where are the bonds reaching across periodic boundaries?

axel.