SPC/E water and 2-Octyl-alcohol Stuck during running

LAMMPS LAMMPS General Discussion
,
1

Hello,

I had some problems running my simulation using the SPC/E water model. The program is blocked at this step and has no output：

Neighbor list info ...
  update every 1 steps, delay 10 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 16
  ghost atom cutoff = 16
  binsize = 8, bins = 11 11 11
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair lj/cut/coul/long, perpetual
      attributes: half, newton on
      pair build: half/bin/newton
      stencil: half/bin/3d/newton
      bin: standard
Setting up Verlet run ...
  Unit style    : real
  Current step  : 0
  Time step     : 1

I use GROMOS_54A7_ATB force field to simulate. The *.in file is built using moltemplate，here is the file I used to run lammps. Firstly, 30 2-Octyl-alcohol molecules were placed at a 8×8×8 nm box before solvating the box with water molecules.Then, a corresponding number of water molecules were replaced by 32 Na ions and 32 Cl ions.

Is it the wrong way I use SPC/E water model? Is there any way to overcome this problem.

Thank you for your help,
Zhu
system.data (3.3 MB)
system.in (1.5 KB)
system.in.settings (222.7 KB)

2

It would be extremely useful to see the output before it gets stuck.
Also, have you visualized the system and its evolution before this step?

3

I’ve looked at your input files. You have some close contacts in your initial geometry.
You are doing a minimization, which can help getting rid of those. However you have a rigid water potential that requires the use of fix shake, and fix shake (or fix rigid) cannot be used during minimization. That now creates a problem as your water molecule atoms can move away from each other. This can be avoided by assigning very large force constants to the water bond and angle interactions.

You should be able to see what I explained from the energies and pressure. First they are very positive, but then become extremely large and negative. That is almost always a bad sign. After that your forces will be bogus and fix shake will fail to satisfy the constraint requirements when you start the MD and the result will lead to a deadlock.

When instead you change the water force constants to, e.g.:

bond_coeff 60 100000.0 1.0 #60=Ow-Hw
angle_coeff 5 10000.0 109.47 #5=Hw-Ow-Hw

The minimization will produce much more reasonable results and the MD can start as expected.

3 Likes
4

Dear Dr. Axel,
Thank you for the great suggestion!
Zhu William