Problem in minimization step: TIP4P model

LAMMPS LAMMPS General Discussion
1

Dear Lammps users/developer,

I am using TIP4P water model in a system including water and silicon.
The problem is that when Lammps starts to minimize the system, it stops. It stays in this step without any progress or error.

The initial part of input file is:
units metal
dimension 3
boundary p p p
atom_style full

restart 10000 restart11 restart22
read_data TIP4P.data

set type 1 charge 0.520
set type 2 charge -1.040
set type 3 charge 0.000

bond_style harmonic
angle_style harmonic
dihedral_style none
pair_style hybrid lj/cut/tip4p/long 2 1 1 1 0.15 10.0 12.0 lj/cut/coul/long 10.0 12.0 meam/c

bond_coeff 1 100.0 0.9572
angle_coeff 1 10.0 104.52

pair_coeff * * meam/c library.meam Si Si.meam NULL NULL Si
pair_coeff 1 1 lj/cut/tip4p/long 0.00 0.00 #H-H
pair_coeff 1 2 lj/cut/tip4p/long 0.00 0.00 #H-O
pair_coeff 2 2 lj/cut/tip4p/long 0.006721 3.1536 #O-O
pair_coeff 2 3 lj/cut/coul/long 0.120878 2.6305 #Si-O
pair_coeff 1 3 lj/cut/coul/long 0.00 1.0475 #H-Si

kspace_style pppm/tip4p 1.0e-5

group SiSi type 3
group water type 1 2

minimize 1.00E-15 1.00E-15 1000000 10000000
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.001

fix 1 water shake 1.00E-06 100 0 b 1 a 1
velocity all create 498.15 4928459 rot yes dist gaussian

The initial part of data file also is:
37658 atoms
19644 bonds
9822 angles
0 dihedrals
0 impropers
3 atom types
1 bond types
1 angle types
0 dihedral types
0 improper types
-86.53004 86.53004 xlo xhi
-10.5651335 10.5651335 ylo yhi
-58.907505 58.907505 zlo zhi

Pair Coeffs

1 H
2 O
3 Si

Bond Coeffs

1 H-O

Angle Coeffs

1 H-O-H

Masses

1 1.008000 # H
2 15.999400 # O
3 28.084999 # Si

Atoms # full

1 1 2 -1.040000 14.856773 6.031134 41.662369 # O O
2 1 1 0.520000 14.075775 5.448134 41.771366 # H H
3 1 1 0.520000 14.515778 6.807135 41.171368 # H H
4 2 2 -1.040000 59.047783 -1.556866 50.970367 # O O
5 2 1 0.520000 59.602776 -2.161866 50.435364 # H H
6 2 1 0.520000 59.327782 -1.719866 51.894363 # H H
7 3 2 -1.040000 -37.135223 -6.559866 31.849365 # O O
…

and the last part of log file is:

Setting up cg style minimization …
Unit style : metal
Current step : 0
Per MPI rank memory allocation (min/avg/max) = 19.48 | 24.78 | 33.17 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 1.1557718e+16 1017.6889 1.1557718e+16 1.7192536e+17

I appreciate any comment to help.

Thank you.

If it is required, the initial pdb file that I used to add water to the system by packmol was:
HETATM 1 O 1 -2.688 2.013 0.000
HETATM 2 H 1 -3.582 1.611 0.000
HETATM 3 H 1 -2.076 1.248 0.000
CONECT 3 1
CONECT 2 1
CONECT 1 3 2
END

2

Please note that you have an extremely high potential energy and pressure.
This is likely due to close contacts or other issues with your input.
It is impossible to say something more specific with so limited information.
For starters, you should always report what version of LAMMPS you are using and platform you are running on.

This is not related to your system stalling, but it is not necessary to use pair style lj/cut/coul/long, in fact it would be incorrect, if you have coulomb interactions between the TIP4P water atoms and Si.

3

Thanks Axel.

The version of Lammps is here:
load modules
module load lammps/lammps-20200505-intel-19.0.4-pk3fyur
I am using mpirun with 64 cpus.

When I use SPC/E model, in the minimization step I have below information and according to your explanation, it also is high but this works without any problem.
Per MPI rank memory allocation (min/avg/max) = 18.23 | 23.23 | 30.66 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 5.1164403e+15 1544.5927 5.1164403e+15 7.6108954e+16
2087 0 -43156.204 80.018742 -43076.186 2233.7386
Loop time of 385.133 on 64 procs for 2087 steps with 37658 atoms

The system with TIP4P model is exactly the same as that of SPC/E model, number of atoms, bonds and so on. Just I changed the numbers in the initial pdb file of water to use in packmol which I have added in the main part of my question, then I corrected the commands in input file and also in the data file O is the first atom and clearly the numbers and orders also have been changed accordingly. No problem I see when I check the file. So I really am confused why it does not work.

4

This is not what I was asking for. Please see: 1.2. What does a LAMMPS version mean — LAMMPS documentation

You just got lucky on that one. You still have a bad initial configuration with possible close contacts, too.
Even after the minimization you have a very high pressure, so I suspect that your initial simulation box is too small.

I cannot read minds, neither of people nor of computers, and there is not enough information here to make a proper assessment. I already mentioned by best guess of a too small box which can cause close contacts due to using periodic boundary conditions.

Due to the different balance between (repulsive) LJ and (attractive) Coulomb for inter-molecule interactions, it is quite possible that your forces for TIP4P may overflow and thus the minimization will not converge. You can check on this by outputting thermo data in every step.

5

Thank you Axel. I take both your comments.

6

Axel,

I took both your comments.
The system is working with increasing a little bit the size of box, no stop, anymore.
I also have removed the extra pair_style as you mentioned between Si and water atoms.
Now the pair_style and pair_coeff are as below. I am using the TIP4P model with a long range Coulombic solver, the last model of 8.4.3. TIP4P water model — LAMMPS documentation.
Now I just wanted to be sure if the cut off of LJ and Coulomb interactions make sense or any part of below code needs correction?
…
set type 1 charge 0.5242
set type 2 charge -1.0484
set type 3 charge 0.000

bond_style harmonic
angle_style harmonic
dihedral_style none
pair_style hybrid lj/cut/tip4p/long 2 1 1 1 0.1250 10.0 9.5 meam/c

bond_coeff 1 100.0 0.9572
angle_coeff 1 10.0 104.52
pair_coeff * * meam/c library.meam Si Si.meam NULL NULL Si
pair_coeff 1 1 lj/cut/tip4p/long 0.00 0.00 #H-H
pair_coeff 1 2 lj/cut/tip4p/long 0.00 0.00 #H-O
pair_coeff 2 2 lj/cut/tip4p/long 0.00706 3.164 #O-O
pair_coeff 2 3 lj/cut/tip4p/long 0.1237 2.630 #Si-O
pair_coeff 1 3 lj/cut/tip4p/long 0.00 1.0475 #H-Si

kspace_style pppm/tip4p 1.0e-5

Here, 8.4.3. TIP4P water model — LAMMPS documentation, it has been mentioned that with the kspace_style pppm/tip4p, which is for a long-range model, the OM distance should be 0.1250 while in the Examples/tip4p, the model seems long range but the OM distance is 0.1546. Also in the bond_coeff and angle_coeff, both Ks are 0. I am asking about the example just to be sure about my model.

I appreciate if there will be any comment.

Thank you.

7

The OM distance for the original rigid TIP4P model is 0.15 \AA.
The page then has a list of variants. 0.125 \AA is for TIP4P-Ew.
Please note that these are all different parameterizations where not only the OM distance differes, but also partial charges and Lennard-Jones epsilon and sigma.
In each case the researchers have looked for the optimal balance to reach the best fit at the desired conditions and simulation boundary conditions.

The TIP4P pair and kspace styles can handle any of these models. To determine which is most suitable, you need to study the corresponding publications and ideally also search for publications that compare the different variants under the conditions that you are interested in.