pairwise repulsion for excluded intra-molecular pairs

LAMMPS LAMMPS Mailing List Mirror
1

Hello everyone,

I am doing MD simulations of triclinic molecular crystal TATB, in which the non-bonded interactions are modeled using the Buckingham potential and Coloumb interactions with fixed partial charges. Bonded interactions are modeled using bonds, angles, dihedrals and impropers. All intra-molecular non-bonded interactions are excluded. Long range Coloumbic interactions are computed using the PPPM solver. The parameters are taken from J. Chem. Phys. 139, 074503 (2013). My LAMMPS version is Dec. 7, 2013.

There are 6 atom types in my molecule. I need to introduce an additional pairwise, repulsive interaction only between the intra-molecular 1-2 pairs. This repulsive interaction should be excluded for all inter-molecular 1-2 pairs. This is what I think I should do:

  1. Tabulate the r^12 repulsion with non-zero values only for the 1-2 interaction.

  2. Use pair_style hybrid/overlay buck/coul/long table

  3. use neigh_modify exclude molecule group36, where group36 is composed of all atoms of type 3-6

  4. modify pair_table.cpp with an ‘if loop’ to calculate the repulsion only for the intra-molecular 1-2 pairs.

I haven’t looked into the details of these steps to be sure that this will work. I was wondering if anyone can suggest an easier way, i.e. using available keywords/commands, to accomplish what I had described.

Thanks in advance for any thoughts/comments.

2

Hello everyone,

I am doing MD simulations of triclinic molecular crystal TATB, in which the
non-bonded interactions are modeled using the Buckingham potential and
Coloumb interactions with fixed partial charges. Bonded interactions are
modeled using bonds, angles, dihedrals and impropers. All intra-molecular
non-bonded interactions are excluded. Long range Coloumbic interactions are
computed using the PPPM solver. The parameters are taken from J. Chem. Phys.
139, 074503 (2013). My LAMMPS version is Dec. 7, 2013.

There are 6 atom types in my molecule. I need to introduce an additional
pairwise, repulsive interaction only between the intra-molecular 1-2 pairs.
This repulsive interaction should be excluded for all inter-molecular 1-2
pairs. This is what I think I should do:

1. Tabulate the r^12 repulsion with non-zero values only for the 1-2
interaction.
2. Use pair_style hybrid/overlay buck/coul/long table
3. use neigh_modify exclude molecule group36, where group36 is composed of
all atoms of type 3-6
4. modify pair_table.cpp with an 'if loop' to calculate the repulsion only
for the intra-molecular 1-2 pairs.

I'm not sure what an "if loop" is in this context. Can you clarify step 4?

I can think of two other ways:
1)
If you know exactly which pairs of atoms of type 1 and 2 are likely to
be close contact, or if there aren't too many of them, you can try
connecting all of these pairs of atomswith explicit bonds and use
bond_style table (or hybrid harmonic table, etc..) to add specific
1/r^12 interactions between these atom pairs.

2)
Alternately, there is a pair_style which allows you to use different
Lennard-Jones coefficients depending on whether the two atoms belong
to the same molecule. Check out:

http://www.moltemplate.org/lammps_code/pair_lj_charmm_coul_charmm_inter.html

You can download the code for that pair style here:
http://www.moltemplate.org/lammps_code/index.html

At least 3 different groups have used it successfully so far.
It was tested in late 2013. (Let me know if LAMMPS has changed and
this code no longer compiles.)

The only problem is that this pair_style assumes you want to use
short-range coulombics/electrostatics. If you want to use long-range
coulombics, you would have to turn off the short range electrostatics
forces, and then use
"pair_style hybrid/overlay" to combine these Lennard Jones forces with
"coul/long".

There's a (long and confusing) discussion of how to combine this style
with long-range coulombics here:
http://lammps.sandia.gov/threads/msg32218.html
Here's an excerpt from that post:

pair_style hybrid/overlay &
lj/charmm/coul/charmm/inter 29.5 30 0 0 &
coul/long 14.0
pair_coeff 1 1 lj/charmm/coul/charmm/inter 0 0 0 0 0.0298 3.4
pair_coeff 1 2 lj/charmm/coul/charmm/inter 0 0 0 0 0.0298 3.5
pair_coeff 2 2 lj/charmm/coul/charmm/inter 0 0 0 0 0.0298 3.6
pair_coeff * * coul/long

Cheers
Andrew

3

There is no need to tabulate a 1/r^12 interaction. Just
use a LJ potential with a short cutoff and you will get only
the repulsive part. I would see if you could formulate

your model with a single pair style, not hybrid, if possible.

I suppose a neigh_modify exclude option could be added
for your case. It already allows for the inverse of what
you want (exclude intra, include inter).

Steve

4

Thanks (!), Andrew and Steve for the comments.

What I meant by the ‘if loop’ is to add an additional check in the force calculations to see if the particular 1-2 pair is an inter or an intra-molecular pair. I haven’t looked enough to see if this can be done easily.

I think I can use Andrew’s suggestion of using a bond_style table. In my case, the intra-molecular 1-2 pairs for which the repulsion has to be introduced is defined by the geometry of the molecule. So I can define an additional bond type and specify the pairs without much trouble. This way, I can also avoid multiple walks through the neighbor list which comes with the hybrid pair_style.

5

Yep.

If you're willing to edit the lammps source code, then that opens up
options for you. Before you calculate the force between atoms i and
j, include this check:

  if (atom->molecule[i] == atom->molecule[j])

Try this with a pair_style of your choice (one that uses long-range
coulombics for example).
(That's what "pair_style lj/charmm/coul/charmm/inter" does. See line
199 of pair_lj_charmm_coul_charmm_inter.cpp. Attached. The complete
code is currently downloadable at http://moltemplate.org/lammps_code/)

As for getting pure 1/r^12 repulsion, (with no 1/r^6 term), you can
approximate this by using a small epsilon and a larger sigma (and then
use a cutoff too).

  # -------------- scratchwork --------------

pair_lj_charmm_coul_charmm_inter.cpp (38.9 KB)

6

Thanks, Andrew. Setting up an extra bond type appears at this point like a safer (than tinkering with the source code), easier and faster (wrt run-time) option. I’ll try that first.