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.