I was trying to simulate fully rigid proteins using HOOMD. For that, I have tried to test a simpler system first. It is a system of three triangular molecules sitting on each other (see picture). My plan was to use LJ potential so I have calculated an analytical solution for this system. So, it was equal to -0.6341327 kcal/mol after excluding intramolecular interactions.
Then, I have created this system with 6 particles in HOOMD. In HOOMD, I have also created 2 central particles with required properties (tag, body, moment of inertia, orientation, mass etc.) and defined rigid chains as below:
rigid = hoomd.md.constrain.rigid();
When I run this fully rigid case, I am getting some discrepancy in potentials between HOOMD and my analytical solution. The pair_lj_energy from HOOMD was -0.258301451 kcal/mole.
Could you please help me with this issue? I have attached my HOOMD simulation file (run_rigid.py). Thanks in advance.