5- and 6-body interactions

Dear LAMMPS Users,

Could you please point me into the right direction which parts of LAMMPS I would have to modify to allow 5- and 6-body interactions. As far as I could see the maximum number of partners in any interaction is 4 for dihedral interactions. Here is some background information.

We developed a coarse grained model for double stranded DNA, which contains dihedral interactions on the single strand. We would like to deactivate some dihedral interactions if the other strand is close by.

The interaction we have in mind is a modified dihedral interaction between 4 atoms, which is informed about the proximity of 1 or 2 other atoms that are defined during setup.

Best regards and many thanks in advance,
Oliver

Dear LAMMPS Users,

Could you please point me into the right direction which parts of LAMMPS I
would have to modify to allow 5- and 6-body interactions. As far as I could
see the maximum number of partners in any interaction is 4 for dihedral
interactions. Here is some background information.

We developed a coarse grained model for double stranded DNA, which contains
dihedral interactions on the single strand. We would like to deactivate some
dihedral interactions if the other strand is close by.

The interaction we have in mind is a modified dihedral interaction between 4
atoms, which is informed about the proximity of 1 or 2 other atoms that are
defined during setup.

modifying the core of LAMMPS to include such kinds of interactions in
a similar fashion to angles and dihedrals would be a major
undertaking.

instead, you are probably best off following the approach that is used
in the (not yet released) support for CHARMM style CMAP interactions,
and implement your special interactions as a fix instead. with current
versions of LAMMPS, all the necessary infrastructure should be in
place (e.g. have a fix read a line in the header and a section of the
data file for topology and/or force field information and have it
write and read its state to and from a restart file).

steve or i can provide more technical details, if needed. if you want
to move the discussion of the mailing list, i would suggest submitting
an issue on github, Issues · lammps/lammps · GitHub and then
the discussion continue there.

axel.

Do you mean that the 5,6 atoms are on

another molecule (strand)? In that case, the CMAP stuff

is not so relevant as it involves dihedral/dihedral

interactions bewteen two overlapping dihedrals

on the same chain (a 5-body interaction).

I think you would need to write a new dihredral

style (or a fix style) that uses a dihedral but

also a neighbor list for the atoms, which

can look at the non-bond neighbors of atoms

in the dihedral to determine which are close

enough to influence the dihedral computation.

If the 2 strands are one molecule (with topological

bonds between them), you might just need to

look at other bond/angle partners of atoms

in the dihedral to find the “influence” atoms.

In that case, it seems you could do that in a new

dihedral style.

Steve

Would a distance based bond order method provide this interaction? It doesn’t really appear you want a 5 and 6 body interaction. You appear to be looking for the interaction of atoms based upon proximity to a particular set of atoms. Is that correct?

Jim

James Kress Ph.D., President

The KressWorks® Institute

An IRS Approved 501 ©(3) Charitable, Nonprofit Corporation

Engineering The Cure” ©

(248) 573-5499

Learn More and Donate At:

Website: http://www.kressworks.org

Facebook: https://www.facebook.com/KressWorks.Institute/

Twitter: @KressWorksFnd

Confidentiality Notice | This e-mail message, including any attachments, is for the sole use of the intended recipient(s) and may contain confidential or proprietary information. Any unauthorized review, use, disclosure or distribution is prohibited. If you are not the intended recipient, immediately contact the sender by reply e-mail and destroy all copies of the original message.

Hi Steve, Alex and Jim,

Many thanks for your quick response.

Atoms 5 and 6 are always on another strand. So it appears the CMAP thing is not relevant.

1 - 2 - 3 - 4 strand 1
      > >
   - 5 - 6 - strand 2

I think you are right that this wouldn’t be a true 5- or 6-body interaction. It would suffice to look at other bond partners of atoms in the dihedral as the modified dihedral would be based upon the proximity of a particular set of other atoms. More specifically, atom 2 in the dihedral (1-2-3-4) would have to identify its partner 5 in a harmonic bond (2-5).

Is there an efficient way to identify atom 5 from inside the dihedral interaction, which doesn’t involve walking through the entire list of harmonic bonds searching for the two partners (2-5)?

Cheers,
Oliver

Hi Steve, Alex and Jim,

Many thanks for your quick response.

Atoms 5 and 6 are always on another strand. So it appears the CMAP thing is not relevant.

there is a misunderstanding here. i didn't mean, that this would be a
CMAP type interaction, but that fixes like fix cmap are a mechanism in
LAMMPS to add custom interactions of *fixed* tuples of atoms, that do
not have to follow the available patterns like bond, angle, dihedral
or improper. at the moment, fix cmap is the only fix of that kind,
that i know.

1 - 2 - 3 - 4 strand 1
      > >
   - 5 - 6 - strand 2

I think you are right that this wouldn’t be a true 5- or 6-body interaction. It would suffice to look at other bond partners of atoms in the dihedral as the modified dihedral would be based upon the proximity of a particular set of other atoms. More specifically, atom 2 in the dihedral (1-2-3-4) would have to identify its partner 5 in a harmonic bond (2-5).

Is there an efficient way to identify atom 5 from inside the dihedral interaction, which doesn’t involve walking through the entire list of harmonic bonds searching for the two partners (2-5)?

yes, but there are some complications. while bonds are stored with
atoms, only with the settings "newton on off" or "newton off off", the
2-5 and 3-6 bonds would be stored with both, atom 2 and 5 or 3 and 6,
respectively. and then you can loop over the bonds of an atom like
this:

for (i = 0; atom->num_bonds[i2]; ++i) {
   // atom->bond_atom[i2][i] will contain the atom ids of the atoms
bonded to atom 2
   [...]
}

without force->newton_bond == 1, things are more complicated. you
would have to make certain, that the 2-5 bond would be stored with
atom 2, and then you would be unable to also have the inverse
interaction on strand 2, as shown below.

       2 - 3 strand 1
       > >
  7 - 5 - 6 - 8 strand 2

so, if i understand this right, you have three options: 1) write a
custom dihedral style and *enforce* that force->newton_bond is off or
2) implement this interaction as a fix, where you simple would
explicitly define the tuples: 1 2 3 4 5 6 (and correspondingly 7 5 6
8 2 3 and so on) and then make this fix compute and apply the distance
based correction to the "normal" dihedral interactions for 1-2-3-4 and
7-5-6-8. this would definitely be the most efficient process. option
3) would be to request a neighbor list with a suitable cutoff and then
searching for the 2-5 and 3-6 bonds on either atoms 2 and 3 or their
neighbors.

please note, walking the bondlist, won't work, since that contains
only "owned" bonds, and for bonds between local and ghost atoms, i.e.
atoms owned by different subdomains, you would see the bond only half
the time, when force->newton_bond is 1 (as is the default). in that
case, you need to enforce force->newton_bond 0 and then you can just
as well follow option 1).

axel.