No response for NVT/NPT simulations after minimisation

Dear LAMMPS users,

I was trying to run a small system of ionic liquids and TIP4P water (built by PACKMOL) with OPLS-AA force field, but there was no “dynamic response” either in NVT ensemble or in NPT ensemble after energy minimization. However it was indeed still running on HPC and did not crash with any error. The job was later cancelled by myself since there was no “response” and still cost computation resource. I also checked there is also no atom overlap in the system. The latest version LAMMPS brought me the same results. The software works well since I could run it well with other jobs or using ReaxFF force field. Could you please advise me about why NVT/NPT simulation did not run in this case? Was this problem caused by bad initial configuration?

Attached please kindly find my input file and data file. Any suggestions would be greatly appreciated!

LAMMPS (22 Aug 2018)

OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (…/comm.cpp:87)

using 1 OpenMP thread(s) per MPI task

Reading data file …

orthogonal box = (-1.0057 -1.00013 -1.00074) to (31.0057 30.9999 30.9993)

5 by 2 by 2 MPI processor grid

reading atoms …

3592 atoms

scanning bonds …

3 = max bonds/atom

scanning angles …

6 = max angles/atom

scanning dihedrals …

10 = max dihedrals/atom

scanning impropers …

2 = max impropers/atom

reading bonds …

2797 bonds

reading angles …

3175 angles

reading dihedrals …

3400 dihedrals

reading impropers …

85 impropers

Finding 1-2 1-3 1-4 neighbors …

special bond factors lj: 0 0 0

special bond factors coul: 0 0 0

4 = max # of 1-2 neighbors

7 = max # of 1-3 neighbors

17 = max # of 1-4 neighbors

18 = max # of special neighbors

2334 atoms in group water

1258 atoms in group ils

Finding 1-2 1-3 1-4 neighbors …

special bond factors lj: 0 0 0.5

special bond factors coul: 0 0 0.5

4 = max # of 1-2 neighbors

7 = max # of 1-3 neighbors

17 = max # of 1-4 neighbors

18 = max # of special neighbors

PPPM initialization …

extracting TIP4P info from pair style

using 12-bit tables for long-range coulomb (…/kspace.cpp:321)

G vector (1/distance) = 0.225743

grid = 15 15 15

stencil order = 5

estimated absolute RMS force accuracy = 0.0213383

estimated relative force accuracy = 6.42596e-05

using single precision FFTs

3d grid and FFT values/proc = 1960 192

Neighbor list info …

update every 1 steps, delay 0 steps, check yes

max neighbors/atom: 2000, page size: 100000

master list distance cutoff = 14.25

ghost atom cutoff = 14.25

binsize = 7.125, bins = 5 5 5

1 neighbor lists, perpetual/occasional/extra = 1 0 0

(1) pair lj/cut/tip4p/long, perpetual

attributes: half, newton on

pair build: half/bin/newton

stencil: half/bin/3d/newton

bin: standard

Setting up sd style minimization …

Unit style : real

Current step : 0

Per MPI rank memory allocation (min/avg/max) = 14.72 | 15.19 | 15.84 Mbytes

Step Temp E_pair E_mol TotEng Press

0 0 26639.624 646.82851 27286.452 247260.62

100 0 4540.3368 -107111.96 -102571.62 383210.94

Loop time of 1.23187 on 20 procs for 100 steps with 3592 atoms

97.6% CPU use with 20 MPI tasks x 1 OpenMP threads

Minimization stats:

Stopping criterion = max iterations

Energy initial, next-to-last, final =

27286.4523604 -101370.450586 -102571.618972

Force two-norm initial, final = 14842.4 26136.1

Force max component initial, final = 3644.63 7794.95

Final line search alpha, max atom move = 5.08344e-06 0.0396252

Iterations, force evaluations = 100 199

MPI task timing breakdown:

Section | min time | avg time | max time |%varavg| %total

in.npt (5.21 KB)

data.new17 (395 KB)

Please note, that your system’s pressure is extremely high with ~250Matm
that means that likely your initial geometry is too dense or your box too small or both.
This will lead to extremely high forces, which possibly keep the shake algorithm from converging resulting in a stalled simulation.

The best approach would be to rebuild the system (and check/adjust the box accordingly).

The second best approach would be to allow the box to relax during minimization and minimize for more steps.
Also you need to have some more tight thermostat coupling in your initial MD run, so that the temperature does not increase too much. it is best to first run with a fixed volume before switching to npt or else, you may end up with too low a density and will have to wait for too long a time until the volume/pressure is converged again (since it is much easier for a system to expand than to compress over a larger amount of volume).

with a better relaxed initial structure, it is also possible to run with a much larger timestep (1fs since you are using fix shake)
i am attaching an accordingly modified input, however, when visualizing the trajectory, it looks like your system would be better off to be started with much more solvent or a different arrangement of the solute molecules. as the solute molecules seem to be pushing the box apart and are only slowly rearranging during the npt run (and the box shrinking again).

bottom line. try a generating a better initial geometry, and the simulation will proceed better and equilibration will take less time and effort.

axel.

in.npt (5.45 KB)