Alternatively, you could turn off ONLY the Lennard-Jones contribution
from atoms in the same molecule by hacking the
pair_lj_cut_coul_long.cpp code (or pair_lj_cut.cpp code)
IE:
if (atoms->molecule[i] == atom->molecule[j]) {
ignore these Lennard-Jones AB parameters:
lj1[itype][jtype]
lj2[itype][jtype]
lj3[itype][jtype]
lj4[itype][jtype]
and treat them as 0
#This will turn off the Lennard-Jones contribution
to the force between atom i and j
}
LAMMPS Molecular Dynamics Simulator
If it helps, I used this trick on line 199 of
pair_lj_charmm_coul_charmm_inter.cpp. You can currently download that
code here:
http://www.moltemplate.org/lammps_code/
currently LAMMPS does not allow you to have
two different rigid bodies which share the same molecule-ID number.
(At least that is my understanding.) So you'd have to hack the LAMMPS
pair_lj_cut_coul_long.cpp code so that it can tell whether two atoms
belong to the same molecule some other way. I'm thinking about
modifying the fix_rigid.cpp code to loosen this restriction. If I
make any progress on this, I'll post again.
Yes this works.
You can get around this particular problem by
1) assigning different molecule-ID numbers to the atoms in the two
different rigid bodies
2) defining a group consisting of all the atoms in either rigid body.
(IE all the atoms in the anizole molecule. In the example below, I
call this group "gRigid")
group gRigid id <> 1 16
group gRigid id <> 17 32
group gRigid id <> 33 48
etc...
3) Using fix rigid with this group (eg "gRigid") and the "molecule" bodystyle
Eg.
fix fRigid gRigid rigid molecule
4) Afterwards, use the "set" command to assign the molecule-IDs of
these atoms back to their original values (ie, one molecule-ID per
anizole molecule).
For a molecule like anizole with 16 atoms, you would use commands like this:
set atom 1*16 mol 1
set atom 17*32 mol 2
set atom 33*48 mol 3
etc...
I included versions of these commands using moltemplate.sh (to assign
the numbers) in the example below.
------- example for moltemplate users: ------
I have not tested this example, and I omitted most of the atoms from
the "Data Atoms" section, but hopefully this gives you an idea. In
this example, I define two types of objects "BenzeneOC" and "H3".
These contain the atoms for the two rigid objects in each anizole
molecule. In each "Data Atoms" section, I set the 2nd column to
"$mol:.", which assigns different molecule-ID numbers to these two
different molecule objects. I also nested the definitions of the
"BenzeneOC" and "H3" objects within the "Anizole" molecule-object.
There are several other ways to do this. You could can define these
separately (as independent molecules), if that makes more sense to
you. More importantly, you don't have to formally break the molecule
into separate objects this way (eg "BenzeneOC" and "H3"). If you
prefer, you can just make a simple molecule and use "set" commands to
set the molecule-IDs for the atoms in the benzene ring and the CH3
group. (In retrospect, perhaps this is the easiest way.)
----- file "anizole.lt" -----
Anizole {
# Define "BenzeneOC"
# This is a benzene ring with one of the hydrogen atoms replaced
# with an oxygen-carbon atom pair.
# (Don't worry about missing or dangling bonds.
# Moltemplate does not care. We will add them later.)
BenzeneOC inherits OPLSAA {
write("Data Atoms") {
$atom:C1 $mol:. @atom:90 0.00 5.274 1.999 -8.568 #AromaticC
$atom:C2 $mol:. @atom:90 0.00 6.627 2.018 -8.209 #AromaticC
etc...
$atom:O $mol:. @atom:121 0.00 4.27 5.22 -8.50 #Anisole-OCH3
$atom:CT $mol:. @atom:123 0.00 4.27 5.22 -8.50 #CT-MethylEtherCH3OR
}
# You don't need to define bonds between atoms in a rigid body
} # BenzeneOC
# Define the "H3" object
# This is a trio of hydrogen atoms
H3 inherits OPLSAA {
write("Data Atoms") {
$atom:H1 $mol:. @atom:89 0.00 4.27 5.22 -8.50 #AlkeneH-C=
$atom:H2 $mol:. @atom:89 0.00 5.22 5.21 -7.53 #AlkeneH-C=
$atom:H3 $mol:. @atom:89 0.00 5.93 5.26 -9.15 #AlkeneH-C=
}
} # H3
# Now "instantiate" (make a single copy of)
# each rigid object:
benzeneoc = new BenzeneOC
h3 = new H3
# ... and link them together with the appropriate bond(s)
write("Data Bond List") {
$bond:CT-H1 $atom:benzeneoc/oC $atom:h3/H1
$bond:CT-H2 $atom:benzeneoc/oC $atom:h3/H2
$bond:CT-H3 $atom:benzeneoc/oC $atom:h3/H3
}
# Add the atoms in the entire Anizole molecule to the "gRigid" group
write("In Settings") {
group gRigid id <> $atom:benzeneoc/C1 $atom:h3/H3
}
# Note: moltemplate assigns atom-ID numbers in the same
# order in which they were instantiated. Since we used the
# "new" command to create "benzeneoc" before "h3"
# $atom:benzeneoc/C1 is the lowest atom-ID in this molecule,
# and $atom:h3/H3 is the highest atom-ID in this molecule,
# Now write a list of "set" commands which change the molecule-ID
# numbers of all the atoms (in either rigid body) so that they share
# the same molecule-ID. (But different for each anizole molecule)
write("molIDs_merged.in") {
set atom \{atom:benzeneoc/C1\}\*{atom:h3/H3} mol $mol:anizole
}
# Note:
# "$mol:anizole" refers to a unique molecule-ID assigned to the entire anizole
# object. (I could have called it "$mol:." too, but that could be confusing.
# Recall "$mol:." refers to a unique molecule-ID for whichever object
# in which it is defined (eg the current instances of "benzeneoc" and "h3".)
} # Anizole