[lammps-users] respa with lj/charmm/coul/long and inner, middle, outer

It looks like there may be a bug in respa. I am trying to use it with the pair_style lj/charmm/coul/long. Using the parameters suggested in the manual for biomolecular simulations with just a box of TIP3P water with or without SHAKE (timestep 4.0 and respa 2 2 inner 1 4.0 5.0 outer 2 or respa 4 2 2 2 inner 2 4.5 6.0 middle 3 8.0 10.0 outer 4), like charges are collapsing onto each other and the energy is huge. This suggests that there is a Coulombic force with the wrong sign. Has anyone else seen this? This happens regardless of the size of the outer time step or whether pppm or ewald is used. The problem is likely in the Coulombic pair forces. Without the inner, middle, outer keywords or with the Verlet integrator there seems to be no problem in NVT or NVE.

Any ideas? I just started using LAMMPS, so maybe I am just doing something stupid.

Thanks,
Brian

Here is my script (about as simple as possible):

Created by charmm2lammps v1.8.1 on Sat Feb 28 19:33:32 CST 2009

units real
neigh_modify delay 2 every 1

atom_style full
bond_style harmonic
angle_style charmm
pair_style lj/charmm/coul/long 8 10
pair_modify mix arithmetic
kspace_style pppm 1e-4

read_data water.data

special_bonds charmm
fix 1 all nvt 300.0 300.0 100.0
#fix 1 all nve
fix 2 all shake 1e-6 500 0 m 1.0 a 1
velocity all create 300.0 29674992 dist gaussian

thermo 10
thermo_style multi
timestep 4.0
#FOR SHAKE
#run_style respa 2 2 bond 1 dihedral 1 angle 1 pair 2
run_style respa 2 2 inner 1 4.0 5.0 outer 2
#NO SHAKE
#run_style respa 4 2 2 2 inner 2 4.5 6.0 middle 3 8.0 10.0 outer 4

restart 200 water.restart1 water.restart2
dump 2 all dcd 200 water.dcd
dump_modify 2 unwrap yes

run 16000

I'll take a look at this

Steve

I'll post a patch for this within a couple days. rRESPA was broken with
an earlier change to neighbor lists. First need to post another
feature addition.

Steve

Just posted the patch for this (7Mar09). However, I noticed
another anomaly. When using dihedral charmm with rRESPA,
the pairwise energy/virial contribution it computes (special 1-4
interactions) is
not tallied correctly with the output energy/pressure.
This is just a small-size bookkeeping issue (doesn't affect dynamics, unless you
are running NPT which relies on the pressure), and isn't an issue with normal
Verlet timestepping (non-rRESPA), but I need to think how
to fix it.

Steve