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