Re: Lipid membrane simulation

[Peter Košovan]

Dear Ksenia,

The "bond broken" error message means that a bond between two particles has been stretched beyond its maximum possible extension. Normally, this happens if a particle experiences a huge force in one integration step and consequently a huge displacement in the next step.

Your python script seems to be based on some outdated TCL script. Back then, the force cap was applied only to short-range interactions but not to bonded interactions. However, the definition of force cap has changed a couple of years ago. Now, the cap is applied to all interactions, including the bonded ones. Consequently, if you use a force cap and have bonded interactions in your system, then your bonds easily break when a rather small force is applied to them.

Recommended solution: do not use the force cap to relax your system if it contains bonded interactions. Instead, you can use the steepest descent energy minimization in order to avoid big overlaps after the initial setup of your system.

Although it's been a couple of years since the redefinition of force cap, I have not encountered a use case in which the new definition would be beneficial. Does anybody have an example of such a use case?

With regards,

Peter

On Thu, Nov 9, 2023 at 7:49 AM Ксения Астахова <astakhova.ksen.762@gmail.com> wrote:

Hello everyone!
I try to simulate lipids self-assembly with this tutorial https://espressomd.org/html/ess2013/Day2/05-1-mbtools.pdf. This article uses tcl and mbtools package but I want to use Python instead. So my code for simulation looks like:

import espressomd

import espressomd.observables

import espressomd.accumulators

import espressomd.analyze

import espressomd.interactions

import numpy as np

required_features = ["DPD", "LENNARD_JONES", "LJCOS2"]

espressomd.assert_features(required_features)

mol_num = 320

box_l = [14.0, 14.0, 14.0]

system = espressomd.System(box_l=box_l)

# Set bonded interactions

fene = espressomd.interactions.FeneBond(k=30, d_r_max=1.5)

hb = espressomd.interactions.HarmonicBond(k=10.0, r_0=4.0)

# bonded_parms = [fene, hb]

system.bonded_inter.add(fene)

system.bonded_inter.add(hb)

# Set non-bonded interactions

lj_eps = 1.0

lj_sigmah = 0.95

lj_sigma = 1.0

lj_offset = 0.0

system.non_bonded_inter[0, 0].lennard_jones.set_params(

epsilon=lj_eps, sigma=lj_sigmah, cutoff=1.1225 * lj_sigmah,

shift=0.25 * lj_eps, offset=lj_offset)

system.non_bonded_inter[0, 1].lennard_jones.set_params(

epsilon=lj_eps, sigma=lj_sigmah, cutoff=1.1225 * lj_sigmah,

shift=0.25 * lj_eps, offset=lj_offset)

system.non_bonded_inter[1, 1].lennard_jones_cos2.set_params(

epsilon=lj_eps, sigma=lj_sigma, offset=lj_offset, width=1.6)

particle_types = [0, 1, 1]

bond_l = 1.0

for i in range(mol_num):

tail_pos = np.random.random(3) * system.box_l

orient = 2 * np.random.random(3) - 1

orient /= np.linalg.norm(orient)

for j in range(len(particle_types)):

cur_part_id = i * len(particle_types) + j

particle_position = tail_pos + (len(particle_types) - j - 1) * bond_l * orient

system.part.add(id=cur_part_id, type=particle_types[j], pos=particle_position, rotation=(True, True, True))

if j > 0:

system.part.by_id(cur_part_id - 1).add_bond((fene, cur_part_id))

if j > 1:

system.part.by_id(cur_part_id - 2).add_bond((hb, cur_part_id))

# Integration parameters

langevin_gamma = 1.0 # not inited

main_time_step = 0.01

verlet_skin = 0.4

systemtemp = 1.1

dpd_gamma = 1.0

dpd_r_cut = 1.12 + 1.8

int_steps = 200

int_n_times = 1000

def warmup(steps, times, startcap = 5, capgoal = 1000, min_dist = 0):

capincr = (capgoal - startcap) / times

cap = startcap

for i in range(times):

act_min_dist = system.analysis.min_dist()

if act_min_dist < min_dist:

break

system.force_cap = cap

system.integrator.run(steps)

cap += capincr

system.force_cap = 0

# Warmup parameters

warm_time_step = 0.002

warmup_temp = 0

warmsteps = 100

warmtimes = 20

free_warmsteps = 0

free_warmtimes = 1

# Warmup

if warmup_temp == 0:

warmup_temp = systemtemp

system.time_step = warm_time_step

system.cell_system.skin = verlet_skin

system.thermostat.set_langevin(kT=warmup_temp, gamma=langevin_gamma, seed=42)

warmup(warmsteps, warmtimes)

# Second warmup

system.time_step = main_time_step

system.thermostat.set_langevin(kT=systemtemp, gamma=langevin_gamma, seed=42)

warmup(free_warmsteps, free_warmtimes, 1000)

system.time = 0

# Integration

system.thermostat.turn_off()

system.thermostat.set_dpd(kT=systemtemp, seed=41)

system.non_bonded_inter[0, 0].dpd.set_params(gamma=dpd_gamma, r_cut=dpd_r_cut)

system.non_bonded_inter[0, 1].dpd.set_params(gamma=dpd_gamma, r_cut=dpd_r_cut)

system.non_bonded_inter[1, 1].dpd.set_params(gamma=dpd_gamma, r_cut=dpd_r_cut)

for i in range(int_n_times):

print("\rrun %d at time=%.0f " % (i, system.time), end='')

system.integrator.run(int_steps)

But I receive an error during the first warmup: while calling method integrate(): ERROR: bond broken between particles 204, 205

I can't understand what I am doing wrong :( Maybe I have missed some details that changed since 3.x versions.

--

Peter Košovan

Associate professor

Department of Physical and Macromolecular Chemistry

Faculty of Science, Charles University, Prague, Czechia

Katedra fyzikální a makromolekulární chemie

Přírodovědecká fakulta Univerzity Karlovy v Praze

www.natur.cuni.cz/chemistry/fyzchem/
Tel. +420221951029

Fax +420224919752

We are constantly searching for talented PhD candidates and postdocs.