Bond Missing Error When Simulating N₂ in Flexible MOF with LAMMPS

LAMMPS LAMMPS General Discussion
, ,
1

I am currently encountering an issue when simulating nitrogen (N₂) molecules inside a flexible MOF using LAMMPS.
The error message I receive is: ERROR on proc 0: Bond atoms 4897 4898 missing on proc 0 at step 1 (src/ntopo_bond_all.cpp:61)
Interestingly, when I apply the same setup (same MOF structure and simulation conditions) to CO₂ molecules, the simulation runs normally without errors. Therefore, I suspect the issue may be related to the N₂ configuration.

I have already tried the following approaches:

  1. Changing the initial configuration of N₂ molecules.
  2. Adjusting the neighbor settings (increasing skin distance).
  3. Disabling the minimization step before dynamics.
    However, none of these attempts resolved the issue.
    Could you please provide advice or suggestions on what might be causing this problem and how I could fix it?
    Thank you very much for your help!

input:
variable RandomSeed equal 720
variable Ndump equal 2000
variable Temp equal 308.0
variable Tdamp equal 100.0

----Initialization----

atom_style full
units real
dimension 3
boundary p p p
box tilt large

----Potential----

pair_style lj/cut/coul/long 12.0
pair_modify shift yes mix arithmetic
bond_style harmonic
angle_style hybrid cosine/periodic harmonic
dihedral_style harmonic
improper_style fourier

dielectric 1.0

special_bonds lj/coul 0.0 0.0 1.0

#flexiblemof
read_data NiMOF.data extra/atom/types 2 extra/bond/types 1 extra/angle/types 1
#N2
read_data N4.data add append offset 4 7 12 3 1 shift 0 0 0 #atom,bond,angle,dihedral,improper
kspace_style pppm 1.0e-6

----------------- Settings Section -----------------

#5 N_N2
#6 X_N2
pair_coeff 5 5 0.0715 3.310
pair_coeff 6 6 0 0

bond_coeff 8 100 0.55 #²Ä¤@­Óµ¹«Ü¤j the second is the bond length

angle_coeff 1 cosine/periodic 132.627486 1 4 # O_2 Ni4+2 O_2
angle_coeff 2 cosine/periodic 82.500443 -1 3 # Ni4+2 O_2 C_R
angle_coeff 3 cosine/periodic 40.148109 -1 3 # Ni4+2 O_2 H_
angle_coeff 4 cosine/periodic 68.268875 -1 3 # C_R O_2 H_
angle_coeff 5 cosine/periodic 87.183054 -1 3 # Ni4+2 O_2 C_R
angle_coeff 6 cosine/periodic 67.154061 -1 3 # Ni4+2 O_2 Ni4+2
angle_coeff 7 cosine/periodic 139.329016 -1 3 # O_2 C_R C_R
angle_coeff 8 cosine/periodic 111.297508 -1 3 # C_R C_R C_R
angle_coeff 9 cosine/periodic 57.289016 -1 3 # C_R C_R H_
angle_coeff 10 cosine/periodic 102.183280 -1 3 # C_R C_R C_R
angle_coeff 11 cosine/periodic 206.778426 -1 3 # O_2 C_R O_2
angle_coeff 12 cosine/periodic 137.714246 -1 3 # O_2 C_R C_R
angle_coeff 13 harmonic 1000 180.0

#group MOF id 1:4896
group MOF type 1 2 3 4
group N2 type 5 6

neighbor 6.0 bin
neigh_modify every 1 delay 0 check yes

timestep 0.25

------------- Setup ----------------

compute myTemp all temp
thermo_style custom step etotal pe ke evdwl ecoul epair ebond eangle edihed eimp elong etail temp
thermo_modify flush yes temp myTemp #print¥X mytemp
thermo 500
#----------------- Minimize Section -----------------
dump whole_minimization all xyz 1000 minimization.xyz # Dump every 1000 steps
dump_modify whole_minimization sort id # Sort atoms by ID for better visualization
min_style cg
minimize 1.0e-4 1.0e-6 10000 100000
minimize 1.0e-6 1.0e-8 10000 100000
undump whole_minimization

----------------- Run Section -----------------

velocity all create {Temp} {RandomSeed} dist gaussian
run 0
velocity all scale ${Temp}

– run at constant volume (Nose-Hoover) –

fix 10 all nvt temp 308.0 308.0 100.0

dump dump_adsorbate N2 custom 1 N2_1.xyz id mol type x y z ix iy iz xs ys zs xu yu zu
dump_modify dump_adsorbate element Ni O C H N_N2 X_N2 sort id

dump dump_mof MOF custom 1 NiMOF_1.xyz id mol type x y z
dump_modify dump_mof element Ni O C H N_N2 X_N2 sort id

run 10 (for testing)
#run 100000000

print “nvt end”

2

Originally, in my N₂ model, there is a massless site whose mass was set to approximately 1e-15.
However, when I changed its mass to 0.1, LAMMPS was able to run without crashing.
I am not sure whether this adjustment is physically appropriate or not.

3

First off, please read our guidelines and suggestions post to learn about how to properly post and specifically how to correct quote text (using triple backticks “```”). This way your input file will not be such a difficult to read mess. Also, the post recommends to search the archives and the error you have seen has been discussed many, many times, so there is a lot of useful information available to you. You just need to look it up.

As has been pointed out in previous discussions (and should be obvious from text books on MD and numerical integration of the EOMs), the smallest mass in the system determines how large a timestep can be. If the mass is as small as 10^{-15}, then the corresponding timestep should be very small. With a massless site it would be impossibly small. Thus massless sites can only be used in a meaningful way if either you have a rigid body, or you have a special pair style can compute the location of the massless site on-the-fly and project the forces on it to the atoms with a mass (like for the TIP4P water model styles).

Whether a simulation run completes without crashing is not any proof that the simulation is correct and physically meaningful.

If you want to do meaningful MD simulations, that is a question that you need to be able to answer yourself.

5

Thanks alot! BTW, I am simulating rigid N2 with massless sites. I have seen several posts talking about this. I tried making the timesteps smaller and assign initial velocity to only those with masses, however, it still doesn’t work if I set the mass of the massless site to be as small as 1e-15.

6

@chenyunhan Since you have not addressed my comments in my previous post, there is nothing else that I can do here. If you don’t make it easy enough for people to help you, you are not likely getting help.

I have run simulations of rigid objects with “fictitious” particles having masses as little as 10^{-100} mass units. So if you are unable to do that, there is probably something wrong in what you do.