Bond atom missing in image check error when using delete_atoms command

Hi all,

I am trying to carve out a droplet from a hydrocarbon simulation. I am defining a region outside the sphere of desired radius, and then deleting the atoms, with complete molecules. The simulation box is 140 nm x 140 nm x 140 nm big. I am able to carve out a sphere of 15 nm radius, but when I try to carve out a bigger liquid sphere, it throws the error, “Bond atom missing in image check”. I have tried increasing the one and page in neigh_modify. Also, tried to use different restart file. I think this might be due to some atoms with non-zero image, but not sure how to rectify this.

Below is the chunk of code to carve out the sphere (This works fine)

region                          NPsphere sphere 250.0 250.0 250.0 157.8
region                          emptySpace sphere 250.0 250.0 250.0 157.8 side out
delete_atoms            region emptySpace compress no mol yes

But

region                          NPsphere sphere 250.0 250.0 250.0 257.8
region                          emptySpace sphere 250.0 250.0 250.0 257.8 side out
delete_atoms            region emptySpace compress no mol yes

does not work.

The complete script is given below for your reference.

units                   real								
dimension               3 									
boundary                p p p 								
atom_style              full 								

comm_modify             vel yes 							
neighbor                1.5 bin 							
neigh_modify            one 300000 page 3000000 			

variable                Nsteps equal 800000					
variable                dt equal 1 						
variable                seed equal 12345 					
variable                T equal 300 						
variable                P equal 7						
pair_style              lj/cut/coul/cut  8.75 25.0 			
#Bonds 														
bond_style              harmonic
#Angles 													
angle_style             harmonic
#Dihedrals 													
dihedral_style          opls 						
#Non-bonded
special_bonds           lj/coul 0.0 0.0 0.5 						

read_restart            ../../step_2_liquid_simulation/liqP7.restart.3000000



pair_coeff    1    1  0.065999     3.500000  # C00 C00
pair_coeff    2    2  0.065999     3.500000  # C01 C01
pair_coeff    3    3  0.065999     3.500000  # C02 C02
pair_coeff    4    4  0.065999     3.500000  # C03 C03
pair_coeff    5    5  0.065999     3.500000  # C04 C04
pair_coeff    6    6  0.065999     3.500000  # C05 C05
pair_coeff    7    7  0.065999     3.500000  # C06 C06
pair_coeff    8    8  0.065999     3.500000  # C07 C07
pair_coeff    9    9  0.065999     3.500000  # C08 C08
pair_coeff   10   10  0.065999     3.500000  # C09 C09
pair_coeff   11   11  0.065999     3.500000  # C0A C0A
pair_coeff   12   12  0.065999     3.500000  # C0B C0B
pair_coeff   13   13  0.030000     2.500000  # H0C H0C
pair_coeff   14   14  0.030000     2.500000  # H0D H0D
pair_coeff   15   15  0.030000     2.500000  # H0E H0E
pair_coeff   16   16  0.026291     2.500000  # H0F H0F
pair_coeff   17   17  0.026291     2.500000  # H0G H0G
pair_coeff   18   18  0.026291     2.500000  # H0H H0H
pair_coeff   19   19  0.026291     2.500000  # H0I H0I
pair_coeff   20   20  0.026291     2.500000  # H0J H0J
pair_coeff   21   21  0.026291     2.500000  # H0K H0K
pair_coeff   22   22  0.026291     2.500000  # H0M H0M
pair_coeff   23   23  0.026291     2.500000  # H0N H0N
pair_coeff   24   24  0.026291     2.500000  # H0O H0O
pair_coeff   25   25  0.026291     2.500000  # H0P H0P
pair_coeff   26   26  0.026291     2.500000  # H0Q H0Q
pair_coeff   27   27  0.026291     2.500000  # H0R H0R
pair_coeff   28   28  0.026291     2.500000  # H0S H0S
pair_coeff   29   29  0.026291     2.500000  # H0T H0T
pair_coeff   30   30  0.026291     2.500000  # H0U H0U
pair_coeff   31   31  0.026291     2.500000  # H0V H0V
pair_coeff   32   32  0.026291     2.500000  # H0W H0W
pair_coeff   33   33  0.026291     2.500000  # H0X H0X
pair_coeff   34   34  0.026291     2.500000  # H0Y H0Y
pair_coeff   35   35  0.026291     2.500000  # H0Z H0Z
pair_coeff   36   36  0.030000     2.500000  # H10 H10
pair_coeff   37   37  0.030000     2.500000  # H11 H11
pair_coeff   38   38  0.030000     2.500000  # H12 H12
pair_coeff   39   39  0.0739378    3.29155   # N   N

group                   Catoms type 1 2 3 4 5 6 7 8 9 10 11 12 	
group                   Hatoms type 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38								

thermo_style            custom step etotal ke pe press vol density temp 		
thermo                  10 												


timestep                ${dt}
fix                     rigidbonds all shake 0.0001 20 0 b 2 3 4 6 7 9 10 12 13 15 16 18 19 21 22 24 25 27 28 30 31 33 34 35 36 37
fix                     liq all npt temp ${T} ${T} $(50*dt) iso ${P} ${P} $(2000*dt)

run                     1 # Ran for one step to see if the simulation is ready to run; it runs fine for this step

region                          NPsphere sphere 250.0 250.0 250.0 257.8
region                          emptySpace sphere 250.0 250.0 250.0 257.8 side out
delete_atoms            region emptySpace compress no mol yes

run 0

write_data              liqP7.50nm.data

Thanks in advance. I am using LAMMPS (2 Aug 2023 - Update 2) in NERSC.

I have two suggestions.

  1. you can try to “minimize” the image flag values with the reset_atoms command
  2. I suspect (but I cannot tell for certain) that your droplet is too large for the simulation box so that the same molecule may be cut from both sides due to periodic boundary conditions. So try to enlarge the system first using “replicate 2 2 2” before trying to cut out a sphere.

Why don’t you use bond yes option? I’m pretty sure bonded atoms don’t have to have the same molID, so mol yes won’t work for you.