Precision loss when writing GSD files

[CalCraven]

Description

When using charge normalization for elementary_charges, floating point errors result in total charges that don't quite cancel out, something close to 1e-17. I think the primary cause for this is the charge factor, 0.0848385920 e, which is pretty much unavoidable.

This small of a rounding issue doesn't seem to cause a problem, unless the snapshot is written to a gsd file, and then loaded back into a simulation using sim.create_state_from_gsd. The total charge jumps from 1e-16 to 1e-7, and this can cause charge issues in the simulation. I especially see this when I run a system with >5000 molecules, which should be charge neutral, but end up with a net charge.

I've tried using the Decimal to improve the float precision to set to 64 bit precision during the initial conversion from the elementary_charge units to the dimensionless units, but saw no change in behavior.

from decimal import Decimal, getcontext

getcontext().prec = 64  # Set desired precision

Maybe there should be a way to charge balance the hoomd state once something is loaded in? Even just shifting every charge by the factor of total_charge/sim.state.N_particles might counteract issues like this?

Script

import hoomd 
import gsd.hoomd
import numpy as np

print(f"Hoomd version: {hoomd.version.version}")
print(f"GSD version: {gsd.version.version}")

# NOTE: first charges fails as well, but has more particles so simplified to just 4 particles.
# charges = [-0.18,  -0.06,  -0.18,  0.51,  -0.43,  -0.33,  0.16,  0.06,  0.06,  0.06,  0.06,  0.06,  0.06,  0.06,  0.03,  0.03,  0.03]
charges = [-0.1, 0.1/3, 0.1/3, 0.1/3]
print(f"Initial float sum for charges with single molecule: {sum(charges)}")
charge_factor = 0.0848385920 # from https://hoomd-blue.readthedocs.io/en/latest/units.html
print(f"Initial normalized sum for charges with single molecule: {sum(charges)/charge_factor}")

snapshot = gsd.hoomd.Frame()
box = [10, 10, 10, 0, 0, 0] # Lx, Ly, Lz, xy, xz, yz
n_mols = 100
snapshot.particles.N = len(charges)*n_mols # replicates highlight charge issue
positions = (np.random.rand(snapshot.particles.N , 3) - 0.5) * 10
snapshot.particles.position = positions
all_charges = np.array(charges*n_mols) /charge_factor
print(f"Final charges plugged in: {sum(all_charges)}")
snapshot.particles.charge = all_charges
print(f"Ending total charge {sum(snapshot.particles.charge)}")

# NOTE: fill in meaningless data
# Define particle types
snapshot.particles.types = ['A']
# Assign type 'A' (typeid 0) to all particles
snapshot.particles.typeid = [0] * snapshot.particles.N 

# save to gsd
device = hoomd.device.CPU()
sim = hoomd.Simulation(device=device, seed=1)
sim.create_state_from_snapshot(snapshot)
hoomd.write.GSD.write(
    state=sim.state, filename="checkcharge.gsd", mode="wb"
)
del sim

# reload gsd
device = hoomd.device.CPU()
sim = hoomd.Simulation(device=device, seed=1)
sim.create_state_from_gsd("checkcharge.gsd")
snap_loaded = sim.state.get_snapshot()
print(f"Reloaded charge: {sum(snap_loaded.particles.charge)}")

Input files

No response

Output

Hoomd version: 6.0.0
GSD version: 4.2.0
Initial float sum for charges with single molecule: -6.938893903907228e-18
Initial normalized sum for charges with single molecule: -8.178935718201485e-17
Final charges plugged in: -1.1102230246251565e-16
Ending total charge -1.1102230246251565e-16
Reloaded charge: -8.940696716308594e-06

Expected output

Hoomd version: 6.0.0
GSD version: 4.2.0
Initial float sum for charges with single molecule: -6.938893903907228e-18
Initial normalized sum for charges with single molecule: -8.178935718201485e-17
Final charges plugged in: -1.1102230246251565e-16
Ending total charge -1.1102230246251565e-16
Reloaded charge: -1.1102230246251565e-16

Platform

macOS, CPU

Installation method

Conda-forge package

HOOMD-blue version

6.0.0

Python version

3.12.8

[joaander]

GSD stores charge (and all other floating point quantities) in single precision: https://gsd.readthedocs.io/en/v4.2.0/schema-hoomd.html#chunk-particles-charge.

This simplified script reproduces the same behavior:

import numpy as np

charges = [-0.1, 0.1/3, 0.1/3, 0.1/3]
charge_factor = 0.0848385920 # from https://hoomd-blue.readthedocs.io/en/latest/units.html

n_mols = 100
a = np.array(charges*n_mols) /charge_factor
b = np.array(a, dtype=np.float32)
c = np.array(b, dtype=np.float64)

print(f"sum(a): {sum(a)}")
print(f"sum(b): {sum(b)}")
print(f"sum(c): {sum(c)}")

sum(a): -1.1102230246251565e-16
sum(b): -1.7881393432617188e-07
sum(c): -8.940696716308594e-06

To achieve the behavior you intend, you could initialize directly using create_state_from_snapshot, or you could store your chargers in a .npy or other format that preserves float64 precision.

Changing the GSD specification to allow either float32 or float64 quantities is another possibility. It would require a breaking change to the spec and coordination with all GSD readers/writers (OVITO, gsd-vmd, HOOMD-blue, and hoomd-rs). I do not have the bandwidth to implement such a feature myself, but would accept contributions from the community. The float64 precision write mode should be optional as it would approximately double the file size.

[CalCraven]

Thanks @joaander, .npy should handle my use case for now, since our simulations are broken up into individual simulation a few runs, so rerunning the snapshot generation using MoSDeF tools would be time intensive for each simulation step.

Changing the GSD specification to allow either float32 or float64 quantities is another possibility.

I'm thinking this would be the most useful overall change. I'll look into making the contribution myself, as it would be a handy option for a lot of the charge conversion in the MoSDeF space. In this case, it's a single Frame, and so doubling the file size wouldn't be a huge deal, my current size is 2.4MB for a large polymer system. I think the question I would need to think about would be if there is reason to initialize a system using float64, but then store the main trajectory gsd as float32 (since this is the main data usage), and then also save the final frame of the system in float64 for future continued simulation use.

[mphoward]

if there is reason to initialize a system using float64, but then store the main trajectory gsd as float32 (since this is the main data usage), and then also save the final frame of the system in float64 for future continued simulation use.

Yes, I think there definitely is for simulation restarts! HOOMD doesn't guarantee binary reproducibility of runs, but the loss of precision from writing and then reading a restart file can still cause some pretty big issues in some cases (e.g., loss of momentum conservation). We work around the momentum issue when it matters by manually re-zeroing via the snapshot (like you could do for charge), but I would be thrilled to write fp64 restart files with fp32 trajectories for analysis.

[mphoward]

(Btw, I am also the developer of gsd-vmd, and I would be happy to support configurable precision in the GSD specification if it were made available.)

[joaander]

The steps to achieve this would be:

1. In https://github.com/glotzerlab/gsd:

• Update the HOOMD schema specification. The simplest approach would be to allow either float32 or float64 data anywhere float32 is used now. Different frames could be stored at different precision, and so could different data chunks in a single frame.
• Support this specification in the reader and writer.
• Unit test the changes.

2. In https://github.com/glotzerlab/hoomd-blue:

• Update GSDReader to support the new specification.
• Add a precision attribute to GSDWriter that controls the precision written to the file.

3. Produce an example file that uses a mix of float32/float64 types.

• Open an issue and request that OVITO add support for the new specification. Provide the example file.
• @mphoward can update gsd-vmd.
• I will implement the new spec in hoomd-rs.

[mphoward]

The simplest approach would be to allow either float32 or float64 data anywhere float32 is used now. Different frames could be stored at different precision, and so could different data chunks in a single frame.

If possible, it would be great to have the precision default to float32 in the file for backward compatibility, then have an option to set the default precision of a frame, then last an option to set the precision of a chunk in a frame.