ERROR on proc 20: Bond atoms 341 342 missing on proc 20 at step 4459285 (../ntopo_bond_all.cpp:63)

Hi all,

When I run LAMMPS using QEQ for ionic liquid systems with 200 ion-pairs, I always encounter this kind of error after several ns. Can you suggest how to deal with this problem? My input file is:

# Number of atoms = 9600; Total charge = 6.953437825529818e-12
units real
boundary p p p

atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls

special_bonds lj/coul 0.0 0.0 0.5

# remove hybrid if not necessary
pair_style lj/cut/coul/long 12.0 12.0
pair_modify mix geometric tail yes
kspace_style ewald 1.0e-5

read_data data-box.lmp

pair_coeff    1    1      0.170000     3.250000  # N1 N1
pair_coeff    2    2      0.065999     3.500000  # C2 C2
pair_coeff    3    3      0.065999     3.500000  # C3 C3
pair_coeff    4    4      0.065999     3.500000  # C4 C4
pair_coeff    5    5      0.065999     3.500000  # C5 C5
pair_coeff    6    6      0.065999     3.500000  # C6 C6
pair_coeff    7    7      0.065999     3.500000  # C7 C7
pair_coeff    8    8      0.065999     3.500000  # C8 C8
pair_coeff    9    9      0.140000     2.900000  # O9 O9
pair_coeff   10   10      0.065999     3.500000  # C10 C10
pair_coeff   11   11      0.065999     3.500000  # C11 C11
pair_coeff   12   12      0.030000     2.500000  # H12 H12
pair_coeff   13   13      0.030000     2.500000  # H13 H13
pair_coeff   14   14      0.030000     2.500000  # H14 H14
pair_coeff   15   15      0.030000     2.500000  # H15 H15
pair_coeff   16   16      0.030000     2.500000  # H16 H16
pair_coeff   17   17      0.030000     2.500000  # H17 H17
pair_coeff   18   18      0.030000     2.500000  # H18 H18
pair_coeff   19   19      0.030000     2.500000  # H19 H19
pair_coeff   20   20      0.030000     2.500000  # H20 H20
pair_coeff   21   21      0.030000     2.500000  # H21 H21
pair_coeff   22   22      0.030000     2.500000  # H22 H22
pair_coeff   23   23      0.030000     2.500000  # H23 H23
pair_coeff   24   24      0.030000     2.500000  # H24 H24
pair_coeff   25   25      0.030000     2.500000  # H25 H25
pair_coeff   26   26      0.030000     2.500000  # H26 H26
pair_coeff   27   27      0.030000     2.500000  # H27 H27
pair_coeff   28   28      0.030000     2.500000  # H28 H28
pair_coeff   29   29      0.030000     2.500000  # H29 H29
pair_coeff   30   30      0.030000     2.500000  # H30 H30
pair_coeff   31   31      0.030000     2.500000  # H31 H31
pair_coeff   32   32      0.030000     2.500000  # H32 H32
pair_coeff   33   33      0.030000     2.500000  # H33 H33
pair_coeff   34   34      0.250000     3.550000  # S1 S1
pair_coeff   35   35      0.170000     3.250000  # N2 N2
pair_coeff   36   36      0.170000     2.960000  # O3 O3
pair_coeff   37   37      0.170000     2.960000  # O4 O4
pair_coeff   38   38      0.060000     2.900000  # F6 F6
pair_coeff   39   39      0.060000     2.900000  # F7 F7
pair_coeff   40   40      0.060000     2.900000  # F8 F8
pair_coeff   41   41      0.250000     3.550000  # S9 S9
pair_coeff   42   42      0.170000     2.960000  # O10 O10
pair_coeff   43   43      0.060000     2.900000  # F12 F12
pair_coeff   44   44      0.060000     2.900000  # F13 F13
pair_coeff   45   45      0.060000     2.900000  # F14 F14
pair_coeff   46   46      0.170000     2.960000  # O15 O15


minimize 1.0e-4 1.0e-6 100 1000
reset_timestep 0

#fix SHAKE all shake 0.0001 20 0

neighbor 2.0 bin
neigh_modify delay 0 every 5 check yes page 1000000 one 20000

timestep 1.0

# QEQ: 2ps * 5000 = 10000 ps=10ns
fix 1 all qeq/point 10 10 1.0e-6 5000 param.qeq1

####################################################
#                                                  #
#            1.       NH:  NPT                     #
#                                                  #
####################################################
variable PBAR equal 1.0

variable TK equal 400.0
velocity all create ${TK} 12345

fix TPSTAT all npt temp ${TK} ${TK} 100 iso ${PBAR} ${PBAR} 1000

thermo_style custom step cpu etotal ke pe evdwl ecoul elong temp press vol density
thermo 5000

# monitor charge
dump TRAJ all custom 10000 eq.lammpstrj id mol type element q xu yu zu
dump_modify TRAJ element N C C C C C C C O C C H H H H H H H H H H H H H H H H H H H H H H S N O O F F F S O F F F O

# NPT 50 ps at 400 K
run 50000

unfix TPSTAT
variable TK equal 298.15
velocity all scale ${TK}

fix TPSTAT all npt temp ${TK} ${TK} 100 iso ${PBAR} ${PBAR} 1000

# NPT 100 ps at 298.15 K
run 100000

unfix TPSTAT
fix TSTAT all nvt temp ${TK} ${TK} 100

# NVT 300 ps
run 300000

write_data eq.data
write_restart eq.restart

####################################################
#                                                  #
#          3. NVT: Viscosity & MSD & RDF           #
#                                                  #
####################################################

variable    T equal 298.15
variable    V equal vol
variable    dt equal 1.0

# Maxdt=1 ns=2000000 fs, nunmber of trajectories = 10 ns/ 1 ns = 10
variable    s equal 10       # sample interval
variable    p equal 500000   # number of correlation time windows
variable    d equal 5000000   # dump interval 2 ns

variable    kB equal 1.3806504e-23    # [J/K] Boltzmann
variable    atm2Pa equal 101325.0
variable    A2m equal 1.0e-10
variable    fs2s equal 1.0e-15
variable    convert equal ${atm2Pa}*${atm2Pa}*${fs2s}*${A2m}*${A2m}*${A2m}

reset_timestep 0

variable     pxy equal pxy
variable     pxz equal pxz
variable     pyz equal pyz
# Accumulate samples & Calculate autocorrelations for every dt
fix          SS all ave/correlate $s $p $d &
             v_pxy v_pxz v_pyz type auto file S0St.dat ave running
# Perform the GK running integral
variable     scale equal ${convert}/(${kB}*$T)*$V*$s*${dt}
variable     v11 equal trap(f_SS[3])*${scale}
variable     v22 equal trap(f_SS[4])*${scale}
variable     v33 equal trap(f_SS[5])*${scale}

variable     v equal (v_v11+v_v22+v_v33)/3.0
variable     ndens equal count(all)/vol

timestep     ${dt}

fix extra1 all ave/time 10 1 10 v_pxy v_pxz v_pyz file pxyz.txt
fix extra2 all ave/time $d 1 $d v_v11 v_v22 v_v33 file v123.txt

thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33 v_v
thermo   $d

restart 2000000 meta.inbetween.restart meta.inbetween.restart
run 10000000    # 10 ns

write_data compute.data
write_restart compute.restart

print        "average viscosity: $v [Pa.s] @ $T K, ${ndens} /A^3"

Thank you very much!

I confess that reading this already in bed at 22h, without my contact lenses, I am not successfully seeing a straightfoward cause for the simution to crash within my not-so-exactly-particular-high knowledge on LAMMPS either.

When exactly does it crash? If it is already during the first NPT? Maybe iit could be worth it to save a microstate after the minimization to see the looks of it, like for example if somehow it is diverging to something that would visually look akward (and then maybe you would find a typo in the potential, idk, although I suppose you checked).

Would it run well without the fix qeq/point command? I dont know this command bu from reading here quicikly in the manual it seems to be challanging to choose parameters for it, so maybe the parameters are not optimum somehow (there seems to be an example of a tutorial using this fix so maybe you can try it and see if it works).

1 Like

These kind of things most commonly happen, when your force field is not suitable for the system or when the simulation settings like size of timestep are not suitable.

For example a timestep of 1fs seems rather aggressive for a system with explicit hydrogen atoms and no bond constraints. Have you checked your energy conservation when running with fix nve? I would suspect it is not very good, which would confirm the timestep length as one issue. At the same time, updating the charges only every 10 steps seems a bit aggressive for such a large timestep. But then again, this may hide the effect of unsuitable force field and polarization parameters to some degree.

I haven’t seen QEQ used with bonded molecules before, and it seems like a point QEQ for ionic liquids would be especially prone to catastrophic polarisation (when a positive charge and a negative charge are near each other and the charge can vary limitlessly, the system energy can be “minimised” by giving both charges infinite magnitude).

You can track this with something like

variable qsq atom q*q
compute qsqmax all reduce max v_qsq

and then adding c_qsqmax to your thermo output. But really, I’m not sure about QEQ in this context. Are you expecting significant charge transfer between the cations and anions?

EDIT: also, that’s an awful lot of atom types with identical parameters that you should be able to consolidate from chemical reasoning. For example, if you’re simulating a triflate-based anion, all of those fluorines are chemically identical and you should only need one F type, not 6. Technically there’s nothing wrong with lots and lots of atom types but the more you have, the more likely you are to make a silly parameter swap (your simulation might run a tiny bit faster with fewer types too).

1 Like

Thank you for the help!

I tried updating the charges every 50, 100, and 200 steps, but WARNING: H matrix size has been exceeded poped up.

I also tried continue the simulation from restart, and sometimes the simulation could finish normally, while sometimes the same error appeared again at some timestep.

I was expecting that you would do the charge equilibration more often, not less often.
What is the setting used in the publication that you took your force field parameters from?

I just tried fix 1 all qeq/point 1 10 1.0e-6 5000000 param.qeq1. The error occurred at step 4756480. Then I read_restart from step 4000000, the simulation could run until the error occurred again at step 8606155. It seems this error occurs by chance.

My force field parameters are from the LigPargen server (http://zarbi.chem.yale.edu/ligpargen/). Unfortunately, I cannot find a publication that use qeq on my system so I am trying to figure out by myself.

So how do you justify the use of charge equilibration in the first place? I don’t think that OPLS-AA parameters are parameterized for that. So your problems may be self-inflicted.

Have you run your system without fix qeq/point? Does it still crash as before?

It would not crash without fix qeq/point.

Charge equilibration works by solving a sparse matrix equation across the entire simulation system (for the targetted electronegativities in terms of the desired charges). Whichever algorithm is used, the result will depend on multiple sums with a large number of terms gathered across all the MPI processes used by LAMMPS, which will induce floating point variation based on the order of summation. Since the charges are then used to calculate forces (as well as the next round of charges), any floating point error accumulates far more quickly than in a comparable fixed-charge simulation.

Again, why exactly are you trying to use QEq for ionic liquids, to say nothing of combining it with OPLS LJ parameters (and therefore making an invalid force field until proven otherwise, because a force field is specified by all its force components at once)?

Thank you very much for the help. Basically, I want to include polarization effects, so I started with QEq. I think it is the most efficient approach compared with polarizable force fields or the Drude model.