ThermodynamicQuantities fails to remove center of mass DOF properly.

[joaander]

Description

There is a problem with this code: https://github.com/glotzerlab/hoomd-blue/blob/03e25a0a17913f55094b1e880d4ed74abc252c1c/hoomd/md/IntegratorTwoStep.cc#L159-L165

The intention is that the number of degrees of freedom of one thermo on all particles is equal to the sum of the degrees of freedom of two thermos that cover all particles with two distinct subsets. Bug 1: The code fails to perform as documented because of the group->getNumMembersGlobal() == m_pdata->getNGlobal() check. To function as documented, it should be replaced with a check against method 0's group size.

However, fixing Bug 1 as described will create another bug. When using rigid bodies, the single integration method is not applied to all particles - but only central and free particles. With ThermodynamicQuantities applied to all particles (including constituents), the ratio based on N including constituents will meet the requirement that multiple ThermodynamicQuantities instances are additive. However, the condition on when to activate the periodic degree of freedom removal needs to check only central and free particles, not constituents. At this time, an efficient solution is not clear to me.

The documentation for the degree of freedom removal that needs to be updated: https://github.com/glotzerlab/hoomd-blue/blob/0b4b6b54f2ec6acbfdc817215e89d9e9e863655c/hoomd/md/compute.py#L211-L226

While my script is built using #1424, the error is present in the most current release. Fixing the second bug is a breaking change.

Script

import hoomd

CUBE_VERTS = [(-0.5, -0.5, -0.5),
              (-0.5, -0.5, 0.5),
              (-0.5, 0.5, -0.5),
              (-0.5, 0.5, 0.5),
              (0.5, -0.5, -0.5),
              (0.5, -0.5, 0.5),
              (0.5, 0.5, -0.5),
              (0.5, 0.5, 0.5),
              ]

rigid = hoomd.md.constrain.Rigid()
rigid.body['R'] = {
    "constituent_types": ['A'] * 8,
    "positions": CUBE_VERTS,
    "orientations": [(1.0, 0.0, 0.0, 0.0)] * 8,
    }

nve = hoomd.md.methods.ConstantVolume(filter = hoomd.filter.Rigid(('center',)))
integrator = hoomd.md.Integrator(dt=0, methods=[nve], integrate_rotational_dof=True)
integrator.rigid = rigid

thermo = hoomd.md.compute.ThermodynamicQuantities(filter = hoomd.filter.Rigid(('center',)))

sim = hoomd.Simulation(device=hoomd.device.CPU())
sim.create_state_from_gsd('init.gsd')

sim.operations.integrator = integrator
sim.operations.computes.append(thermo)

sim.run(0)
print(thermo.translational_degrees_of_freedom)

Input files

init.gsd.zip

Output

1536.0

Expected output

1533.0

Platform

CPU, GPU, Linux, macOS

Installation method

Compiled from source

HOOMD-blue version

feature/modern-integrators

Python version

3.10.6

[joaander]

The fix in #1488 is not breaking.

[joaander]

Fixed in #1488.