Bonded Atoms flying off in simulations

GROMACS version: 2023
GROMACS modification: No

Using the PPP topology file from https://www.sciencedirect.com/science/article/abs/pii/S0892687521005094. Is there any reason why the oxygens are flying away from the bonded molecule.

Using PDB:

REMARK 1 File created by GaussView 6.0.16
HETATM 1 P1 PPP 1 -2.819 2.431 0.890 P
HETATM 2 O1 PPP 1 -2.916 4.018 0.184 O
HETATM 3 P2 PPP 1 -3.507 5.573 0.321 P
HETATM 4 O2 PPP 1 -2.848 6.250 -1.055 O
HETATM 5 P3 PPP 1 -2.703 7.670 -2.050 P
HETATM 6 O3 PPP 1 -1.292 2.200 1.159 O
HETATM 7 O4 PPP 1 -2.911 6.247 1.566 O
HETATM 8 O5 PPP 1 -3.447 8.863 -1.373 O
HETATM 9 O6 PPP 1 -3.366 1.469 -0.220 O
HETATM 10 O7 PPP 1 -3.672 2.385 2.196 O
HETATM 11 O8 PPP 1 -5.039 5.556 0.209 O
HETATM 12 O9 PPP 1 -3.335 7.269 -3.427 O
HETATM 13 O10 PPP 1 -1.158 7.913 -2.168 O
END

and

Topology:

[ atomtypes ]
; PHOSFATE
P3 15 30.97 0.0000 A 3.74177e-01 8.36800e-01
OP2 8 15.99940 0.0000 A 3.11690e-01 8.78640e-01 ;O2 p3
OP3 8 15.99940 0.0000 A 3.11690e-01 8.78640e-01 ;O3 p3
OR 8 15.99940 0.0000 A 3.15699e-01 7.11280e-01 ;OS

[ bondtypes ]
; i j func b0 kb
OP2 P3 3 0.1552 432.15350 22.110 ; GQE
OP3 P3 3 0.1552 372.15350 21.870 ; GQE
OR P3 3 0.1739 128.20180 23.119 ; GQE

[ angletypes ]
; i j k func th0 cth
P3 OR P3 1 146.639 446.488 ; GQE
OR P3 OR 1 99.813 1414.380 ; GQE
OR P3 OP2 1 109.700 1197.660 ; GQE
OR P3 OP3 1 109.700 1437.520 ; GQE
OP2 P3 OP2 1 111.746 1289.850 ; GQE
OP3 P3 OP3 1 116.822 1333.210 ; GQE

[ dihedraltypes ]
;i j k l func phase kd pn
OP3 P3 OR P3 9 180.0 4.11810 1 ; GQE
OP3 P3 OR P3 9 180.0 2.22870 3 ; GQE
OP2 P3 OR P3 9 245.0 47.38680 1 ; GQE
OP2 P3 OR P3 9 13.0 1.93320 3 ; GQE
OR P3 OR P3 9 0.0 45.77210 1 ; GQE
OR P3 OR P3 9 0.0 3.09440 2 ; GQE
OR P3 OR P3 9 0.0 1.91100 3 ; GQE

[ moleculetype ]
; Name nrexcl
PPP 3
[ atoms ]
; nr type resnr resid atom cgnr charge mass
1 P3 1 PPP P 1 +1.269450 30.9700
2 OP3 1 PPP O 1 -0.995975 15.9994
3 OP3 1 PPP O 1 -0.995975 15.9994
4 OP3 1 PPP O 1 -0.965300 15.9994
5 OR 1 PPP O 1 -0.560050 15.9994
6 P3 1 PPP P 1 +1.333700 30.9700
7 OP2 1 PPP O 1 -0.919000 15.9994
8 OP2 1 PPP O 1 -0.919000 15.9994
9 OR 1 PPP O 1 -0.560050 15.9994
10 P3 1 PPP P 1 +1.269450 30.9700
11 OP3 1 PPP O 1 -0.995975 15.9994
12 OP3 1 PPP O 1 -0.995975 15.9994
13 OP3 1 PPP O 1 -0.965300 15.9994
;total charge -5.000
[ bonds ]
; ai aj funct c0 c1
1 2 3 ; P3 OP3
1 3 3 ; P3 OP3
1 4 3 ; P3 OP3
1 5 3 ; P3 OR
5 6 3 ; OR P3
6 7 3 ; P3 OP2
6 8 3 ; P3 OP2
6 9 3 ; P3 OR
9 10 3 ; OR P3
10 11 3 ; P3 OP3
10 12 3 ; P3 OP3
10 13 3 ; P3 OP3
[ pairs ]
; ai aj funct
1 7 1
1 8 1
1 9 1
2 6 1
3 6 1
4 6 1
5 10 1
6 11 1
6 12 1
6 13 1
7 10 1
8 10 1
[ angles ]
; ai aj ak funct angle fc
2 1 3 1
2 1 4 1
2 1 5 1
3 1 4 1
3 1 5 1
4 1 5 1
1 5 6 1
5 6 7 1
5 6 8 1
5 6 9 1
7 6 8 1
7 6 9 1
8 6 9 1
6 9 10 1
9 10 11 1
9 10 12 1
9 10 13 1
11 10 12 1
11 10 13 1
12 10 13 1
[ dihedrals ]
; ai aj ak al funct ph0 cp mult
1 5 6 7 9 ; P3 OR P3 OP2
1 5 6 8 9 ; P3 OR P3 OP2
1 5 6 9 9 ; P3 OR P3 OR
2 1 5 6 9 ; OP3 P3 OR P3
3 1 5 6 9 ; OP3 P3 OR P3
4 1 5 6 9 ; OP3 P3 OR P3
5 6 9 10 9 ; OR P3 OR P3
6 9 10 11 9 ; P3 OR P3 OP3
6 9 10 12 9 ; P3 OR P3 OP3
6 9 10 13 9 ; P3 OR P3 OP3
7 6 9 10 9 ; OP2 P3 OR P3
8 6 9 10 9 ; OP2 P3 OR P3
[ dihedrals ]
; ai aj ak al funct angle fc
#ifdef POSRES1
[ position_restraints ]
; atom type fx fy fz
1 1 10000 10000 10000
2 1 10000 10000 10000
3 1 10000 10000 10000
4 1 10000 10000 10000
5 1 10000 10000 10000
6 1 10000 10000 10000
7 1 10000 10000 10000
8 1 10000 10000 10000
9 1 10000 10000 10000
10 1 10000 10000 10000
11 1 10000 10000 10000
12 1 10000 10000 10000
13 1 10000 10000 10000
#endif

Are you sure the atoms are actually flying off and that it is not only a visualization artifact?

You should also make sure that the atoms in the coordinate input file (pdb file) are in the same order as in the topology. I guess you would have gotten a warning about that from gmx grompp.