How to exclude a type of interaction while keeping another type of interaction on the same group of atoms

Dear Lammps users

I want to simulate a system containing graphite which includes 30 layers of graphene. For some reasons, I need to use Tersoff and LJ interaction styles simultaneously by using hybrid/overlay pair style on all atoms but exclude the LJ interactions on atoms within layers. I have defined 30 molecules as for 30 layers of graphene in order to use command “neigh_modify exclude molecule”. The problem is that I can not exclude LJ interactions while keeping Tersoff interactions within those layers.
Any help is appreciated.
Thanks in advance,
Hasan

There isn't a simple way to do this, that I can think of,
since the neigh_modify exclusions apply to both neighbor
lists (for Terssoff and LJ). Here are 2 options, neither ideal.

a) define 30 atom types, one per layer, the define
the epsilon for each I = J interaction to be 0.0. I != J is
normal non-zero LJ. However, this will require your Tersoff
potential to recognize 30 types, diffo for the Tersoff potential
file - which would be tedious to construct

b) use the molecule approach you describe, but don't use
neigh-modify exclude. Add some code to pair_lj_cut.cpp
to check if the molecule IDs of the 2 atoms are the same
and skip the interaction if so.

Steve

There isn't a simple way to do this, that I can think of,
since the neigh_modify exclusions apply to both neighbor
lists (for Terssoff and LJ). Here are 2 options, neither ideal.

a) define 30 atom types, one per layer, the define
the epsilon for each I = J interaction to be 0.0. I != J is
normal non-zero LJ. However, this will require your Tersoff
potential to recognize 30 types, diffo for the Tersoff potential
file - which would be tedious to construct

b) use the molecule approach you describe, but don't use
neigh-modify exclude. Add some code to pair_lj_cut.cpp
to check if the molecule IDs of the 2 atoms are the same
and skip the interaction if so.

If you only want to turn off the LJ interactions, then there's a new
pair style "lj/coul/charmm/coul/inter" which does what Steve
described. I attached the source code files to this message. Copy
these files to the "src" directory of your LAMMPS folder and recompile
LAMMPS, (I'll try to keep this pair style up to date. In the future
you should be able to download updated versions of this pair style
from moltemplate.org.) For now, the documentation and examples for
this pair style is located here:

In your case, since you are also using tersoff, I would try something like this:

pair_style hybrid/overlay tersoff lj/charmm/coul/charmm/inter 7.5 8.0
(This should all fit on one line).

The pair_coeff command which turns off LJ interactions between atoms
in the -same- molecule would be this:

pair_coeff 1 1 0.0 3.6 0.0 3.6 0.40 3.6

  --- explanation of the pair_coeff command ---
Following the atom type pair ("1 1" in this example),

-the next two columns are epsilon and sigma parameters for atoms in
the same molecule.

-the next two columns are the epsilon and sigma parameters for 1-4 interactions
   (which are usually in the same molecule).

-the last two columns are the epsilon and sigma parameters for atoms
in different molecules.

---- repulsive interactions ---

If you only want to turn off the attractive portion of the LJ
interaction, you could use:

pair_style hybrid/overlay tersoff lj/charmm/coul/charmm/inter es4k4l 7.5 8.0
    (again this should fit on one line)

pair_coeff 1 1 0.40 3.6 1 0 0.40 3.6 10 0.40 3.6 1 -1

     The extra parameters ("1 0" or "1 -1") are the dimensionless
coefficients which are inserted in front of the 1/r^12 and 1/r^6
terms. Mixing rules work for this pair style as well, so you don't
have to define coeffs for every atom pair. (See documentation in link
above.)

     Let me know if you have any trouble using this.

    If you want to create your own pair style, the relevant section of
the code are on lines 135 and following line 199 of
"pair_lj_charmm_coul_charmm_inter.cpp" (attached file).

int *molecule = atom->molecule;
:
and
:
if (molecule[i] == molecule[j]) {
  do something
}
else {
  do something else
}

  Unfortunately for users browsing the mail-list-archives, the the
source files I attach end up getting embedded in the message body, so
nobody can easily compile them. Hopefully one day these files will be
included with LAMMPS. Until then, people can get these files from
moltemplate.org)

I do hope this is helpful to someone.
Andrew

pair_lj_charmm_coul_charmm_inter.cpp (36.3 KB)

pair_lj_charmm_coul_charmm_inter.h (3.55 KB)

pair_lj_remix.cpp (15.9 KB)

pair_lj_remix.h (8.98 KB)

pair_lj_softcore2.cpp (9.06 KB)

pair_lj_softcore2.h (2.61 KB)

Sorry. Typo. (That should be "1 0", not "10".)

Thanks Andrew,
Could you let me know to which version of lammps these files should be added.

Hasan

I have been using this pair style with the lammps-1Jul12 version of LAMMPS.
Please let me know if it does not work with the latest version of
LAMMPS. (I don't always have time to keep my LAMMPS code up to date.)

Also let me know if you find a bug. I want to know and I will
definitely fix it quickly.

I made a web page with some updated documentation for this pair style
Go to http://www.moltemplate.org and click on "LAMMPS Code" on the left.
In the future, the most up-to-date version of this code will be
located there (unless Steve adds it to LAMMPS).
Documentation is here for now:
http://www.moltemplate.org/lammps_code/pair_lj_charmm_coul_charmm_inter.html)

Let me know if you have any difficulties.

Cheers

Andrew

Hi Andrew

I just tried this pair_style. It does not work out as the graphite atoms dissolve in the solvent. I have also tried the same simulation by defining different atom types for each sheet (molecule), and the problem did not arise and simulation went as expected. So there should be a problem with the way of using this pai_style. Below I have attached part of my input file. Could you take a look and let me know if you see any mistake in using it. The atom type for the graphite carbon atoms in which I want to exclude LJ for atoms in the smae molecules is “3”.

units real
atom_style full
dimension 3
neighbor 2 bin
neigh_modify every 1 delay 3 check yes
boundary p p p

pair_style hybrid/overlay tersoff lj/charmm/coul/charmm/inter 10 12 lj/cut 10.0
bond_style harmonic
angle_style harmonic
dihedral_style opls
read_data data.graphite

pair_coeff * * tersoff SiC.tersoff NULL NULL C
pair_coeff 1 1 lj/cut 0.0907 3.93
pair_coeff 2 2 lj/cut 0.2059 3.91
pair_coeff 3 3 lj/charmm/coul/charmm/inter 0.0 3.4 0.0 3.4 0.05 3.4
pair_coeff 1 3 lj/cut 0.06 3.665
pair_coeff 2 3 lj/cut 0.1 3.655
pair_coeff 1 2 lj/cut 0.1366 3.92

bond_coeff 1 95.535 1.53

angle_coeff 1 62.094 114

dihedral_coeff 1 1.411 -0.2710 3.145 0

special_bonds lj 0.0 0.0 0.0

Thanks,
Hasan

Dear Hasan

   I'm not happy to hear that, but thanks for the report.

Not to offend you, but I just to be sure of one thing:
   Did you give the atoms in each layer of graphine a different molecule-ID?
   (This is the number in the -second- column of the "Atoms" section
of your data file.
    Do all the graphite layers have different numbers in this column?)
   If not, then that could cause exactly the problem you are observing.

Otherwise, your commands look okay to me. They should work.
If not, then try this:
Replace these two commands avove with these instead:

Thanks!
I checked my data file. The atom IDs and Molecule IDs are correct. Therefor I will go ehead and try these to new lines. Before checking replacing these two lines, should I replace the last “1” with “-1”? I checked the formula and I think it should be “-1” for L.

Hasan

You are correct. My apologies!