Memory Leaks during Simulations with increasing timesteps

[Soulchemist97]

Hello everyone,

I'm currently running HOOMD-blue 5.1.1 simulations using Python 3.12 (in Windows 11, WSL2, Linux backend, Ubuntu) and
I have a Problem with my simulation regarding scaling in size and particles due to memory leaks.
This code works perfectly fine for smaller Simulations , so my System and hardware limitations are not the problem.
But my Memory usage increases with each timestep monitored with print_memory_usage().
After a couple of steps my simulation hits multiple GB of RAM and ends with the word: Killed.
I also do not think it is regarding the additional Packages (Loguru for external logging and tqdm for progress bar).

What I’ve already tried:

• Reducing the chunk size (i.e. how many steps I run at once) → no improvement
• Increasing the frequency of gsd_writer.flush() → no effect
• Disabling logging (logger=None) and minimizing output
• Adding explicit gc.collect() calls in the loop → no reduction in memory usage

What could I do to free my memory during the simulation to scale it for bigger and longer simulations ?

Thanks in advance and best regards

# ... simulation setup, integrator and so on

def print_memory_usage():
    process = psutil.Process(os.getpid())
    mem = process.memory_info().rss / 1024 / 1024  # in MB
    print(f"[MEM] RAM usage: {mem:.2f} MB")

gsd_trajectory_writer = hoomd.write.GSD(
    filename="Trajectory.gsd",
    trigger=hoomd.trigger.Periodic(write_step),   
    mode="wb",                             
    filter=hoomd.filter.All(),            
    logger=logger
)
gsd_trajectory_writer.write_diameter = True
simulation.operations.writers.append(gsd_trajectory_writer)

if __name__ == "__main__":
    chunk: int = int(1e3)
    n_steps: int = int(1e5)

    flush_interval = 10
    for Sim_step in tqdm(range(n_steps // chunk),desc="Simulation"):
        try:
            simulation.run(chunk) 

            # Save every N chunks  
            if Sim_step % flush_interval == 0:
                gsd_trajectory_writer.flush()
                logging.info(f"Saved Trajectory at step:{simulation.timestep:_}/{n_steps:_}")

        except (KeyboardInterrupt,Exception) as e:
            logging.info(f"Simulation error: {e}")
            gsd_trajectory_writer.flush()
            logging.info(f"Stopped at Simulation step: {simulation.timestep}")
            break
        finally: 
            print_memory_usage()

[joaander]

HOOMD-blue regularly sees production use on billion time step simulations. There is nothing in the package that inherently increases memory as a function of step. There are, however, ways you might use HOOMD-blue outside normal parameters that could cause increasing memory usage. You don't give enough information to even guess if that is the case, or if it is something else in your script.

I suggest that you use memray (https://github.com/bloomberg/memray) to what library/class/function is consuming so much memory. If that consumer is in HOOMD-blue, then I can help you identify ways to mitigate the usage.

[Soulchemist97]

The problem seems to be my custom potential.
Without the line integrator.forces.append(Solvophobicity_potential)
everything runs fine.
Here is my code regarding my custom potential:

def V_alpha(r,alpha,epsilon=epsilon,sigma=sigma,R0=1.5):
    gamma = pi/(2.25-2**(1/3))
    beta = 2*pi-2.25*gamma
    V_Alpha = 0.5*alpha*epsilon*(cos(gamma*(r/sigma)**2+beta)-1)

    V_Alpha = np.where(r <= 2**(1/6)*sigma, V_Alpha, -epsilon * alpha)
    V_Alpha = np.where(R0*sigma <= r, V_Alpha, 0)
    return V_Alpha

rmin_alpha=0.1
rcrit_alpha = 2

alpha=0

r_alpha = np.linspace(start=rmin_alpha,stop=rcrit_alpha,num=500) #np.arange(rmin,rcrit,0.005)
U_alpha = V_alpha(r_alpha,alpha)
F_alpha = -np.gradient(U_alpha,r_alpha) #Force as derivative

Solvophobicity_potential = hoomd.md.pair.Table(nlist=Neighbor_list,default_r_cut=rcrit_alpha)
Solvophobicity_potential.params[["A","B"],["B","B"],] = dict(r_min=rmin_alpha,U=U_alpha,F=F_alpha)
Solvophobicity_potential.params["A","A"] = dict(r_min=rmin_alpha,U=U_alpha,F=F_alpha)

stepsize= 5e-4
integrator = hoomd.md.Integrator(dt=stepsize,integrate_rotational_dof=True)
simulation.operations.integrator = integrator 

# --- Forces --- #
integrator.forces.append(fenewca)
integrator.forces.append(Solvophobicity_potential)

According to memray, the problem are in

• get_available_devices from site-packages.hoomd.device (Own Bytes 40 GB)
• run from site-packages.hoomd.simulation (Own Bytes 7.9 GB)

If the reason could be found in my device logic I used the following code to distinguish between WSL and native Linux. :
If I understand it correctly, HOOMD GPU Support is not available under Windows because of something named united Memory right ?

# ---- build simulation ----- #
GPU_List = hoomd.device.GPU.get_available_devices()
print("GPU List:",GPU_List)
cpu = hoomd.device.CPU()

def is_wsl()->bool:
    """
    Checks if script is running in WSL (Windows Subsystem for Linux).
    Returns:
        bool: True if running in WSL, False otherwise.
    """

    if "microsoft" in platform.uname().release.lower():
        return True
    try:
        with open("/proc/version", "r") as f:
            version = f.read().lower()
            if "microsoft" in version:
                return True
    except FileNotFoundError:
        pass
    return False

Operating_System:str = platform.system()
print("OS:",Operating_System)
isWSL:bool = is_wsl()

print("Setting up HOOMD-blue simulation on CPU or GPU")
# Check if GPU is available
# GPU Support only on native LINUX
if len(GPU_List) > 0 and not isWSL:
    print("Running on GPU")
    gpu = hoomd.device.GPU()
    gpu.gpu_error_checking = True
    simulation = hoomd.Simulation(device=gpu, seed=0)
elif isWSL: 
    print("Running in WSL, using CPU")
    simulation = hoomd.Simulation(device=cpu, seed=0)
else:
    simulation = hoomd.Simulation(device=cpu, seed=0)

simulation.create_state_from_gsd(filename="Gel_Particle.gsd")

I hope this could explain my problem better and you are able to help me now.

Thanks for your help and your time.

[joaander]

I would ignore any memory usage reported by get_available_devices. The CUDA driver maps a lot of memory, but doesn't allocate it or consume it. You are correct that HOOMD-blue does not support GPUs under WSL due to limitations in the driver. You can install a CPU-only build of HOOMD-blue to avoid seeing this in the memray profile.

I'm not aware of any memory leaks in the Table potential. It is possible that the memory utilization comes from a data structure that Table uses. Do you get the same memory utilization if you replace Table with a LJ potential that has the same r_cut?

Please also run memray with the native mode enabled and share the full flamegraph html file so that I can see exactly which C++ objects are allocating the memory.

memray run --native script.py options

memray flamegraph memray-*.bin

[Soulchemist97]

I tried to switch my custom potential with a Lennard-Jones Potential with the same r_cut and get the same problem with increased memory usage.:

Neighbor_list = hoomd.md.nlist.Cell(buffer=0.2,exclusions=[]) 

sigma = 1 
epsilon = 1
rmin_alpha=0.1
rcrit_alpha = 2

lj = hoomd.md.pair.LJ(nlist=Neighbor_list, default_r_cut=rcrit_alpha)
lj.params[("A", "A")] = {"sigma": 1.0, "epsilon": 1.0}

fenewca = hoomd.md.bond.FENEWCA()
fenewca.params['A'] = dict(k=kf, r0=R0, epsilon=epsilon, sigma=sigma, delta=0.0)

stepsize= 5e-4
integrator = hoomd.md.Integrator(dt=stepsize,integrate_rotational_dof=True)
simulation.operations.integrator = integrator # assign Integrator

# --- Forces --- #
integrator.forces.append(fenewca)
integrator.forces.append(lj)

Here is the suggested memray Flamegraph with my Custom potential.
I hope that by this everything will be clear.

memray-flamegraph-Microgel_simulation.py.16475.zip

[joaander]

As you can see in that flame graph, most of the memory is allocated by CellList (nlist.Cell in Python). The cell list allocates memory proportional to the volume of the simulation box, inversely proportional to r_cut cubed, and proportional to the local density of the densest region of your simulation. There is a note in the documentation to warn users of this behavior.

You don't provide a complete minimal working example, so I cannot tell you which of those factors is driving the memory usage. You can change your simulation parameters to mitigate any or all of those factors, or you can use nlist.Tree.

[Soulchemist97]

I do not think that nlist.tree would be a good idea, because I want to scale the number of particles in further simulations with the same volume.

Here is a minimal working example.
I do not understand why it is that way, even when my custom interaction is 0.
I have this problem not only for this structure but also for random particles I put into the simulation box.

import gsd.hoomd
import hoomd 
import numpy as np
from numpy import exp,pi,cos

simulation = hoomd.Simulation(device=hoomd.device.CPU(), seed=0)
simulation.create_state_from_gsd(filename="Gel_Particle.gsd")

Neighbor_list = hoomd.md.nlist.Cell(buffer=0.2,exclusions=[]) 

sigma = 1 
epsilon = 1
R0 = 1.5*sigma #maximal Bond extension
kf = 30

def V_alpha(r,alpha,epsilon=epsilon,sigma=sigma,R0=1.5):
    gamma = pi/(2.25-2**(1/3))
    beta = 2*pi-2.25*gamma
    V_Alpha = 0.5*alpha*epsilon*(cos(gamma*(r/sigma)**2+beta)-1)
    V_Alpha = np.where(r <= 2**(1/6)*sigma, V_Alpha, -epsilon * alpha)
    V_Alpha = np.where(R0*sigma <= r, V_Alpha, 0)
    return V_Alpha

rmin_alpha=0.1
rcrit_alpha = 2
alpha=0

r_alpha = np.linspace(start=rmin_alpha,stop=rcrit_alpha,num=500)
U_alpha = V_alpha(r_alpha,alpha)
F_alpha = -np.gradient(U_alpha,r_alpha) #Force as derivative

Solvophobicity_potential = hoomd.md.pair.Table(nlist=Neighbor_list,default_r_cut=rcrit_alpha)
Solvophobicity_potential.params["A","A"] = dict(r_min=rmin_alpha,U=U_alpha,F=F_alpha)

fenewca = hoomd.md.bond.FENEWCA()
fenewca.params['A'] = dict(k=kf, r0=R0, epsilon=epsilon, sigma=sigma,
                             delta=0.0)

stepsize= 5e-4
integrator = hoomd.md.Integrator(dt=stepsize,integrate_rotational_dof=True)
simulation.operations.integrator = integrator # Assign Integrator

# --- Forces --- #
integrator.forces.append(fenewca)
integrator.forces.append(Solvophobicity_potential)

# --- Thermostat ---#
nose_hoover = hoomd.md.methods.thermostats.MTTK(kT=1,tau=100*stepsize)
nvt = hoomd.md.methods.ConstantVolume(
    filter=hoomd.filter.All(),
    thermostat=nose_hoover,
)
nvt.anisotropic = True
integrator.methods.append(nvt) 

gsd_trajectory_writer = hoomd.write.GSD(
    filename="Trajectory.gsd",trigger=hoomd.trigger.Periodic(int(1e4)), 
    mode="wb", filter=hoomd.filter.All(),logger=None)

gsd_trajectory_writer.write_diameter = True
simulation.operations.writers.append(gsd_trajectory_writer)

simulation.run(int(1e6))

[joaander]

It is not a minimum working example if others cannot execute it. I get RuntimeError: GSD: No such file or directory - Gel_Particle.gsd when I run your script.

You don't even mention your simulation box dimensions.

If you don't want to use nlist.Tree, that is your choice. It will solve your problem. You can instead execute your simulations on a system with more RAM. Or you can reduce the size of your simulation box. There really isn't anything more that anyone can do to help you. I already explained all the variables that influence the memory usage of nlist.Cell.

[Soulchemist97]

Ok nlist.Tree fixed my problem.
I do not really understand why the memory scales for nlist.cell by the number of timesteps, but when it works now I am fine with it.
I could not share the gsd file I used for the code.
If this is a more common problem, I would try to reproduce my problem with a really minimum working example, but when not I will go with nlist.Tree and am really Thankful for your help