Repeated Errors and Truncation of simulation during Production run

LAMMPS
1

Hello, I really need an help in my simulation (currently I am doing an MD simulation related to Gel polymer electrolyte and I am studying the ion transport and Structural dynamics of Gel polymer electrolyte). I am facing the repeated problems of “bond atom missing” and “out of range atoms - cannot compute pppm” during my production run even when my both the equilibration (NVT and NPT) has been completely successfully , can you give me some solutions so as to avoid these errors.
I am pasting my job script with which I am running the simulation:-

# ============================================================
# Polymer + EMIM + TFSI
# NVT (200 ps) → NPT equil (500 ps) → NPT production (100 ns)
# With density, MSD, RDF, conductivity, Rg
# ============================================================

units           real
atom_style      full
dimension       3
boundary        p p p

bond_style      harmonic
angle_style     harmonic
dihedral_style  opls
improper_style  cvff

# ---- GPU (pair) + CPU (PPPM) ------------------------------------
package         gpu 1
suffix          off
pair_style      lj/cut/coul/long/gpu 12.0 12.0
kspace_style    pppm 1.0e-6
suffix          gpu

read_data       mylammps.data
include         parameters_new.inc

# ---- Neighbor / communication -----------------------------------
neighbor        4.0 bin
neigh_modify    delay 0 every 1 check yes one 8000
comm_modify     cutoff 20.0 vel yes

# ---- Groups ------------------------------------------------------
group polymer type 1 2 5 6 8 12 14
group EMIM    type 4 7 9 10 11 15 16 17 18 19 20
group TFSI    type 3 13 21 22 23
group ions subtract all polymer

# ---- Thermo ------------------------------------------------------
thermo          1000
thermo_style    custom step time temp press vol density pe ke etotal
thermo_modify   flush yes

restart         100000 final_restart.*

# ============================================================
# (1) Energy minimization
# ============================================================
min_style       sd
thermo          500
minimize        1.0e-4 1.0e-5 100000 1000000

min_style       cg
thermo          500
minimize        1.0e-8 1.0e-10 100000 1000000
write_data      minimize.data

# ============================================================
# (2) NVT equilibration: 200 ps @ 1 fs
# ============================================================
timestep        1.0
velocity        all create 298.0 12345 mom yes rot yes dist gaussian loop geom

fix             mynvt all nvt temp 298.0 298.0 100.0
thermo          1000
run             200000
unfix           mynvt

# ============================================================
# (3) NPT equilibration: 500 ps @ 1 fs
# ============================================================
fix             mynpt_eq all npt temp 298.0 298.0 100.0 iso 1.0 1.0 5000.0
thermo          1000
run             500000
unfix           mynpt_eq

# ============================================================
# (4) Production NPT: 100 ns @ 2 fs
# ============================================================

reset_timestep 0

# ---- Dumps -------------------------------------------------------------------------------
# Wrapped coordinates
dump            d_wrap all custom 5000 dump_wrapped.lammpstrj id mol type q x y z
dump_modify     d_wrap sort id

# Unwrapped coordinates (for MSD, diffusion)
dump d_unwrap all custom 5000 dump_unwrapped.lammpstrj id mol type q xu yu zu ix iy iz
dump_modify     d_unwrap sort id
#-------------------------------------------------------------------------------------------

timestep 2.0
fix mynpt_prod all npt temp 298.0 298.0 100.0 iso 1.0 1.0 5000.0

fix nodrift all momentum 100 linear 1 1 1

# ---- STATIC / INSTANTANEOUS COMPUTES -----------------------

variable rho equal density
fix f_rho all ave/time 1000 100 100000 v_rho file density_ave.dat

compute c_rg polymer gyration
fix f_rg all ave/time 1000 100 100000 c_c_rg file rg_polymer.dat

# ---- TIME-DEPENDENT COMPUTES (AFTER FIX) -------------------
set group all image 0 0 0

compute c_msd_poly polymer msd 
compute c_msd_emim EMIM msd
compute c_msd_tfsi TFSI msd

fix f_msd all ave/time 1000 100 100000 c_c_msd_poly[4] file msd_poly.dat
fix f_msd_1 all ave/time 1000 100 100000 c_c_msd_emim[4] file msd_emim.dat
fix f_msd_2 all ave/time 1000 100 100000 c_c_msd_tfsi[4] file msd_tfsi.dat

# ---- VACF --------------------------------------------------------

compute c_vacf_poly polymer vacf
compute c_vacf_emim EMIM vacf
compute c_vacf_tfsi TFSI vacf

fix f_vacf all ave/time 10 1000 10000 c_c_vacf_poly[4] file vacf_poly.dat
fix f_vacf_1 all ave/time 10 1000 10000 c_c_vacf_emim[4] file vacf_emim.dat
fix f_vacf_2 all ave/time 10 1000 10000 c_c_vacf_tfsi[4] file vacf_tfsi.dat

# ---- CONDUCTIVITY (Green–Kubo) -----------------------------

compute c_q all property/atom q
variable jx atom q*vx
variable jy atom q*vy
variable jz atom q*vz

compute c_j all reduce sum v_jx v_jy v_jz

variable Jx equal c_c_j[1]
variable Jy equal c_c_j[2]
variable Jz equal c_c_j[3]

fix f_cond all ave/correlate 10 1000 10000 v_Jx v_Jy v_Jz type auto file charge_current_acf.dat

# ---- SEGMENTAL RELAXATION ----------------------------------

compute c_dih_poly polymer dihedral/local phi

dump d_dih polymer local 5000 dihedral_dump.dat index c_c_dih_poly


thermo 5000
run 50000000


# ============================================================
# ---------------------- END --------------------------------
# ============================================================

So is there any way I can avoid that…

2

Have you looked up previous discussions on these errors and checked the LAMMPS documentation what it suggests to look for with these errors?

If not, please do so, if yes, what suggested remedies have you tried and what was the outcome?

Three additional questions:

  • Is your system an all-atom simulation and thus have explicit hydrogen atoms?
  • Have you checked the energy conservation without a thermostat?
  • How did you compile the GPU package (single, mixed, or double precision)?
3
  1. Yes this is all-atom simulation and have explicit hydrogen atoms in the co-polymer and EMIM+
  2. No, I have not done it.
  3. I have compiled the GPU package with Double Precision one.

I have done following remedies like :-
a) changing bin size
b) changing coulombic radius and van der waals radius { pair_style lj/cut/coul/long/gpu 12.0 12.0}
c) changing the parameters of the forcefield related to dihedrals , angles and pairs.

So, is there any other thing that can be done??

4

It is a bad idea to use a 2fs timestep with explicit hydrogens without constraining those bonds to be rigid. Since hydrogen is so much lighter than the other atoms, its bond stretch vibrations are of much higher frequency and you need a shorter time step to accommodate that.
Something like 0.25fs is a conservative choice, people often get away with 0.5fs but beyond that, energy conservation should be unacceptable and you can get artifacts like you have noticed. That is why testing for energy conservation is a good thing to do early. It can tell you if your choice of timestep, neighbor list updates, and cutoffs is too aggressive.

Neither of these have any impact on the symptoms you describe. Rather more on the efficiency of your calculation. Speaking of that:

Neither of these settings make any sense for your current settings. They only make your calculation slower. If you choose to reduce the timestep instead of adding bond constraints, then you can increase delay to, say, 4, but the savings of that are generally small with check yes.

I don’t understand this. You cannot just make up parameters or change them as you like.

5

Sorry for misunderstanding , I was talking about the changing the parameters of bond coeff, angle coeff and dihedral coeff of all the bond type, angle type and dihedral type in the parameter file with respect to OPLS forcefield I have taken, with checking the forcefield file.

6

Sorry, I don’t understand what you are saying here. Both, grammatically and from the science point of view. There is only one setting in a force field for a given atom, bond, angle, or dihedral type, so there is no changing of parameters possible.