SHAKE TIP3P/CHARMM issues

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I’m having an issue currently where I cannot apply fix SHAKE to water molecules in an nve/limit simulation to minimize systems. If I’m using pair_style lj/charmm/coul/long I won’t get any steps as it crashes instantly (Out of range atoms - cannot compute PPPM), however if I use lj/charmm/coul/charmm I can get a few steps in before it crashes. My SHAKE stats are the following:

fix 1 TIP3 shake 0.001 20 1 b 56 a 129

fix 2 system nve/limit 0.1

(system group contains TIP3 molecules as well)

SHAKE stats (type/ave/delta) on step 0
56 0.957222 0.00199436
129 104.516 0.19229

SHAKE stats (type/ave/delta) on step 1
56 -nan 1.39977e-05
129 -nan 0.0011278

I can get the system to run fine if I don’t use SHAKE, however as I’m using CHARMM I would rather be able to use SHAKE.

2

What is the problem? Just run without shake for a bit until your system is somewhat relaxed and then turn shake on and relax some more.

Perhaps it would be helpful to check for overlaps and tweak the box size.

Axel

3

Dear Kamron

It is always a good idea to relax and equilibrate your system before
starting a simulation. If disabling "fix shake" helps you do this,
that's great.

People often use the "minimize" command as a first step towards
equilbrating the system. Here's an example of a mixture of tip3p
water:

# --- atom definitions ---
read_data your_data_file

# --- settings: ---
group tip3p type 1 2
fix fShake tip3p shake 0.0001 10 100 b 1 a 1
# You may have to change the atom, bond, and angle type numbers
# depending on your system. (Here I assume they are all "1" and "2").

#... do some other stuff ...

# --- run section: ---

# Minimization protocol:

# If you already called "fix shake", you must "unfix" it before invoking
# the "minimize" command. (minimize is incompatible with fix shake)
# http://lammps.sandia.gov/doc/fix_shake.html

unfix fShake

# Now minimize the system using a command like this:
minimize 1.0e-5 1.0e-7 100000 400000

# Now turn "fix shake" back on again:
fix fShake tip3p shake 0.0001 10 100 b 1 a 1

# ...and run the simulation

fix fxnpt all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 1.0

# if you want to be really careful, you can start with small
timesteps, and increase them
timestep 0.5
run 2000
timestep 1.0
run 2000
timestep 2.0

# ...and run the simulation

thermo 100
dump 1 all custom 1000 traj_npt.lammpstrj id mol type x y z ix iy iz

run 50000

# optional:
write_data system_after_npt.data
write_restart system_after_npt.rst

Cheers

Andrew

4

I don’t think you can use fix nve/limit with fix shake (which we should document
as a restriction on the fix nve/limit doc page).

The reason is that fix shake is adjusting forces based on the assumption
that the next time step will move atoms to a particular position. It is correcting
the forces so the new positions will end up satisfying the SHAKE constraints.

But nve/limit is chaning the timestep so atoms won’t move to where fix shake
thought they would.

Axel’s suggestion is much better - don’t apply SHAKE until your system
is unoverlapped.

Steve