ERROR on proc 29: Bond atoms 56578 59390 missing on proc 29 at step 35369

LAMMPS LAMMPS Beginners
1

I’ve been trying to equilibrate my system, but I keep getting an error saying “bond missing.” I’ve checked the data file multiple times, and everything looks fine. The problem occurs when the bond breaks during the dynamics simulation. I’ve tried using both 1 fs and 0.5 fs time steps, but I’m still getting the same error. I’ve been stuck on this for a while now, so any help or suggestions would be really helpfull my input is # ============================================================

FINAL PRODUCTION SCRIPT (FIXED SHAKE)

============================================================

— 1. Initialization —

units real
atom_style full
boundary p p p
newton on

— 2. Force Field Styles —

pair_style lj/charmmfsw/coul/long 10.0 12.0
pair_modify mix arithmetic
special_bonds charmm

bond_style harmonic
angle_style hybrid charmm cosine/periodic fourier
dihedral_style hybrid charmmfsw harmonic
improper_style hybrid fourier harmonic

— 3. Read Data —

read_data final_minimized.data
include combined_params_cleaned.in
kspace_style pppm 1e-6

— 4. Group Definitions —

group enzyme id 1:5604
group sol id 5605:52320
group mof id 52321:59394

fix safe_shake all shake 0.0001 20 0 b 60 110 a 1

Create a group for the solid parts (MOF + Enzyme)

#group nonsol union enzyme mof

— 5. System Settings —

neigh_modify delay 0 every 1 check yes
comm_style tiled
fix bal all balance 1000 1.1 rcb

— Restraints —

variable K_enzyme equal 2.0
variable K_mof equal 0.5
fix restrain_enzyme enzyme spring/self {K_enzyme} fix restrain_mof mof spring/self {K_mof}

— Velocity Init —

velocity all create 10.0 54321 mom yes rot yes dist gaussian
fix mom_remove all momentum 1000 linear 1 1 1

============================================================

STAGE 0: Soft Warmup (NVE/Limit)

============================================================

timestep 0.5
fix limiter all nve/limit 0.1
fix lang_warm all langevin 10.0 50.0 100.0 699483

thermo 100
run 5000

unfix limiter
unfix lang_warm

============================================================

STAGE A: Heating (NVT)

============================================================

timestep 1.0

Resume from 50K

fix stage_A all nvt temp 50.0 300.0 100.0

thermo 1000
thermo_style custom step temp pe etotal press vol density
dump dump_A all custom 5000 dump_eq_A_heat.lammpstrj id type xu yu zu

run 200000 # 200 ps

write_restart restart.stage_A
unfix stage_A
undump dump_A

============================================================

STAGE B: Relaxation

============================================================

fix stage_B all nvt temp 300.0 300.0 100.0
dump dump_B all custom 5000 dump_eq_B_relax.lammpstrj id type xu yu zu

run 500000 # 0.5 ns

write_restart restart.stage_B
unfix stage_B
undump dump_B

============================================================

STAGE C: Production

============================================================

unfix restrain_enzyme
unfix restrain_mof

fix stage_C all nvt temp 300.0 300.0 100.0
dump dump_C all custom 10000 dump_production.lammpstrj id type xu yu zu

run 1000000 # 1 ns

write_restart final_production.restart
write_data final_production.data

2

Hi @droy,

It is impossible to provide meaningful help with so few information. Some comments though:

  1. Please read this post. It will help you formatting your question and give you some steps into understanding where your error comes from.

  2. When you say:

    The problem occurs when the bond breaks during the dynamics simulation.

    This means that either the bond break is expected and then your model is wrongly set (harmonic type bonds are not supposed to break), either your not expecting such breaks and then you simulation parameters and methodology are likely wrong.

3

Thank you for your time. I am working on a project where I’m inserting an enzyme into a MOF pore. To prepare the system, I generated the force field for the MOF using the LAMMPS-Interface and the force field for the enzyme using CHARMM-GUI. Afterward, I used VMD to insert the enzyme into the MOF and renumbered the atoms: atoms 1-57 are assigned to the enzyme, and atoms 58-66 correspond to the MOF. I also updated the coefficients accordingly. I’ve reviewed everything, and it seems correct, but as a new user, I am unable to upload the files here for further inspection. Is there any way I could share the files with you so you can help me identify any potential issues?

4

I don’t know any tool called “LAMMPS-interface” that allow to generate a force field so I don’t know what you are talking about.

There are always many services that allow to do so and from which you can post a link to your files. But this is something that would gain more from discussing with your advisor, colleagues or people in your lab who actually know about your scientific problem and topic. You are more likely to gain long term learning about your problem at hand this way than looking for answers from strangers on the internet.

5

Did you also update the atom IDs for bonds, angles, dihedrals, etc. accordingly?

6

Yes, I made all the changes. I have uploaded all the relevant files here:

https://drive.google.com/drive/folders/1IC0lI3McNv1V9N5GWqZ6pqlThjUdWb2C?usp=drive_link

I also modified one line in the test.in LAMMPS input file under the System Settings section:

neigh_modify delay 0 every 1 check yes page 100000 one 10000
neighbor     2.0 bin

However, after this change I am now getting the following error:

ERROR on proc 35: Dihedral atoms 3649 3651 3657 3659 missing on proc 35 at step 38469 (src/ntopo_dihedral_all.cpp:64)
Last command: run 200000   # 200 ps

I am unable to determine the cause of this issue. Could you please help me figure out what might be happening?

7

As far as I can tell the issue originates from fix shake or rather not using it.
Without fix shake, a timestep of 0.5fs and especially 1.0fs is not stable.
Also, your system appears to be of high initial potential energy. Using fix nve/limit is not the best option in that case, better is to do a minimization first. Using fix shake is allowed with minimization provided you have a sufficiently recent LAMMPS version. It will apply some strong harmonic restraint forces during minimization. I modified your input accordingly and also noted that you may have set the shake angle type incorrectly. The water angle type appears to be 148 and not 1. Please find below the updated first part of the test.in file for evaluation and further testing:

# --- 1. Initialization ---
units           real
atom_style      full
boundary        p p p
newton          on

# --- 2. Force Field Styles ---
pair_style      lj/charmmfsw/coul/long 10.0 12.0
pair_modify     mix arithmetic
special_bonds   charmm

bond_style      harmonic
angle_style     hybrid charmm cosine/periodic fourier
dihedral_style  hybrid charmmfsw harmonic
improper_style  hybrid fourier harmonic

# --- 3. Read Data ---
read_data       final_minimized_again.data
include       combined_params_cleaned.in  
kspace_style    pppm 1e-6

# --- 4. Group Definitions ---
group enzyme id 1:5604
group sol    id 5605:52320
group mof    id 52321:59394

fix safe_shake all shake 0.0001 20 1000 b 60 110 116 a 148

# Create a group for the solid parts (MOF + Enzyme)
#group nonsol union enzyme mof

# --- 5. System Settings ---
neigh_modify    delay 0 every 1 check yes
neighbor        2.0 bin

# --- Restraints ---
variable K_enzyme equal 2.0
variable K_mof    equal 0.5
fix restrain_enzyme enzyme spring/self ${K_enzyme}
fix restrain_mof    mof    spring/self ${K_mof}

# --- Velocity Init ---
velocity all create 10.0 54321 mom yes rot yes dist gaussian
fix mom_remove all momentum 1000 linear 1 1 1

thermo 10
minimize 0.00001 0.000001 100 1000

reset_timestep 0

# ============================================================
#   STAGE 0: Soft Warmup
# ============================================================
timestep 0.5
fix limiter all nve
fix lang_warm all langevin 10.0 50.0 100.0 699483

thermo 100
run 20000    # <--- Increased duration

unfix limiter
unfix lang_warm
8

yes I tried this and it running simulation but only in 0.5 fs if i go to 1fs that will fail…thank you for your time.