2D Gay-Berne work on CPU but not on GPU

[Syrocco]

I just started to use HOOMD this morning, so far, it's really intuitive and easy to work with.
I assembled a bunch of example code to make a simulation of a Gay-Berne fluid in 2D:

import gsd.hoomd
import hoomd
import numpy as np
import itertools
import datetime

    
class StatusF:
    def __init__(self, potential):
        self.potential = potential
        
    @property
    def zforce(self):
        return np.max(np.abs(self.potential.forces[:, 2]))
    @property
    def xytorque(self):
        return np.max(np.abs(self.potential.torques[:, :-1]))
    
L = 20
spacing = 1
x = np.arange(-L/2, L/2 - spacing, spacing)
pos = list(itertools.product(x, repeat=2))
position = []
for x, y in pos:
    position.append((x, y, 0))
N_particles = len(position)
frame = gsd.hoomd.Frame()
frame.particles.N = N_particles
frame.particles.position = position[0:N_particles]
frame.particles.typeid = [0] * N_particles
frame.configuration.box = [L, L, 0, 0, 0, 0]
frame.particles.types = ['A']
orientation = []
for i in range(N_particles):
    angle = np.random.uniform(0, 2 * np.pi)
    orientation.append([np.sqrt(2)/2,  np.sqrt(2)/2,  0,  0]) 

frame.particles.orientation = orientation
frame.particles.moment_inertia = [(0, 1, 0) for _ in range(N_particles)]  


with gsd.hoomd.open(name='lattice.gsd', mode='w') as f:
    f.append(frame)

device = hoomd.device.GPU()
simulation = hoomd.Simulation(device=device, seed=1)
simulation.create_state_from_gsd(filename='lattice.gsd')



integrator = hoomd.md.Integrator(dt=0.005, integrate_rotational_dof=True)
cell = hoomd.md.nlist.Cell(buffer=0.4)
gay_berne = hoomd.md.pair.aniso.GayBerne(nlist=cell, default_r_cut=2.5)
gay_berne.params[('A', 'A')] = dict(epsilon=1, lperp=0.2, lpar=0.5)
gay_berne.r_cut[('A', 'A')] = 2 ** (1.0 / 6.0)
integrator.forces.append(gay_berne)

nvt = hoomd.md.methods.Langevin(hoomd.filter.All(), kT = 1)
integrator.methods.append(nvt)
simulation.operations.integrator = integrator
simulation.run(0)

logger = hoomd.logging.Logger(categories=['scalar'])
status = StatusF(gay_berne)
logger[('StatusF', 'zforce')] = (status, 'zforce', 'scalar')
logger[('StatusF', 'xytorque')] = (status, 'xytorque', 'scalar')

table = hoomd.write.Table(trigger=hoomd.trigger.Periodic(period = 1), logger=logger)

simulation.operations.writers.append(table)
logger2 = hoomd.logging.Logger()
logger2.add(gay_berne, quantities=['type_shapes'])


gsd = hoomd.write.GSD(trigger = 20, filename='random.gsd', mode = 'wb', dynamic =  ["particles/orientation","attribute", "property", "momentum"], logger = logger2)
simulation.operations.writers.append(gsd)


simulation.run(1000)
gsd.flush()

The simulation works well when I'm using the CPU (device = hoomd.device.CPU()). But with the GPU (device = hoomd.device.GPU()), the particles after a little bit of time start to go out of the box (in the z direction):

GPU:

CPU:

This does not seem to depend on the method of integration.

---

I struggled a bit with the orientations and the moment of inertia to have all my particles lying on the xy plane and correctly rotate along the z axis. Perhaps, the issue is here?

[joaander]

It would appear that the original author of the Gay-Berne code did not consider 2D simulations when implementing it - the parallel axis is pointing in the z direction (in particle local space). This is why you needed a rotation about x to place the long axis in the plane and set $I_{yy} \ne 0$. It would have made more sense if HOOMD-blue used x or y for the parallel axis. Then your quaternions will all be rotations about z and $I_{zz} \ne 0$.

HOOMD-blue computes all simulations in a 3D space. 2D simulations assume that $F_z = 0$ and that $\tau_x = \tau_y = 0$. This works generally well, but not in your case. It appears that some $F_z \ne 0$ has been added to the particles which Langevin faithfully integrated. The compiled GPU code may be rounding results differently than on the CPU (it is the same source code for GB on both). Because of the rotation about x, the particles need to be rotated into some local reference frame where the force is computed and then the force must be rotated back. Though, it is possible that the bug does not result from the GB code, but somewhere else.

Thanks for the complete script. I'll take a closer look in a few days. If you have time to try more in the meantime, here are some things I would check:

• Starting in an initial condition where all particles are in the plane (e.g. the last frame of the CPU simulation) - check the output gay_berne.forces and see if F_z is non-zero for any particles. Also check gay_berne.torque
• Double check that a LJ pair force (in place of the GB) stays in the plane.
• Try GB with I = [0, 0, 1] and orientation = [1, 0, 0, 0].

A non-zero $F_z$ from GB in any of these cases would confirm my suspicion. Assuming that the value is very small, it may not be possible to fix the numerical problem in GB. We may need to modify all MD integrators to effectively set $F_z = 0$ in 2D simulations. They already do this for torques based on the value of $I$ (this is why your particles don't rotate out of the frame), so expending that to $F_z$ would be consistent:

hoomd-blue/hoomd/md/TwoStepLangevin.cc

Lines 107 to 122 in 2d9d4df

	// rotate torque into principal frame
	t = rotate(conj(q), t);
	// check for zero moment of inertia
	bool x_zero, y_zero, z_zero;
	x_zero = (I.x == 0);
	y_zero = (I.y == 0);
	z_zero = (I.z == 0);
	// ignore torque component along an axis for which the moment of inertia zero
	if (x_zero)
	t.x = 0;
	if (y_zero)
	t.y = 0;
	if (z_zero)
	t.z = 0;

[joaander]

Also, be aware of this reported issue with GB: #1598. I welcome pull requests to fix the problem, preferably from a user familiar with knowledge of the correct procedure for the cutoff application.

[Syrocco]

Thanks for the incredibly fast answer.
I edited the main post to include the logger in the program. Indeed, with CPU, max($F_z$) and max($\tau_{xy}$) are 0 with time. On the GPU, it increases with time, starting from a small but non-zero value:
StatusF.zforce StatusF.xytorque
1.29579e-14 1.34392e-14
8.47241e-15 7.55947e-15
8.15968e-15 8.35148e-15
9.01403e-15 8.10719e-15
............... ...............
............... ...............
............... ...............
23.95944 6.36492
26.41064 6.53523
26.01478 6.15156
21.88178 10.55814
30.64415 15.52885

I compiled HOOMD myself, maybe I messed up some configurations on the floating point precision on the GPU? Indeed, with I = [0, 0, 1] and orientation = [1, 0, 0, 0], my particles stay on the plane; fz and torquexy start and stay at 0 on GPU and CPU. Is there a way to create a 2D box along x and z or y and z?

---

LJ stays in plane on CPU and GPU with max fz and torquexy = 0 at all time.

[Syrocco]

On an unrelated note (sorry I maybe should create an other discussion?). I'm just trying to simulate non-attracting (WCA like) ellipsoidal particles. Maybe the pair potential ALJ is more suited for this kind of thing? I did not manage to understand how to create ellipses with this pair potential though, so I went with Gay-Berne.

[joaander]

A ~1e-14 force in the 0th frame is enough to start pushing the particles out of the plane. I think the best solution at this point would be to modify all the integrators to force F_z = 0 when in 2D simulations. With that change, your forces should stay bounded and very small. Would you be interested in contributing this change? If not, I can get to it at some point in the next few weeks.

To make an ellipsoid in ALJ, use a single vertex shape with rounding_radii = [a, b, c]. Note that ALJ computes a single interaction at the closest point of contact. I'm not familiar enough with GB to know whether it can implement purely repulsive interactions.

HPMC's ellipsoid shape is by far the easiest way to model volume excluding ellipsoid shapes: https://hoomd-blue.readthedocs.io/en/v4.8.2/module-hpmc-integrate.html#hoomd.hpmc.integrate.Ellipsoid

[Syrocco]

Ok thank you very much! I think GB works fine for purely repulsive with the cutoff à la WCA provided by the pair interaction.
Unfortunately, I'm not comfortable enough (at all) with C++ and the codebase to contribute myself to this change :/.

But I think I will keep learning about HOOMD (it's really a breeze of fresh air after struggling with LAMMPS for 2 years..) and hopefully, someday be able to contribute myself to the project.

Thanks again, for your answer and for this very nice tool!

[joaander]

I opened #1907 to track the implementation of this fix.