Coding a new fix for a triangulated membrane model

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I’m trying to simulate a membrane model using LAMMPS.

This triangulated membrane model consists of beads connected through hard bonds. Hard repulsion at close range is employed to prevent bead overlapping. There is a harmonic dihedral associated with the two neighboring triangles of each bond to simulate the bending energy of the membrane. The above three are all the interactions in this membrane model.

What’s special is, to add fluidity into this model, bond flipping needs to be performed. Please have a look at the attached image, which demonstrates one bond flipping. Each bond flipping requires changing the data of one bond and several related dihedrals. Every N time steps, the entire membrane needs to be scanned, and all eligible bonds need to be flipped.

This model is described in this paper: http://iopscience.iop.org/0295-5075/12/4/002

I suppose LAMMPS can’t do the above bond flipping, right? So I’m thinking about coding a new fix for this. I have some experience in C/C++ but I’ve never modified LAMMPS. I read the manual, developer guide, and dived into the source code. I think I’m clear about the data structure and how to change bonds, dihedrals. However, to flip a bond near a processor boundary, several processors are involved and have to work together, because a processor can only change its owned bonds, dihedrals (correct?), and a bond flipping requires changing all neighboring dihedrals that can possibly cross the processor boundary. So this requires communication among processors, which I don’t quite understand. I vaguely know I need to implement the pack_comm, unpack_comm, pack_reverse_comm, unpack_reverse_comm methods, and possibly pack_exchange and unpack_exchange?

What’s the difference between forward communication and reverse communication?

I’m really lost in the complexity of LAMMPS. Maybe what I have thought so far is entirely garbage. Please point me to the right direction.

Any help, suggestion, advice, insight is welcome.

Thanks!

Sincerely,

example.png

2

Have you looked at bond/swap, bond/break, and bond/create commands? I suppose bond/swap is not something you want, but it seems to me that what you wish to achieve can be done with a combination of bond/break and bond/create commands.

Forward communication copies information of owned local atoms to ghost atoms, reverse communication does the opposite.

Ray

3

There are two more options to look at: the BODY package and manybody pair styles like airebo, they compute the topology and bond, angle and dihedral contributions on the fly from post processed neighbor lists.

Axel

4

Hi,

Thanks for you guys’ help. Really appreciate it! I read the code you guys suggested, and have more understanding of LAMMPS. I still have a few questions. Sorry this may be a bit tedious but it’s really not something I can do alone.

- Bond flipping criteria

My bond flipping criteria are:

  1. The new bond length will be between 0.7 and 1.3.
  2. The two atoms forming the old bond both have more than 2 bonds before the flipping. (This is to prevent dangling atoms.)
  3. The new bond isn’t already there before the flipping.

So as you see those are a bit complicated and I don’t think it can be done using bond/break and bond/create commands, right?

And for now, I’ll just try to implement non-parallel flipping.

- Forward and reverse communication

Is the following understand of forward and reverse communication correct?

  1. If a processor makes changes to its owned atoms, it should initiate a forward communication to reflect the changes to ghost atoms on other processors.
  2. If a processor makes changes to its ghosts atoms, it should initiate a reverse communication to reflect the changes back to owned atoms.

I saw a possibility of race conditions. One atom can be ghost atoms on two processors. What if the two processors made different changes to the two ghost atoms referring to the very same real atom and called a reverse communication?

- Calculating distance in periodic boundary condition

As you see in the bond flipping criteria, I need to calculate the distance between the two atoms that are going to form a new bond. Because I do use PBC to achieve a virtually infinite membrane, I need to be careful since one of the atom can be on the other side of the box, a ghost atom. In pair_lj_cut.cpp and fix_bond_create.cpp, I saw the distance was calculated by directly calculating the difference of coordinates, which didn’t seem to take into account of PBC. Could it be that the coordinates of the ghost atoms are already prepared so that the distance can be calculated that way?

Again, thanks guys! I can’t be more grateful.

Sincerely,

– Cong, the anxious Ph.D. student :slight_smile:

5

Hi,

Thanks for you guys' help. Really appreciate it! I read the code you guys
suggested, and have more understanding of LAMMPS. I still have a few
questions. Sorry this may be a bit tedious but it's really not something I
can do alone.

*- Bond flipping criteria*

My bond flipping criteria are:

   1. The new bond length will be between 0.7 and 1.3.

Fix bond/create and bond/break have Rmin and Rmax keywords.

   1. The two atoms forming the old bond both have more than 2 bonds
   before the flipping. (This is to prevent dangling atoms.)

This does not matter.

   1. The new bond isn't already there before the flipping.

Of course, hence the "create".

So as you see those are a bit complicated and I don't think it can be done

using bond/break and bond/create commands, right?

As such, I don't see why you can't achieve what you want to achieve with a
combination of bond/break and bond/create.

I suggest you, instead of speculating, try to realize what you want to
achieve with bond/break and bond/create, then post your questions and input
files.

Of course, please feel free to implement the new "bond/flip" if that is
easier.

And for now, I'll just try to implement non-parallel flipping.

*- Forward and reverse communication*

Is the following understand of forward and reverse communication correct?

   1. If a processor makes changes to its owned atoms, it should initiate
   a forward communication to reflect the changes to ghost atoms on other
   processors.
   2. If a processor makes changes to its ghosts atoms, it should
   initiate a reverse communication to reflect the changes back to owned atoms.

Correct, but these changes are summed up eventually so it does not matter

which processor adds first.

I saw a possibility of race conditions. One atom can be ghost atoms on two
processors. What if the two processors made different changes to the two
ghost atoms referring to the very same real atom and called a reverse
communication?

*- Calculating distance in periodic boundary condition*

As you see in the bond flipping criteria, I need to calculate the distance
between the two atoms that are going to form a new bond. Because I do use
PBC to achieve a virtually infinite membrane, I need to be careful since
one of the atom can be on the other side of the box, a ghost atom. In
pair_lj_cut.cpp and fix_bond_create.cpp, I saw the distance was
calculated by directly calculating the difference of coordinates, which
didn't seem to take into account of PBC. Could it be that the coordinates
of the ghost atoms are already prepared so that the distance can be
calculated that way?

Yes, ghost atoms are PBC images of real atoms.

Ray

6

Hi,

Thanks for you guys' help. Really appreciate it! I read the code you guys
suggested, and have more understanding of LAMMPS. I still have a few
questions. Sorry this may be a bit tedious but it's really not something I
can do alone.

it only seems too complicated to you, because you want to solve the
entire problem in one big step. that almost never works and is indeed
tedious and intimidating.

why don't you just start small, build a test system with a few atoms
and try out using fix bond/create and bond/break and see how far this
will take you. the logic to correctly do all the things you need to do
in parallel correctly is quite tricky and over the last few months
different people spent a significant amount of time to track this
down. if i would be in your position, i would not assume that i would
get this right immediately, so i'd rather try and reuse what is
already known to work and then see, if this leads to some limitations
that require me to build a new fix styles.

and - again - let me remind you, that from my current understanding,
it may be easier to assemble the connectivity data on the fly and
implement your CG membrane model as a pair style only. the fact that
you build a network of bonds with dihedral interactions along the bond
topology doesn't mean that it must be implemented using explicit
bonds/dihedrals in LAMMPS. those are mostly meant for cases where they
don't change. changing the settings and exclusions and so on is fairly
inefficient and thus should only be used where this would not be time
critical.

- Bond flipping criteria

My bond flipping criteria are:

The new bond length will be between 0.7 and 1.3.
The two atoms forming the old bond both have more than 2 bonds before the
flipping. (This is to prevent dangling atoms.)
The new bond isn't already there before the flipping.

So as you see those are a bit complicated and I don't think it can be done
using bond/break and bond/create commands, right?

i disagree. those fixes allow very complex constraints. i think all
that is primarily needed would be to create fix bond/break before fix
bond/create.

And for now, I'll just try to implement non-parallel flipping.

what is that?

- Forward and reverse communication

Is the following understand of forward and reverse communication correct?

If a processor makes changes to its owned atoms, it should initiate a
forward communication to reflect the changes to ghost atoms on other
processors.
If a processor makes changes to its ghosts atoms, it should initiate a
reverse communication to reflect the changes back to owned atoms.

I saw a possibility of race conditions. One atom can be ghost atoms on two
processors. What if the two processors made different changes to the two
ghost atoms referring to the very same real atom and called a reverse
communication?

those race conditions are avoided since only pairs of atoms where at
least one of the two atoms is a local atom can form a new bond and
then any atom can only be a local atom in one subdomain. this has been
sorted out to the best of my knowledge in fix bond/create.

- Calculating distance in periodic boundary condition

As you see in the bond flipping criteria, I need to calculate the distance
between the two atoms that are going to form a new bond. Because I do use
PBC to achieve a virtually infinite membrane, I need to be careful since one
of the atom can be on the other side of the box, a ghost atom. In
pair_lj_cut.cpp and fix_bond_create.cpp, I saw the distance was calculated
by directly calculating the difference of coordinates, which didn't seem to
take into account of PBC. Could it be that the coordinates of the ghost
atoms are already prepared so that the distance can be calculated that way?

yes. ghost atoms are *copies*. LAMMPS is *not* subject to minimum
image conventions. a single local atoms can have many ghost atoms
depending on the cutoff used for building the neighbor list.

axel.

7

well, here is one limitation to consider with the fix break and create commands. fix bond/create only really counts the number of bonds per atom at the beginning of the run in its setup. it then increments the count as it adds new bonds. that means that fix bond/create won’t understand that bonds have been broken by another fix (like bond/break) until the next run.