Out of memory error

[dpo854]

Hi everyone, I made some modifications to the grand canonical code to simulate ions under cylindrical confinement. However, it throws up the following error.

/var/spool/slurmd/job6198668/slurm_script: line 19: 29195 Killed $ESPRESSO_PATH/pypresso MD-cell.py 0.0000048184 -0.06561 slurmstepd: error: Detected 1 oom-kill event(s) in step 6198668.batch cgroup. Some of your processes may have been killed by the cgroup out-of-memory handler.

All I have done is put a cylindrical wall in the center of the box. The corresponding lines are given below:

cylinder = shapes.Cylinder(center=[boxl/2.0, boxl/2.0, boxl/2.0], axis=[0, 0, 1],direction=-1,radius=boxl/6.0,length=boxl/2.0)
myConstraint = system.constraints.add(shape=cylinder, particle_type=4)

[KaiSzuttor]

Could you please check if you script works without a scheduling system?

[fweik]

I made the following minimal working script from the code you posted:

import espressomd
from espressomd import shapes, System

boxl=10.
system = System(box_l=3*[boxl])

cylinder = shapes.Cylinder(center=[boxl/2.0, boxl/2.0, boxl/2.0], axis=[0, 0, 1],direction=-1,radius=boxl/6.0,length=boxl/2.0)
myConstraint = system.constraints.add(shape=cylinder, particle_type=4)

according to valgrind's massif tool, this has a peak memory usage of 8.406 MB, which sounds reasonable. It could be possible that you have to request memory with your submission script to slurm.

[dpo854]

Requesting memory solved the problem for the time being. However, another problem came up, which unless resolved can make it difficult to predict whether the current problem has been resolved. The error is due probably to large forces experienced by particles, but I had made sure to put all particles inside the cylinder. This actually causes the electrostatic solver to fail as shown below.

ERROR: Constraint violated by particle

File "MD-cell.py", line 113, in system.actors.add(p3m) File "espressomd/actors.pyx", line 215, in espressomd.actors.Actors.add (/home/dpo854/mycode/espresso/build/src/python/espressomd/actors.cpp:5586) File "espressomd/actors.pyx", line 70, in espressomd.actors.Actor._activate (/home/dpo854/mycode/espresso/build/src/python/espressomd/actors.cpp:2365) File "espressomd/electrostatics.pyx", line 343, in espressomd.electrostatics.P3M._activate_method (/home/dpo854/mycode/espresso/build/src/python/espressomd/electrostatics.cpp:6925) File "espressomd/electrostatics.pyx", line 334, in espressomd.electrostatics.P3M._tune (/home/dpo854/mycode/espresso/build/src/python/espressomd/electrostatics.cpp:6664) File "espressomd/electrostatics.pxd", line 136, in espressomd.electrostatics.python_p3m_adaptive_tune (/home/dpo854/mycode/espresso/build/src/python/espressomd/electrostatics.cpp:13360) File "espressomd/utils.pyx", line 264, in espressomd.utils.handle_errors (/home/dpo854/mycode/espresso/build/src/python/espressomd/utils.cpp:7717) Exception: Error in p3m_adaptive_tune: ERROR: Constraint violated by particle 33

[dpo854]

I am sorry for rushing with the program. Later thinking a bit more deeply about the particle addition/deletion scheme in espresso, I realized that perhaps the constraint due to the wall does not automatically limit the positions of the created particles in grand canonical ensemble to the inside of the cylinder, which is why they might be colliding with the wall and hence, getting kicked out of the box. Is there a way to constrain the reaction to occur only in the inner part of the cylinder?

[jonaslandsgesell]

I am sorry for rushing with the program. Later thinking a bit more deeply about the particle addition/deletion scheme in espresso, I realized that perhaps the constraint due to the wall does not automatically limit the positions of the created particles in grand canonical ensemble to the inside of the cylinder, which is why they might be colliding with the wall and hence, getting kicked out of the box.

The situation you mention will happen, unless you restrict the volume to which you want to insert the ions. This is because it is not clear what the default behavior should be. The program cannot know which system should be simulated (i.e. whether the region of interest is inside or outside the cylindrical constraint, or both?).

A cylindrical constraint in z-direction can be set for the reaction_ensemble module. With the following command you restrict proposed insertions to a cylinder centered around(center_x, center_y) with radius radius_of_cylinder :

RE = reaction_ensemble.ReactionEnsemble(
        temperature=1,
        exclusion_radius=1,
        seed=77)
RE.set_cylindrical_constraint_in_z_direction(center_x, center_y, radius_of_cylinder)
RE.set_volume(volume_of_cylinder)

There is a similar constraint for a slab-like situation called set_wall_constraints_in_z_direction. Please note that you also need to provide the volume of your region of interest via set_volume, otherwise you will encounter wrong results: the reaction ensemble and the grand canonical ensemble are sensitive to the volume V.

In order to avoid the "explosion" of the simulation due to too high forces, you need to prevent particle insertions too close to the wall. Therefore, you will have to choose radius_of_cylinder about 0.9 sigma smaller than the WCA-cylinder wall which is set for the MD particles. Otherwise, you will get too high forces, due to insertions too close to the wall.

[KaiSzuttor]

btw why isnt there just a “respect constraints” flag?

[dpo854]

Hi Jonas, I added the above lines to my code and put a gap of at least 2 sigma between radius_of_cylinder and WCA cylinder wall. Unfortunately, the error still remains.

[dpo854]

The code is attached below. It's the same code as grandcanonical.py except for the lines on the constraint. I have tried to simulate this cell model with hoomd, which gave me memory error.

"""
Perform a grand canonical simulation of a system in contact
with a salt reservoir while maintaining a constant chemical potential.
"""
epilog = """
Takes two command line arguments as input: 1) the reservoir salt
concentration in units of 1/sigma^3 and 2) the excess chemical
potential of the reservoir in units of kT.
The excess chemical potential of the reservoir needs to be determined
prior to running the grand canonical simulation using the script called
widom_insertion.py which simulates a part of the reservoir at the
prescribed salt concentration. Be aware that the reservoir excess
chemical potential depends on all interactions in the reservoir system.
"""
import numpy as np
import argparse

import espressomd
from espressomd import reaction_ensemble
from espressomd import electrostatics
from espressomd import shapes
from espressomd import constraints
from espressomd.io.writer import vtf

required_features = ["P3M", "EXTERNAL_FORCES", "WCA"]
espressomd.assert_features(required_features)

parser = argparse.ArgumentParser(epilog=__doc__ + epilog)
parser.add_argument('cs_bulk', type=float,
                    help="bulk salt concentration [1/sigma^3]")
parser.add_argument('excess_chemical_potential', type=float,
                    help="excess chemical potential [kT] "
                         "(obtained from Widom's insertion method)")
args = parser.parse_args()

# System parameters
#############################################################

cs_bulk = args.cs_bulk
excess_chemical_potential_pair = args.excess_chemical_potential
box_l = 100.0
boxl=100.0
print("csbulk",cs_bulk,"excess", excess_chemical_potential_pair)
# Integration parameters
#############################################################
system = espressomd.System(box_l=[box_l, box_l, box_l])
np.random.seed(seed=42)

system.time_step = 0.01
system.cell_system.skin = 0.4
temperature = 2.4935

print(system.box_l)
#############################################################
#  Setup System                                             #
#############################################################

# Particle setup
#############################################################
# type 0 = HA
# type 1 = A-
# type 2 = H+
# type 3 = MONO

for i in range(int(cs_bulk * box_l**3)):
    system.part.add(pos=system.box_l/2.0+(0.5-np.random.random(3)) * system.box_l/3.0, type=1, q=-2.9761)
    system.part.add(pos=system.box_l/2.0+(0.5-np.random.random(3)) * system.box_l/3.0, type=2, q=2.9761)
    print("inside loop",system.box_l/2.0,(0.5-np.random.random(3)) * system.box_l/3.0) 
wca_eps = 2.4935
wca_sig = 1.0
types = [0, 1, 2, 3,4]
for type_1 in [0,1,2,3]:
    for type_2 in types:
        system.non_bonded_inter[type_1, type_2].wca.set_params(
            epsilon=wca_eps, sigma=wca_sig)

RE = reaction_ensemble.ReactionEnsemble(
    temperature=temperature, exclusion_radius=1.0, seed=3)
RE.add_reaction(
    gamma=cs_bulk**2 * np.exp(excess_chemical_potential_pair / temperature),
    reactant_types=[], reactant_coefficients=[], product_types=[1, 2],
    product_coefficients=[1, 1], default_charges={1: -1, 2: +1})
print(RE.get_status())
system.setup_type_map([0, 1, 2])

cylinder = shapes.Cylinder(center=[boxl/2.0, boxl/2.0, boxl/2.0], axis=[0, 0, 1],direction=-1,radius=boxl/3.0,length=boxl/2.0)
myConstraint = system.constraints.add(shape=cylinder, particle_type=4)
RE.set_cylindrical_constraint_in_z_direction(boxl/2.0, boxl/2.0, boxl/3.0-2.0)
volume=3.14*((boxl/3.0-2.0)**2)*(boxl)
RE.set_volume(volume)
print("Hello")

RE.reaction(10000)

p3m = electrostatics.P3M(prefactor=2.0, accuracy=1e-3)
system.actors.add(p3m)
p3m_params = p3m.get_params()
for key, value in p3m_params.items():
    print("{} = {}".format(key, value))

# Warmup
#############################################################
# warmup integration (steepest descent)
warm_steps = 20
warm_n_times = 20
min_dist = 0.9 * wca_sig

# minimize energy using min_dist as the convergence criterion
system.integrator.set_steepest_descent(f_max=0, gamma=1e-3,
                                       max_displacement=0.01)
i = 0
while system.analysis.min_dist() < min_dist and i < warm_n_times:
    print("minimization: {:+.2e}".format(system.analysis.energy()["total"]))
    system.integrator.run(warm_steps)
    i += 1

print("minimization: {:+.2e}".format(system.analysis.energy()["total"]))
print()
system.integrator.set_vv()

# activate thermostat
system.thermostat.set_langevin(kT=temperature, gamma=.5, seed=42)

# MC warmup
RE.reaction(1000)
fp = open('trajectory.vtf', mode='w+t')
vtf.writevsf(system, fp)
vtf.writevcf(system, fp)

n_int_cycles = 10
n_int_steps = 60
num_As = []
deviation = None
for i in range(n_int_cycles):
    RE.reaction(10)
    system.integrator.run(steps=n_int_steps)
    num_As.append(system.number_of_particles(type=1))
    vtf.writevcf(system, fp)
    if i > 2 and i % 50 == 0:
        print("HA", system.number_of_particles(type=0), "A-",
              system.number_of_particles(type=1), "H+",
              system.number_of_particles(type=2))
        concentration_in_box = np.mean(num_As) / box_l**3
        deviation = (concentration_in_box - cs_bulk) / cs_bulk * 100
        print("average num A {:.1f} +/- {:.1f}, average concentration {:.8f}, "
              "deviation to target concentration {:.2f}%".format(
                  np.mean(num_As),
                  np.sqrt(np.var(num_As, ddof=1) / len(num_As)),
                  concentration_in_box, deviation))
fp.close()

[dpo854]

There is no way I can upload my file here. Please let me know if you need the file.

[jonaslandsgesell]

Hello Debadutta,

from a quick look at how the cylindrical constraint is setup, I guess the problem lies there. If you want to debug your script, please check the positions of all particles and check whether "unexpected" positions occur. Good candidates for possible blow ups of your simulation are also these parameters in the cylindrical constraint part:

• direction=-1 (please check)
• length=boxl/2.0 (better choose more, e.g. add 20 , otherwise the cylinder could have interacting caps)

I will send you my single-chain CGM simulation script today in the evening so you can perform the simulations.

[dpo854]

Thanks Jonas for offering to help. I have indeed set the direction to -1 at the time of defining the cylinder and made sure that all particles initially respect the constraint by explicitly printing out their positions while they are being created. Of course, this happens before RE.reaction(), meaning there is no way I could track the positions of added particles. Regarding the length, does it take the total length of the cylinder as the input or the half of total length. The reason why I ask this is in one of the illustrations of available shapes in your documentation, the length was shown to be half the actual length of the cylinder. Thanks, Debadutta Prusty

[jonaslandsgesell]

I sent you an email with my simulation setup of the single-chain CGM (also to @kosovan, @KaiSzuttor , @fweik ).

[jonaslandsgesell]

Can this be closed here? The wrongly documented parameters for the cylinder could have resulted in the explosion of the system.