Re: Espressomd-users Digest, Vol 150, Issue 3

[Ксения Астахова]

Hello everyone!

I have managed to make the warmup with steepest descent and integration is working now but the membrane has formed just 2 times out of 10 I suppose. I have also tried the dpd thermostat (commented code) but it doesn't help. Do you have any idea how to increase the probability of its forming? Or maybe someone has done lipid membrane forming and can share with me.

I use this article https://espressomd.org/html/ess2013/Day2/05-1-mbtools.pdf. It says that I can increase one dimension of the box but it doesn't seem to help me.

Best wishes,

Kseniia

чт, 9 нояб. 2023 г. в 13:23, <espressomd-users-request@nongnu.org>:

Send Espressomd-users mailing list submissions to
        espressomd-users@nongnu.org

To subscribe or unsubscribe via the World Wide Web, visit
        https://lists.nongnu.org/mailman/listinfo/espressomd-users
or, via email, send a message with subject or body 'help' to
        espressomd-users-request@nongnu.org

You can reach the person managing the list at
        espressomd-users-owner@nongnu.org

When replying, please edit your Subject line so it is more specific
than "Re: Contents of Espressomd-users digest..."


Today's Topics:

   1. Re: Lipid membrane simulation (Peter Košovan)


----------------------------------------------------------------------

Message: 1
Date: Thu, 9 Nov 2023 09:22:38 +0100
From: Peter Košovan <peter.kosovan@natur.cuni.cz>
To: Ксения Астахова  <astakhova.ksen.762@gmail.com>
Cc: espressomd-users@nongnu.org
Subject: Re: Lipid membrane simulation
Message-ID:
        <CAHG7JguHmbwBsYx4Rja2TF9_13HdLJ8J02WFCjQcmn54eWVdnQ@mail.gmail.com" target="_blank">CAHG7JguHmbwBsYx4Rja2TF9_13HdLJ8J02WFCjQcmn54eWVdnQ@mail.gmail.com>
Content-Type: text/plain; charset="utf-8"

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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.nongnu.org/archive/html/espressomd-users/attachments/20231109/3a741811/attachment.htm>

------------------------------

Subject: Digest Footer

_______________________________________________
Espressomd-users mailing list
Espressomd-users@nongnu.org
https://lists.nongnu.org/mailman/listinfo/espressomd-users


------------------------------

End of Espressomd-users Digest, Vol 150, Issue 3
************************************************

membrane-2.py
Description: Text Data