MPI Differences with Rotational Degrees of Freedom in ALJ

[Haoyuan-Shi]

Hello,

I am working on a cubes system using ALJ without thermostat in HOOMD 5.2.0. The codes in md2.py are:

import os
import numpy as np
import hoomd

cpu = hoomd.device.CPU()
simulation = hoomd.Simulation(device=cpu, seed=0)
simulation.create_state_from_gsd(filename='lattice.gsd')
simulation.state.replicate(2, 2, 2)
integrator = hoomd.md.Integrator(dt=1.0E-4, integrate_rotational_dof=True)

cell = hoomd.md.nlist.Cell(buffer=0.4)
alj = hoomd.md.pair.aniso.ALJ(nlist=cell, default_r_cut=2.5)
alj.params[('cube', 'cube')] = dict(
    epsilon=1, 
    sigma_i=2, 
    sigma_j=2, 
    alpha=3
)

alj.r_cut[('cube', 'cube')] = 2.5

vertices = [
    (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),
]
faces = [
    (2, 1, 0),
    (0, 3, 2),
    (5, 6, 4),
    (4, 6, 7),
    (1, 2, 5),
    (5, 2, 6),
    (4, 7, 3),
    (3, 0, 4),
    (0, 1, 5),
    (5, 4, 0),
    (3, 7, 6),
    (6, 2, 3)
]

alj.shape['cube'] = dict(
    vertices=vertices,
    rounding_radii=(0.15, 0.15, 0.15),
    faces=faces
)

integrator.forces.append(alj)

nvt = hoomd.md.methods.ConstantVolume(
    filter=hoomd.filter.All(),
    thermostat=None
)

integrator.methods.append(nvt)

simulation.operations.integrator = integrator

thermodynamic_properties = hoomd.md.compute.ThermodynamicQuantities(
    filter=hoomd.filter.All()
)
simulation.operations.computes.append(thermodynamic_properties)

logger = hoomd.logging.Logger()
logger.add(
    thermodynamic_properties,
    quantities=[
        'kinetic_temperature',
        'kinetic_energy',
        'translational_kinetic_energy',
        'rotational_kinetic_energy',
        'potential_energy'
    ]
)
logger.add(alj, quantities=['type_shapes', 'energies', 'forces', 'torques'])

gsd_writer = hoomd.write.GSD(
    filename='trajectoryCubes.gsd', 
    trigger=hoomd.trigger.Periodic(100), 
    mode='xb', 
    filter=hoomd.filter.All(), 
    dynamic=['property', 'momentum', 'attribute'],
    logger=logger
)

simulation.operations.writers.append(gsd_writer)
simulation.run(2.0E4)
gsd_writer.flush()
print(thermodynamic_properties.kinetic_energy)

When I ran the simulation with different MPI process counts:

mpirun -n 1 python md2.py & mpirun -n 2 python md2.py

mpirun -n 4 python md2.py & mpirun -n 8 python md2.py

The results differ depending on the number of MPI processes. However, when I set integrate_rotational_dof=False, all four MPI runs produced exactly the same energy results:

This leads me to suspect that MPI might not properly handle the rotational degrees of freedom in ALJ, but it could also be due to incorrect settings on my part. Any comments or suggestions would be greatly appreciated.

Attached are md2.py, lattice.gsd, and plot.py.
ALJcubes.zip

Thank you,
Haoyuan Shi

[joaander]

Your r_cut is too small, leading to the evaluation of discontinuous forces that destabilize the integrator. See the "Choosing r_cut" section in the documentation: https://hoomd-blue.readthedocs.io/en/v5.2.0/hoomd/md/pair/aniso/alj.html .

Also, you should not use rounding_radii in conjunction with a polyhedron. The sigma values already round the polyhedron. rounding_radii is intended for the creation of ellipsoids via ALJ.

[Haoyuan-Shi]

Thank you, Joshua! I now understand the issue with discontinuous forces, and I’ve increased r_cut and removed rounding_radii as suggested.

However, my main concern is that the total energy (translational + rotational + potential energy) results vary depending on the number of MPI ranks, as shown below:

When I set integrate_rotational_dof=False, the total energy (translational + potential energy) remains exactly the same across different MPI configurations:

Thank you again for your help!
Haoyuan

[joaander]

If you want the highest accuracy for simulations of shapes, use HPMC: https://hoomd-blue.readthedocs.io/en/v5.2.0/hoomd/hpmc/integrate/convexspheropolyhedron.html . HPMC will also perform much faster than ALJ with very small dt values. Read on to learn about some of the difficulties with ALJ.

You should not expect identical trajectories. MD simulations are chaotic. Even a single bit changed in one of the force calculations will cause an exponential divergence in the following trajectory. When you change the number of ranks, you change the order in which forces are summed which causes the net force to be rounded differently. If using the GPU, simulations even on a single GPU will not be identical from one execution to the next. What you can expect is that simulations with thermostats produce the same average thermodynamic quantities, once they are equilibrated.

The fact that the total energy is not conserved in ALJ is concerning, however I see similar behavior in my validation tests. The ALJ developer did everything possible to make the potential conservative - however there are fundamental issues with discontinuities in the torque. In practice, these are relatively small and can be absorbed by a thermostat.

You might also try making dt smaller. Most users of ALJ find that they need a surprisingly small value of dt to stabilize their simulations. The tips of a rotating body have a faster linear velocity than an equivalent point particle with the same potential and thus a smaller dt is needed to accurately handle any collisions.

Here are the energy conservation tests from the hoomd-validation tests. These are a little more stable than yours, but still show several discontinuous jumps in the energy:

The momentum jumps are also concerning. One user has an ALJ simulation where the whole system picks up an appreciable velocity in one direction, but we haven't been able to isolate the cause of the bug.

[Haoyuan-Shi]

Thanks for your reply, Joshua! We'll continue testing more cases and studying the code to deepen our understanding.