Modeling a Partially Rigid Molecule

Hello,

I’m doing a study on substituted benzenes using the OPLS potentials as given in Jorgensen 1992. One of these molecules is anisole, which is a benzene ring with one carbon connected to an oxygen which is connected to a united CH3 group. The molecule is supposed to be rigid, with the one exception that the CH3 group can rotate around its dihedral angle.

Since this is a kind of hybrid molecule in which there are both rigid and nonrigid components, an obvious manner in which this can be modeled hasn’t yet occurred to me. I’ve been tinkering with the implementation of this in LAMMPS. A few thoughts have occurred to me:

  1. As per the “fix rigid” documentation, I tried using a “fix rigid/NPT” on all atoms except the CH3 group, and then used a “fix NVT” on the CH3 group. However, the problem with this is that there appears to be no manner in which LAMMPS will treat the bond length/angle between the CH3 group and the rest of the molecule. Therefore, the bond length and bond angle begin changing. I see no way to remedy this treatment, and I don’t really think this method of using 2 separate “fixes” was meant to be used on a single molecule anyway.

  2. As per the “peptide” example provided in the LAMMPS distribution, one could use a “fix NPT” on the entire molecule, but constrain certain bonds and angles via “fix shake.” However, the problem with this method is that due to the limitations of fix shake, one can’t fully constrain the entire benzene ring; when I ran the peptide example, I observed the C-C bond distance in the aromatic group fluctuating between 1.34A and 1.44A since it could not be set (the author of the example chose to “fix shake” the C-H bond distance in the aromatic group instead).

  3. As per multiple suggestions on the mailing list, one could simply use “fix NPT” on the whole molecule, and set pseudo bonds, angles, and dihedrals throughout the rest of the molecule with very large harmonic constants. This should work if I use a short enough timestep, but it’s very ugly and inelegant.

I’m wondering whether anyone knows of a solution to this “partially rigid molecule” problem, as I’m sure it’s a class of problem widely encountered, yet I couldn’t find anything about it on the mailing list.

As always, my gratitude to this fantastic mailing list.

Efrem Braun

Graduate Research Assistant
University of California, Berkeley

I am not aware of any restriction in fix rigid which prevents atoms
from two different rigid bodies from interacting with each other via
bonded interactions (bonds, angle, dihedral, or improper
interactions).
As far as LAMMPS is concerned, I would guess that it just treats these
as other forces acting on these particles (just as it would treat
pairwise repulsive forces when they bump into each other). It seems
you are aware of this.

In your case, the first rigid body could be the
benzene-ring+oxygen+carbon, and the second rigid body could be the 3
remaining hydrogens. (Or you could make the second body the entire
CH3 group.)

However if, in addition to this, you want the length of the bond(s)
connecting them to be fixed (and/ the bond-angle to be frozen), then
that could be a problem. I don't know if you can combine fix shake
with fix rigid to constrain the distance of a bond connected two rigid
bodies (for example). It might not work, but have you tried it?

But the easiest thing to do is us fix rigid, and use a fairly stiff
bond (and bond-angle) connecting the two rigid bodies. You thought
this was inelegant, but how stiff do you really need it be? You can
probably get away with a bond which is 4 times stiffer (K is 4 times
larger) than you normally would use. This does not cost you anything.
This is because fix rigid forces you to use smaller timesteps (0.5fs
instead of 1fs). If you need something stiffer, lower the timestep.

This isn't the answer you wanted, but I doubt LAMMPS can do it another
way (unless you were somehow able to combine fix rigid with shake).

If I'm wrong, hopefully somebody will correct me
Good luck!

andrew

I am not aware of any restriction in fix rigid which prevents atoms
from two different rigid bodies from interacting with each other via
bonded interactions (bonds, angle, dihedral, or improper
interactions).

This is correct. You can have 2 (or more) rigid bodies, connected
by bonds, angles, dihedrals, etc. From the perspective of
the rigid body, and bonds, angles, etc that include atoms outside
the rigid body, just induce external forces on the body, same
as pairwise forces due to an external atom.

Steve

Thanks for the replies. It seems that the recommendation is that I use method 1 (doing 2 separate fixes on the 2 different parts of the molecules). I tried increasing the bond and angle harmonic constants for the bond and angle connecting the 2 parts of the molecule, and it indeed worked to keep them approximately rigid, so the method appears valid. As only 1 bond and 1 angle are using large harmonic constant, I don’t think this is too inelegant a method. (As a side note, using fix shake on the bond and angle did not work, as the screen output was “WARNING: Shake determinant < 0.0”, and eventually the program crashed with a segmentation fault, but this is not something I’m concerned about.)

However, this method of a treating a partially flexible molecule introduces a separate problem. The dihedral angle in my molecule (the only internal degree of freedom the molecule has that is meant to be flexible) does not fluctuate around its equilibrium value of 0 degrees, but rather around a value of 90 degrees. The reason for this is that the LJ and Coulombic forces acting on the united atom CH3 group (the terminal atom in the dihedral) from WITHIN the molecule dominate the dihedral_style energy. Of course, excluding pairwise interactions within the molecule would fix this, but as I’m using a long-range kspace solver, I’m unable to do so.

This issue appears unique to the case of a partially flexible molecule with a long-range solver; in a fully rigid molecule, one need not worry about having to turn off pairwise interactions, while in a fully flexible molecule, most force fields have one taking into account LJ and Coulomb forces within the molecule as long as the atoms are separated by more than 3 bonds.

Any ideas? The only thing I can think to fix this is to go to method 3 (making the molecule fully flexible but with large harmonic constants for all bonds, angles, and dihedrals except the “flexible” dihedral). However, as stated earlier, I’d prefer to avoid this method.

Does the force field you are using allow you the freedom to turn these
pair interactions off using the "special_bonds" command? For example:

special_bonds lj/coul 0 0 0
This will turn off all pairwise interactions between pairs of atoms
connected by 1, 2, and 3 bonds (also known as "1-2", "1-3", and "1-4"
interactions, respectively)
http://lammps.sandia.gov/doc/special_bonds.html

If this excludes too many pairwise interactions, you can always them
back again (one pair at a time) by adding explicit bonds connecting
effected pairs of atoms using bond_style hybrid table. You can adjust
the table so that the force matches the lennard-jones interaction
which would have been present had you not excluded it. (Although it's
a pain to do this for every pair, I admit.)

It's too bad "neigh_modify ... exclude" does not work in this case.
Just turning off lennard-jones repulsion might be enough to solve your
problem. (I'm running into that issue too with a different
simulation.) I wonder if using neigh_modify exclude in this case
would be equivalent to turning off pairwise short-range lennard-jones
(and coulombic?) interactions, whilst leaving the long-range coulombic
interactions in place. If "neigh_modify exclude" effectively does
this, then perhaps it's okay. Is this a stupid idea? I'm curious to
hear back from someone who knows more about how long-range solvers
work than I do.

Cheers!

Andrew

Does the force field you are using allow you the freedom to turn these
pair interactions off using the "special_bonds" command? For example:

special_bonds lj/coul 0 0 0
This will turn off all pairwise interactions between pairs of atoms
connected by 1, 2, and 3 bonds (also known as "1-2", "1-3", and "1-4"
interactions, respectively)
http://lammps.sandia.gov/doc/special_bonds.html

If this excludes too many pairwise interactions, you can always them

​please note that special_bonds lj/coul 0 0 0 is the default and thus ​all
pair interactions between particles connected via bonds, angles or
dihedrals are excluded from the neighbor list.

back again (one pair at a time) by adding explicit bonds connecting
effected pairs of atoms using bond_style hybrid table. You can adjust
the table so that the force matches the lennard-jones interaction
which would have been present had you not excluded it. (Although it's
a pain to do this for every pair, I admit.)

It's too bad "neigh_modify ... exclude" does not work in this case.
Just turning off lennard-jones repulsion might be enough to solve your
problem. (I'm running into that issue too with a different
simulation.) I wonder if using neigh_modify exclude in this case
would be equivalent to turning off pairwise short-range lennard-jones
(and coulombic?) interactions, whilst leaving the long-range coulombic

​neigh_modify exclude will remove pairs from the neighbor list​ and thus
all short range non-bonded interactions are excluded.

interactions in place. If "neigh_modify exclude" effectively does
this, then perhaps it's okay. Is this a stupid idea? I'm curious to
hear back from someone who knows more about how long-range solvers
work than I do.

​long-range solvers pay no attention to exclusions. thus even if you
exclude a pair i,j all interactions with all periodic images of i,j for
both i and j will be included. in the long-range solver you compute the
interaction of each particles with an infinite lattice of *all* particles
in the system (including itself).

in the ewald style this is done directly via a lattice sum. in pppm this is
approximated by collecting a charge density distribution on a 3d grid and
then computing the potential in reciprocal space and ​fourier transforming
it back to real space. so by its very nature, you compute a manybody
interaction and not a pairwise additive potential and thus you cannot apply
the logic used in real space.

axel.​

1 Like

My forcefield actually “wants” me to turn off ALL intramolecular pair interactions. This is where the problem lies, as its the atoms beyond 3 bonds away that still have a strong enough effect such that the dihedral angle ends up being 90 degrees rather than 0 degrees.

2 thoughts:

  1. I actually can’t use the method 3 from above (making a fully flexible molecule but with high harmonic constants on all bonds, angles, and dihedrals except the one dihedral I’m interested in), because I would still have the same problem of atoms from beyond 3 bonds away having an overpowering effect on the dihedral.

Thus, to my eye, it appears that there is no way for LAMMPS to handle this type of molecule (I’d be very happy if somebody could prove me wrong).

  1. However, I’ve come up with the following “hack,” which I’d like your thoughts on. Before this hack, I had 1 bond (CH3-O), 1 angle (CH3-O-C), and 1 dihedral (CH3-O-C-C), with the first 2 having high harmonic constants and the dihedral having the form given by OPLS. The “hack” involves adding 11 additional bonds, in which CH3 gets a “bond” to EVERY atom in the molecule that it wasn’t already bound to (the 6 Cs and the 5 Hs), each of which have a harmonic constant of 0.0. This way, “special_bonds lj/coul 0 0 0” turns off intermolecular interactions between CH3 and all other atoms in the molecule, making the only interactions that get counted the original 3 (the bond and angle with high harmonic constants and the dihedral).

I’ve tried this method out, and it sort of appears to do the trick; the dihedral angle oscillates around 0 degrees as it should and the rest of the molecule remains rigid, BUT my NPT simulation comes to a density rather far off from Jorgensen’s (1037 kg/m^3 instead of Jorgensen’s simulated 992). I’m curious as to whether this hack won’t severely mess up the long-range kspace solver (or something else) in some unanticipated manner; I don’t think it will, since I’m using “special_bonds” instead of “neigh_modify exclude,” but since this is a really strange hack, I can’t be sure that something is going on that I don’t see.

If anyone has a more proper way that LAMMPS should handle this type of molecule or has a comment on this hack, I’d much appreciate it.

CH3
/
O

H C H
\ / \ /
C C

C C
/ \ /
H C H

H

My forcefield actually "wants" me to turn off ALL intramolecular pair
interactions.

​that is easy to do. make sure that all molecules have unique molecule ids,
make a group for these molecules and use:

neigh_modify ​exclude molecule group-id

axel.

Hi Axel,

In an earlier message I had written that I’m using a long-range kspace solver, so I don’t believe that using “neigh_modify exclude molecule” is possible.

Hi Axel,

In an earlier message I had written that I’m using a long-range kspace solver, so I don’t believe that using “neigh_modify exclude molecule” is possible.

Of course it is possible.

Stan - does this cause some problem I’m not seeing? I.e. if
pairwise interactions for all pairs within a single molecule
are excluded, does that mess up something for the pair + KSpace
solver in terms of what the total eng of the system is?

Steve

Hello,

A brief update on this problem. I used the hack I mentioned to model both
anisole and benzene (benzene just as a test case) using the OPLS force
field given in Jorgensen 1992 ("Monte Carlo Simulations of Pure Liquid
Substituted Benzenes with OPLS Potential Functions"). The densities,
intermolecular potential energy, and intramolecular potential energy
(comparing my dihedral energy to Jorgensen's intramolecular potential to
neglect the effects of the stiff bond and stiff angle) all match what
Jorgensen reported in Tables IV and VI of his paper. This gives me some
confidence that the hack works. I've attached the benzene.lt moltemplate
file, the system.lt moltemplate file, and the system.in file that I used.

Steve, can I interpret your last message to mean that you see no problem
with the hack, but that you're just not 100% sure? I'd like to be a bit
more confident that this method isn't messing up the KSpace solver, though
it doesn't appear to be...

Efrem Braun

anisole.lt (5.5 KB)

system.in.old (1.18 KB)

system.lt (690 Bytes)

Your approach seems to work. I guess the only trouble is that it
might be difficult to explain to potential paper referees exactly what
effect "neigh_modify exclude" has on the forces between atoms in
different rigid objects (but the same molecule).

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)

http://lammps.sandia.gov/threads/msg27855.html

This is easier to explain in your paper. (I don't know if it would
work. You'd have to try it.)

The only problem is that 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.

Cheers

Andrew

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
      }

http://lammps.sandia.gov/threads/msg27855.html

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

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.

​please don't. ​

as far as LAMMPS is concerned, a ​molecule id​ is just a number that is
attached to a specific atom. there is no need for it to coincide with
physical molecules. so just use the molecule id as needed to define your
rigid bodies. if you still need an additional identifier for individual
molecules, this can be easily added via fix property/atom.

​modifying fix rigid is a tricky business with many subtle and sometimes
unexpected issues.​

Yes this works.

You can get around this particular problem by

​this is all far too complex and not really necessary.​

​[...]​

Hi Axel,

Do you have a best-practices suggestion for how to implement this kind of partially flexible molecule in LAMMPS? I’ve been using my “hacked” version of anisole since my last post (in intramolecular energies are eliminated by placing bonds between ALL atoms in the molecule with an epsilon value of 0), and I’ve been finding that my intramolecular and intermolecular energies from LAMMPS match the values from Jorgensen’s paper, so I believes that this method works (though it’s clearly ugly). If you have a better method, I’d love to use it.