Tabulated pair potentials not working properly in single precision mode

[ara137]

Description

The tabulated pair potential does not work properly when compiling HOOMD-blue v.4.7.0 in single precision mode. We tested this issue by creating a tabulated version of the Lennard-Jones potential, placing two particles in a simulation box, and iteratively computing their potential energy as a function of their distance. This procedure reproduces the correct result when using double precision, but fails when using single precision.

Script

import numpy as np
import gsd
import hoomd
import sys

# Read self-made potentials
data_array = np.loadtxt('potential.csv', delimiter=' ')
potential = data_array[:,1]
force = data_array[:,2]

# Initialize device and read starting configuration from file
dev = hoomd.device.GPU() #CPU
simulation = hoomd.Simulation(device=dev, seed=42)
simulation.timestep = 0
simulation.create_state_from_gsd(f'start.gsd', frame=-1)
box = simulation.state.box

#Create cell list
cell = hoomd.md.nlist.Cell(buffer=0.4)

#Non-bonded potential
self_made = hoomd.md.pair.Table(default_r_cut=5.0, nlist=cell);
self_made.params[('A', 'A')]=dict(r_min=0.0, U=potential ,F=force);

#Group
stationary_particles = hoomd.filter.Type(['A'])
mobile_particles = hoomd.filter.SetDifference(hoomd.filter.All(),  stationary_particles)
all_particles = hoomd.filter.All()

# Create intergator & add interactions to the integrator
langevin = hoomd.md.methods.Langevin(filter=mobile_particles, kT=1)
integrator = hoomd.md.Integrator(dt=0.005)
integrator.forces.append(self_made)
integrator.methods.append(langevin)
simulation.operations.integrator = integrator

# Prepare trajectory output
gsd_writer = hoomd.write.GSD(filename=f'traj_tabular.gsd', trigger=hoomd.trigger.Periodic(1), mode='wb')
simulation.operations.writers.append(gsd_writer)

# Prepare logger output
thermodynamic_properties = hoomd.md.compute.ThermodynamicQuantities(filter=all_particles)
simulation.operations.computes.append(thermodynamic_properties)

# Create logger of thermodynamic quantities to file
logger = hoomd.logging.Logger(categories=['scalar'])
logger.add(simulation, quantities=['timestep'])
logger.add(thermodynamic_properties, quantities=['kinetic_temperature', 'potential_energy'])
logger.add(obj=self_made, quantities=['energy'])
file = open(f'log_tabular.dat', mode='w', newline='\n')
table_file = hoomd.write.Table(output=file, trigger=hoomd.trigger.Periodic(period=1), logger=logger)
simulation.operations.writers.append(table_file)

# Create logger of performance to screen
logger_perf = hoomd.logging.Logger(categories=['scalar'])
logger_perf.add(simulation, quantities=['timestep', 'tps', 'walltime'])
table_perf = hoomd.write.Table(trigger=hoomd.trigger.Periodic(period=1), logger=logger_perf)
simulation.operations.writers.append(table_perf)

#Move Particle
class MoveParticle(hoomd.custom.Action):
    def act(self, timestep):
        snapshot = simulation.state.get_snapshot()
        snapshot.particles.position[1,0] += 0.01
        simulation.state.set_snapshot(snapshot)

action0 = MoveParticle()
moving_particle = hoomd.update.CustomUpdater(action=action0, trigger=hoomd.trigger.Periodic(1))
simulation.operations.updaters.append(moving_particle)

# Integrate the equations of motion of the mobile particles.
try:
    simulation.run(500)
finally:
    gsd_writer.flush()

Input files

potential.csv start.zip

Output

Erroneous potential energy, which does not resemble the tabulated pair potential at all.

Expected output

We expected to get the input pair potential (cf. the curves for "LJ pot single precision", "LJ pot double precision", and "Tabular LJ pot double precision" in attached graph). bugCheck.pdf

Platform

GPU, Linux

Installation method

Compiled from source

HOOMD-blue version

4.7.0

Python version

3.11.8

[joaander]

Many operations work incorrectly when you build HOOMD-blue in single precision, so I do not offer any support for this configuration.

You are welcome to submit a pull request that fixes the table potential bug. I recommend that first check that the table values are translated correctly from C++ to Python - the bug may be in that conversion.

[ara137]

That's understandable. We'll try to have a more detailed look into it and work on a solution.