Thank you for you prompt reply.
In both models, I didn’t specify the special_bonds command.
In fix shake model, bonds and angles are constrained by fix shake command. In the fix rigid model, bonds and angles are not specified in data file, since I think they are unnecessary for rigid bodies.
The thermo output for ebond and eanlge are all zero during the simulations. evdwl on step 0 are almost the same (9503.941 kcal/mol vs. 9503.9436 kcal/mol).
Special bonds are on and exclude neighbor interactions by
default for bonded atoms.
If you don’t turn off pairwise interactions within the rigid
bodies (see neigh modify exclude molecule) you
will compute them for rigid bodies (no bonds), but not
for fix shake (bonds with exclusions). That could
be the difference in your energies. The intra-body
forces cancel, so it shouldn’t affect dynamics.