Please check that you intended to have the C-C bond flexible instead of constrained-length and that its bond constant is correct.
Take out the line with xy yz xz if you don’t intend to deform the box. That line is forcing LAMMPS into triclinic instead of orthogonal box mode, which slows down the calculation for no benefit (and might even be the cause of the discrepancy, which would be valuable to know if true)
It would’ve been good if we could have seen all this earlier!
I don’t intend to keep the C-C bond flexible. But, in lammps using shake, I can at most keep 3 bonds rigid. That is why bond #02 is missing. Besides, I thought that bond constants values that I have chosen (600 kcal/mol) would be high enough to constrain the bond anyways. I did not know how else to use fix shake to constrain all 4 bonds.
GOMC doesn’t sample bond length changes (…why?) but your LAMMPS script has flexible C-C bonds. I imagine your GROMACS script has fixed all bond lengths (which is easier there than in LAMMPS). But given the very high bond strength I can’t imagine this being significant.
More importantly both GOMC and GROMACS use full 1-4 non-bonded interactions by default, while LAMMPS ignores them by default – that is, the default special_bonds setting is 0 0 0 in LAMMPS and (equivalent to) 0 0 1 in GOMC and GROMACS.
All the best! Since you know GROMACS reproduces the GOMC result with all bond lengths fixed, you might get the clearest comparison by running LAMMPS with 1-4 interactions turned off and GROMACS with the C-C bond length changing instead of held constant.
Also bear in mind that the model with fixed C-C bond is technically a different model from having a flexible C-C bond – but with the GROMACS runs you shouls be able to bound the difference.