Important bug with RESPA code

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS users and developers,

I discovered a fairly important bug in the RESPA code, and I have an easy way to fix it, too.
This is all for the latest release of LAMMPS (14 August, 2012).

I first ran into the bug when I was attempting to compare conservation of energy with RESPA
versus with Verlet integration. I found that the total energies were vastly different between
the two integrators. So, I did some digging in the code, and this is what I found:

The bug appears when you run with "run_style respa" and use "dihedral charmm". The bug is that
the pairwise Coulomb and van der Waals energies (and forces) are wrong as compared to running
with "run_style verlet". This happens at the very first time step, which *should*
be the same energies and forces in the RESPA and Verlet integrators because the first step
is an outer RESPA loop step.

I realized that the cause of the bug is that the 1-4 scaled Coulomb
and van der Waals interactions are evaluated *not* in the usual pair routines but in
dihedral_charmm.cpp. The scaled 1-4 energies and forces are then added onto those from
the usual pair routine.

So, the problem is that respa.cpp calls to evaluate the dihedral energy/force ***before***
it evaluates the pair energy/force. And, when the pair energy/force is called, it
sets the energy/force to zero. Thereby, the energy/force that was evaluated in
dihedral_charmm.cpp is lost.

To fix the bug, I simply moved the call to evaluate the dihedral energy/force in respa.cpp
after the call to evaluate the pair energy/force. This needs to be done in the setup, setup_minimal,
and recurse member functions of respa.cpp. I think that this shouldn't cause problems for
other dihedral potential energy function definitions, but I'm not 100% sure. I suppose some testing
is due for that. On top of that, you must make RESPA compute the dihedral energy/force
on the same RESPA loop level as that for the pair energy/force, otherwise its wrong again.

You can confirm this bug, if you like, with the rhodopsin benchmark in the $LAMMPS/bench directory
(you'll have to turn off the "npt" fix if using SHAKE).

Has anyone else experienced this bug? It's clearly something that needs to be fixed if people
want to run simulations with the CHARMM force field. The AMBER force field also uses 1-4 scaling,
right? So, yeah, for anyone trying to use RESPA in LAMMPS with these force fields, they are guaranteed
to have problems.

Best,
Adrian

2

Good catch on figuring out the problem here.

I think dihedral charmm is the only dihedral style that will
have this problem, b/c it is the only one that calculates
contributions to the pairwise energy/virial. Other styles,
like those used by the AMBER ff, do have 1/4 weigtings
but those are computed in the pair styles. CHARMM is
somewhat unique in that it assigns new sigma and epsilon
specifically to 1/4 atoms in a dihedral, so the dihedral
style is where those coeffs are stored and known.

I can think of a couple possible solutions. One is yours,
which does as you say, force you to put the dihedral
and pair calcluation at the same rRESPA level. Is that
a drawback to the way you would normally use rRESPA
for these models?

Another option might be to tally the energy/virial with
the dihedral instead of the pair, when rRESPA is being
used. It's really just a bookkeeping issue as I see it
or maybe have a flag on the dihedral style charmm command
for that. If you're not
running NPT (using the pressure), then it shouldn't even
effect dynamics, correct? You could argue conceptually
for rRESPA that those terms should be contributing at
the level they are looped over, i.e. with the dihedral.

Suggestions?

Steve

3

Good catch on figuring out the problem here.

I think dihedral charmm is the only dihedral style that will
have this problem, b/c it is the only one that calculates
contributions to the pairwise energy/virial. Other styles,
like those used by the AMBER ff, do have 1/4 weigtings
but those are computed in the pair styles. CHARMM is
somewhat unique in that it assigns new sigma and epsilon
specifically to 1/4 atoms in a dihedral, so the dihedral
style is where those coeffs are stored and known.

actually, no. the 1-4 lj parameters are read-in and stored
with the pair style. this is how it works in CHARMM: for
most pairs of non-bonded interactions, you use the regular
lennard jones, only in a few select cases, you have a
different set of (weaker) 1-4 interactions.

[...]

Another option might be to tally the energy/virial with
the dihedral instead of the pair, when rRESPA is being
used. It's really just a bookkeeping issue as I see it
or maybe have a flag on the dihedral style charmm command
for that. If you're not
running NPT (using the pressure), then it shouldn't even
effect dynamics, correct? You could argue conceptually
for rRESPA that those terms should be contributing at
the level they are looped over, i.e. with the dihedral.

Suggestions?

i think the *real* solution is to remove the 1-4 pair
contributions from the dihedral style and put it into
the pair interactions, where it belongs. i really don't
see a technical reason why it should not be done,
and can only assume that this split is for historical
reasons and long been obsoleted by improvements
in the exclusion handling in lammps.

you can determine from the exclusion bit, whether a
specific pair is a 1-4 pair and then use the 1-4 lj
parameters instead of the regular lj parameters.
you'll have to change the special_bonds alias for
charmm as well, but that should be it.

this change would also solve another nasty and
annoying issue. right now, one needs the weight
flag to make certain, that the 1-4 pair interaction
is calculated only once, since CHARMM allows
to have multiple dihedral parameter sets for a given
quartet of atoms and those are additive, thus
in lammps one would have to define multiple
dihedral interactions, but then set the weight to
0.0 for all but the first.

there may be a little breakage with backward
compatibility in a few cases, but i think in the
long run, removing this inconsistency is a good
thing. that it would automatically fix r-RESPA
is a bonus. :wink:

cheers,
     axel.

4

Comments below.

Steve

Good catch on figuring out the problem here.

I think dihedral charmm is the only dihedral style that will
have this problem, b/c it is the only one that calculates
contributions to the pairwise energy/virial. Other styles,
like those used by the AMBER ff, do have 1/4 weigtings
but those are computed in the pair styles. CHARMM is
somewhat unique in that it assigns new sigma and epsilon
specifically to 1/4 atoms in a dihedral, so the dihedral
style is where those coeffs are stored and known.

actually, no. the 1-4 lj parameters are read-in and stored
with the pair style. this is how it works in CHARMM: for
most pairs of non-bonded interactions, you use the regular
lennard jones, only in a few select cases, you have a
different set of (weaker) 1-4 interactions.

yes, you're right, the extra coeffs are read and stored by
the pair style, but dihedral_charmm is where the code
is that uses them

[...]

Another option might be to tally the energy/virial with
the dihedral instead of the pair, when rRESPA is being
used. It's really just a bookkeeping issue as I see it
or maybe have a flag on the dihedral style charmm command
for that. If you're not
running NPT (using the pressure), then it shouldn't even
effect dynamics, correct? You could argue conceptually
for rRESPA that those terms should be contributing at
the level they are looped over, i.e. with the dihedral.

Suggestions?

i think the *real* solution is to remove the 1-4 pair
contributions from the dihedral style and put it into
the pair interactions, where it belongs. i really don't
see a technical reason why it should not be done,
and can only assume that this split is for historical
reasons and long been obsoleted by improvements
in the exclusion handling in lammps.

you can determine from the exclusion bit, whether a
specific pair is a 1-4 pair and then use the 1-4 lj
parameters instead of the regular lj parameters.
you'll have to change the special_bonds alias for
charmm as well, but that should be it.

I don't see how this would work. The 1-4 special wt
factor is 0.0 for charmm, which means no 1/4 neighbors
will be in the neighbor list at all. So the pair style
doesn't see any 1/4 pairs.

Even if they were, how would the pair style determine
that those 2 atoms were in a dihedral? It is valid
for the user not to list a dihedral involving those
2 atoms as end points if they choose to. In which
case, I don't think the 1-4 LJ interaction should be
computed.

this change would also solve another nasty and
annoying issue. right now, one needs the weight
flag to make certain, that the 1-4 pair interaction
is calculated only once, since CHARMM allows
to have multiple dihedral parameter sets for a given
quartet of atoms and those are additive, thus
in lammps one would have to define multiple
dihedral interactions, but then set the weight to
0.0 for all but the first.

Are you saying the weight factor (input with the dihedral
coeffs) is something LAMMPS added? I don't
recall that. I thought it was part of the CHARMM method
of defining these extra 1-4 interactions. If so,
again, how would the pair style know what wt factor
to apply to a 1/4, since it doesn't know what dihedral
interaction(s) it is part of?

5

I don't see how this would work. The 1-4 special wt
factor is 0.0 for charmm, which means no 1/4 neighbors
will be in the neighbor list at all. So the pair style
doesn't see any 1/4 pairs.

but the 1-4 factor is *only* set to 0.0 because of the
LAMMPS version of the charmm dihedral. if you look
up the docs about CHARMM force fields, the philosophy
is that 1-4 interactions should *always* be considered,
unless in special cases and there - rather than using
a constant scale factor, you use the listed 1-4 LJ parameters.

so for charmm you should have
special_bonds lj/coul 0.0 0.0 1.0
that is how it is meant to be.

Even if they were, how would the pair style determine
that those 2 atoms were in a dihedral? It is valid
for the user not to list a dihedral involving those
2 atoms as end points if they choose to. In which
case, I don't think the 1-4 LJ interaction should be
computed.

here is a chunk for a CHARMM 27 parameter file:

NONBONDED nbxmod 5 atom cdiel shift vatom vdistance vswitch -
cutnb 14.0 ctofnb 12.0 ctonnb 10.0 eps 1.0 e14fac 1.0 wmin 1.5
                !adm jr., 5/08/91, suggested cutoff scheme

this part defines how the non-bonded interactions are computed.
i.e. you have a neighbor list skin of 2.0 angstrom, a cutoff of 12
and start switching (of LJ) starting at 10. you have dielectric
1.0 and special_bonds coul 0.0 0.0 1.0

moving on to the parameter lists.

!V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6]
!
!epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j)
!Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j
!
!atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4
!
!carbons
C 0.000000 -0.110000 2.000000 ! ALLOW PEP POL ARO
                ! NMA pure solvent, adm jr., 3/3/93
CA 0.000000 -0.070000 1.992400 ! ALLOW ARO
                ! benzene (JES)
CC 0.000000 -0.070000 2.000000 ! ALLOW PEP POL ARO
                ! adm jr. 3/3/92, acetic acid heat of solvation
CD 0.000000 -0.070000 2.000000 ! ALLOW POL
                ! adm jr. 3/19/92, acetate a.i. and dH of solvation
CE1 0.000000 -0.068000 2.090000 !
                ! for propene, yin/adm jr., 12/95
CE2 0.000000 -0.064000 2.080000 !
                ! for ethene, yin/adm jr., 12/95
CM 0.000000 -0.110000 2.100000 ! ALLOW HEM
                ! Heme (6-liganded): CO ligand carbon (KK 05/13/91)
CP1 0.000000 -0.020000 2.275000 0.000000 -0.010000
1.900000 ! ALLOW ALI
                ! alkane update, adm jr., 3/2/92
CP2 0.000000 -0.055000 2.175000 0.000000 -0.010000
1.900000 ! ALLOW ALI
                ! alkane update, adm jr., 3/2/92
CP3 0.000000 -0.055000 2.175000 0.000000 -0.010000
1.900000 ! ALLOW ALI
                ! alkane update, adm jr., 3/2/92
:

here you see, some interactions only have regular LJ,
some also have 1-4 LJ. those are the ones that should
should be scaled, but there is no constant scaling factor,
so special_bond lj 0.0 0.0 1.0 and for 1-4 interactions,
you use the 1-4 parameters.

now if you look into a "topology" file, i.e. where the
topology data for all known residues is defined, you
will see that only atom properties (atom type and charge),
bonds, impropers and CMAP terms are defined. *no*
angles or dihedrals. those are *all* inferred from the
bond topology.

thus, what you would need for charmm is to set the
1-4 lj scaling factor some some value different from
0.0 and 1.0 and then have the charmm pair styles
*ignore* that scaling factor but *use* the flag that
identifies an interaction as a 1-4 LJ interaction and
use the LJ14 parameters (which are copied from the
regular LJ, if not given on input) to compute the LJ
part of the charmm pair potentials.

[...]

Are you saying the weight factor (input with the dihedral
coeffs) is something LAMMPS added? I don't
recall that. I thought it was part of the CHARMM method
of defining these extra 1-4 interactions. If so,
again, how would the pair style know what wt factor
to apply to a 1/4, since it doesn't know what dihedral
interaction(s) it is part of?

the only explanation for the weight factor is to not double
count the non-bonded terms in the dihedral as they must
not be double counted. this term does not appear anywhere
in the CHARMM parameter or topology files and not in a
PSF file.

getting rid of this mess would make a lot of things easier
and should also make CHARMM for LAMMPS a bit faster.

axel.

6

I have to agree with Axel here, based on my experience of coding the AMBER and
CHARMM force fields in Q-Chem. CHARMM scales neither
the 1-4 Coulomb nor the 1-4 LJ interactions (i.e. weight=1), but it
*does sometimes* have certain special 1-4 LJ interactions which
should get defined by the force field parameter set or by the user somehow.
That appears to be the only reason for the existence of the
weight factor in LAMMPS.

Setting "special_bonds lj/coul 0.0 0.0 0.0" is just a device to
avoid double counting in LAMMPS, but it also incurs unnecessary overhead
checking for 1-4 interactions (or multiplying by zero) during
the pair computation when using the CHARMM force field.

AMBER is an exception, in that it does uniformly scale 1-4
interactions (LJ by 0.5 and Coulomb by 5/6), which is
automatically taken care of in force.cpp around line 556.
This is done using the special_lj[] and special_coul[]
arrays. The weight factor in AMBER is set to zero,
which is unnecessary overhead in the dihedral computation.

So, that means that LAMMPS must know that, for any given
i-j pair, whether or not it is a 1-4 interaction
during the pair computation. If that's true, there is no good
reason to do this computation in dihedral_charmm.cpp at all.
I think LAMMPS just "invented" the notion of the weight factor as
a way to do the special 1-4 LJ interactions in CHARMM. But, I don't see
why it has to be applied during the dihedral computation if
it is known during the pair computation what is and is not a 1-4 pair.

I think a solution to the problem (which I think Axel is getting at)
is to delete the pair computation entirely from dihedral_charmm.cpp
and then introduce code into the appropriate pair styles
(like pair_lj_charmm_coul_long.cpp)
that, when checking for a 1-4 pairwise interaction, applies
the weight factor if and only if using the CHARMM force field.
Perhaps we should make separate pair routines specifically for the
CHARMM force field. C++ templates might be good for that to avoid
redundant code.

Like Axel said, doing all of this will automatically correct
the issues with rRESPA and should also likely result in some speedup.

As a side note, I am looking into the rRESPA code because I want
to have rRESPA working with USER-OMP as well as the partitioning
of the verlet/split integrator. I have a preliminary version of
all of this working, including the new respa/split integrator.
However, I have had to reduce the forces
across all OpenMP threads at every RESPA level, which is
inelegant, I know, but is the easiest way to get it to go.

I'm hoping this will boost scaling and productivity.
I want to see LAMMPS be more competitive with the likes of
NAMD, where people regularly use multiple time step
integrators to boost productivity. Then, they report
how much "faster" NAMD is than LAMMPS (with its regular
old Verlet integrator) without even realizing that
it's an unfair comparison. Anyhow, I've been encountering
this kind of unfair comparison a lot lately, and I hope
to quell it to some extent. I believe in LAMMPS!

Adrian

7

I am very confused about why there is a weighting factor in
dihedral_coeff for use with dihedral_style charmm. Shouldn't
the epsilon14 and sigma14 of the CHARMM force field take care
of this? Why would there be an additional uniform weighting on top
of it?

If so, I don't see why the weighting is present in
dihedral_charmm.cpp. We can very easily do the 1-4 interactions
in the appropriate pair routine.

Adrian

8

I am very confused about why there is a weighting factor in
dihedral_coeff for use with dihedral_style charmm. Shouldn't
the epsilon14 and sigma14 of the CHARMM force field take care
of this? Why would there be an additional uniform weighting on top
of it?

yes.

If so, I don't see why the weighting is present in
dihedral_charmm.cpp. We can very easily do the 1-4 interactions
in the appropriate pair routine.

charmm allows to have multiple sets of
parameters for the same dihedral. this is
like a /multi variant, but in lammps those
are defined as explicit dihedrals, and
now you need the weight variable, to
adjust the lj and coulomb contributions
from the pair part in the dihedral calc
so they are not counted twice.

this all is extremely convoluted and
needlessly complex. if you want to
get *really* scared, have a look at
charmm2lammps.pl. i thought i knew
how to read perl code until i saw that.

as i wrote before i can only assume
that this was conceived to work around
an issue in a historical version of LAMMPS
and that somebody just got carried away
and tried to solve a problem that
doesn't exist.

axel.

9

Ah, that clears up my confusion.

I'm going to try out writing a new pair routine
for the CHARMM force field, no weight factors involved.

Adrian

10

So if I understand these various suggestions,
we should:

a) set the "charmm" default for speical bonds to be
    something like 0.0 0.0 0.5, so that 1/4s are included
    in the neighbor list
b) remove the wt flag on the dihedral charmm coeffs,
    and all 1/4 pairwise code from dihedral charmm
c) add logic to pair lj/charmm to apply the correct
    1/4 specific sigma/epsilon to each 1/4 pair it encounters
    (this works b/c that is the CHARMM-compatible
     definition of a valid interaction, independent
     of whether the 2 atoms are in 0,1,or N dihedrals?)
    does it also need the weight value that was specified
     with the matching dihedral, or is that now moot?
d) these changes get rid of the problem with rRESPA
    that Adrian noted
e) does this break the charmm2lmp converter we
    have, or mean it can be simplified to ignore
    multiplicity factors in enumerating dihedrals?
    note there is a NOTE on the dihedral style
    doc page about making sure the prefactor K
    accounts for multiplicity, which I think is there
    b/c of issues with how LAMMPS vs CHARMM
    treats K - not clear if that is related to what
    we are discussing now or not
f) all of this is unrelated to AMBER or other force
    fields that define 1/4 weightings - I don't think
    they have any issues

If this is all good, then I agree we should go ahead.
Paul - you worked on the initial CHARMM implementation
in LAMMPS including charmm2lmp with Pieter - does
this all sound good to you?

Steve

11

Yes, I'm in agreement with nearly everything that Axel and Adrian are suggesting here. It will be somewhat painful to make the needed changes to LAMMPS and the converter script, and prove to everyone's satisfaction that we still get correct energies and forces, but I agree with Axel that cleaning up the "mess" will pay off in the long run. The main thing is that before any changes to LAMMPS are released, adequate verification testing will be needed to make sure nothing gets broken.

Paul