ERROR on proc 0: Shake determinant = 0.0 (../fix_shake.cpp:2315)

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I am trying to diffuse water molecule SPC/E model inside MOF. I modeled my MOF as rigid (no bonds, no angles, only coordinates and pair coeff ) and excluded pairwise interaction between MOF atoms by “neigh_modify exclude MOF MOF”.

First i created 10 water molecules randomly inside MOF. Then, to separate overlapped atoms i run fix nve/limit for 50000 steps with timestep of 1 fmsc (I think it is enough step to separate overlapped atoms because etotal is almost constant after around 10000 step). I used nve/limit because if use minimization command, it removes atoms far apart from each other and ends with “lost atoms error”. Here i observed that there is ±1500 atm fluctuation in pressure (Maybe this is the problem?). Finally i run molecules in nvt ensemble by holding them rigid with fix shake. But i got an error message : ERROR on proc 0: Shake determinant = 0.0 (…/fix_shake.cpp:2315).

If i remove fix shake command and run in normal nvt ensemble, simulation works fine. I run same simulation (with rigid water and rigid MOF) in RASPA without problems. I looked several mail archives and manual but i couldn’t solve the problem. Can u give some ideas?
thanks

Here my input:
clear
#Initialization
dimension 3

periodic boundries

boundary p p p
units real
atom_style full

pair_style lj/cut/coul/long 14 14
pair_modify mix geometric tail yes
kspace_style ewald 1.0E-6
bond_style harmonic
angle_style harmonic
special_bonds lj/coul 0.0 0.0 0.5

#Read atom cinfigurations
read_data framework.data extra/atom/types 2 &
extra/bond/types 1 &
extra/angle/types 1 &
extra/bond/per/atom 2 &
extra/angle/per/atom 1 &
extra/special/per/atom 2

molecule h2omol H2O.txt toff 8
create_atoms 0 random 10 25143 NULL mol h2omol 464563

mass 9 15.9994
mass 10 1.0

pair_coeff 1 1 0.01499 2.69007
pair_coeff 2 2 0.169958 3.12
pair_coeff 3 3 0.2099498 2.96
pair_coeff 4 4 0.1478 3.61815
pair_coeff 5 5 0.1478 3.61815
pair_coeff 6 6 0.1478 3.61815
pair_coeff 7 7 0.0382887 2.44833
pair_coeff 8 8 0.169959 3.12
pair_coeff 9 9 0.15535 3.166
pair_coeff 10 10 0.0000 0.0000

bond_coeff 1 1000.0 1.0
angle_coeff 1 100.0 109.47

group sorbate type 9:10
group mof type 1:8

neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes exclude group mof mof

velocity sorbate create 300.0 54654
fix overlapRemoved sorbate nve/limit 0.1

thermo_style custom step temp press fmax ke etotal
thermo 100
#dump w2 sorbate image 1000 image.*.jpg type type zoom 1.6 adiam 1.2

timestep 1
run 50000

unfix overlapRemoved
reset_timestep 0

fix 2uu sorbate nvt temp 300.0 300.0 $(dt*100)
fix wshake sorbate shake 0.0001 20 0 b 1 a 1 mol h2omol

thermo_style custom step temp press fmax ke etotal
thermo 100

#dump w2A sorbate image 1000 imageAfter.*.jpg type type zoom 1.6 adiam 1.2

run 50000

log.lammps (42.9 KB)

2

Your procedure of unoverlapping the added molecules can be best described as “brutal”.
I favor using a gentler approach starting with a soft (or soft core) potential which can be gradually raised using fix adapt and then switching to the intended potential.

the shake error can be due to a largely deformed water molecule, so the solver cannot solve the constraint equations.
the chance for this happening can be largely reduced by increasing the force constants for bond and angle during minimization
and also sometimes by doing a “run 0” to before starting a regular run. the latter “trick” is to work around the issue that the shake solver in LAMMPS assumes that you are starting from a reasonable initial configuration. the “run 0” can help to enforce the constraints without also doing time stepping and thus improve the initial guess for subsequent MD time stepping. if your initial structure is still of very high potential energy, that may not be sufficient and a gentler method to do unoverlapping may be required.

when comparing LAMMPS to other MD codes, you also need to factor in efficiency in parallelization. LAMMPS favors parallelization over stability and when starting from extreme conditions it may be more likely to terminate due to instablity than other codes, who have more options to recover from bad geometries (at the expense of worse parallel scaling). thus if your approach works better with some other MD code, why not use that MD code to generate an initial, resonably well equilibrated structure and then transfer the configuration to LAMMPS afterwards?

axel.