delete_atoms works 3 times but not the 4th time

Dear all

I have here a slightly strange and annoying problem. I want to add the same cluster 4 times to a system consisting of a layer of organic molecules. I have set it up and tried it with one cluster for which it works nicely; all the overlapping atoms from the layer are deleted. It works with 2 clusters, it works with 3 clusters. But suddenly, when I want to delete the molecules that overlap with the 4th cluster, it does not work anymore:

Total # of neighbors = 61245534
Ave neighs/atom = 330.635
Ave special neighs/atom = 11.0702
Neighbor list builds = 0
Dangerous builds = 0
Displacing atoms ...
System init for delete_atoms ...
Neighbor list info ...
   2 neighbor list requests
   update every 5 steps, delay 0 steps, check yes
   max neighbors/atom: 20000, page size: 5000000
   master list distance cutoff = 14
   ghost atom cutoff = 14
   binsize = 7, bins = 58 58 58
Deleted 2816 atoms, new total = 182420
Deleted 2752 bonds, new total = 178093
Deleted 5376 angles, new total = 347604
Deleted 7488 dihedrals, new total = 483987
Deleted 0 impropers, new total = 240
System init for delete_atoms ...
Deleted 3300 atoms, new total = 179120
Deleted 3225 bonds, new total = 174868
Deleted 6300 angles, new total = 341304
Deleted 8775 dihedrals, new total = 475212
Deleted 0 impropers, new total = 240
System init for delete_atoms ...
Deleted 3388 atoms, new total = 175732
Deleted 3311 bonds, new total = 171557
Deleted 6468 angles, new total = 334836
Deleted 9009 dihedrals, new total = 466203
Deleted 0 impropers, new total = 240
System init for delete_atoms ...
Deleted 0 atoms, new total = 175732
Deleted 0 bonds, new total = 171557
Deleted 0 angles, new total = 334836
Deleted 0 dihedrals, new total = 466203
Deleted 0 impropers, new total = 240
Neighbor list info ...
   1 neighbor list requests
   update every 5 steps, delay 0 steps, check yes
   max neighbors/atom: 20000, page size: 5000000
   master list distance cutoff = 14
   ghost atom cutoff = 14
   binsize = 7, bins = 58 58 58

As you can see, the 4th time nothing is deleted even though the cluster does overlap with a lot of molecules.

The input script I'm using:

clear

# ------------------------ INITIALIZATION ----------------------------
units real
dimension 3
atom_style full #charge
boundary p p p

#region supcell block 0 1 0 1 0 1
#replicate 1 1 1

# ------------------------ FORCE FIELDS ------------------------------
pair_style lj/cut/coul/cut 12.0
#kspace_style pppm 10e-6
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style cvff

# ---------- Create Atoms ---------------------
# Box is 400x400x400 Ang^3
read_data 15nmSlab.data group slab extra/atom/types 28 extra/bond/types 28 extra/angle/types 52 extra/dihedral/types 60 extra/improper/types 4
read_data au.data group cluster01 add append offset 3 4 7 7 0
read_data au.data group cluster02 add append offset 10 11 20 22 1
read_data au.data group cluster03 add append offset 17 18 33 37 2
read_data au.data group cluster04 add append offset 24 25 46 52 3

# Set variable for element names of the cluster
variable elmts string "S C H Au Au Au C"

# Move the cluster(s) and the slab to the middle of the box...
compute cofm_slab slab com
compute cofm_cluster01 cluster01 com
compute cofm_cluster02 cluster02 com
compute cofm_cluster03 cluster03 com
compute cofm_cluster04 cluster04 com

# Move the slab's com to the center of the box:
variable trans_slab_x equal "(xhi-xlo)/2-c_cofm_slab[1]"
variable slab_x equal v_trans_slab_x
variable trans_slab_y equal "(yhi-ylo)/2 - c_cofm_slab[2]"
variable slab_y equal v_trans_slab_y
variable trans_slab_z equal "(zhi-zlo)/2 - c_cofm_slab[3]"
variable slab_z equal v_trans_slab_z

# Move cluster #1 to its position:
variable trans_cluster01_x equal "(xhi-xlo)/2 - c_cofm_cluster01[1] + 35"
variable cluster01_x equal v_trans_cluster01_x
variable trans_cluster01_y equal "(yhi-ylo)/2 - c_cofm_cluster01[2]"
variable cluster01_y equal v_trans_cluster01_y
variable trans_cluster01_z equal "(zhi-zlo)/2 - c_cofm_cluster01[3] - 35"
variable cluster01_z equal v_trans_cluster01_z

# Move cluster #2 to its position:
variable trans_cluster02_x equal "(xhi-xlo)/2 - c_cofm_cluster02[1] - 35"
variable cluster02_x equal v_trans_cluster02_x
variable trans_cluster02_y equal "(yhi-ylo)/2 - c_cofm_cluster02[2]"
variable cluster02_y equal v_trans_cluster02_y
variable trans_cluster02_z equal "(zhi-zlo)/2 - c_cofm_cluster02[3] - 35"
variable cluster02_z equal v_trans_cluster02_z

# Move cluster #3 to its position:
variable trans_cluster03_x equal "(xhi-xlo)/2 - c_cofm_cluster03[1]"
variable cluster03_x equal v_trans_cluster03_x
variable trans_cluster03_y equal "(yhi-ylo)/2 - c_cofm_cluster03[2] + 35"
variable cluster03_y equal v_trans_cluster03_y
variable trans_cluster03_z equal "(zhi-zlo)/2 - c_cofm_cluster03[3] + 35"
variable cluster03_z equal v_trans_cluster03_z

# Move cluster #4 to its position:
variable trans_cluster04_x equal "(xhi-xlo)/2 - c_cofm_cluster04[1]"
variable cluster04_x equal v_trans_cluster04_x
variable trans_cluster04_y equal "(yhi-ylo)/2 - c_cofm_cluster04[2] - 35"
variable cluster04_y equal v_trans_cluster04_y
variable trans_cluster04_z equal "(zhi-zlo)/2 - c_cofm_cluster04[3] + 35"
variable cluster04_z equal v_trans_cluster04_z

# ------------------------- SETTINGS ---------------------------------
restart 1000 restart.1.Au-in-C14 restart.2.Au-in-C14
neighbor 2.0 bin
neigh_modify delay 0 every 5 check yes one 20000 page 5000000

# ---------------- Evaluate the displacement vector for all systems

Whatever the problem is, it is not subtle and the cause should be easy to spot once you get a good look at it. I don’t see anything obviously wrong with the script, except for the fact that the system size is so gigantic. The quickest way to track down this problem (and probably a few others too), would be to reduce the calculation down to a few hundred atoms. That would allow you to do a lot of quick tests, which would eliminate a lot of guesswork. But if I had to make a guess, I would say that cluster 4 somehow has the same position as cluster 1, 2, or 3.

Aidan

Well, subtle is only the first name... I had to change this part:

# ---------------- Equilibrate the slab first ---------------------

fix slab_eq slab nvt temp 50.0 298.0 100.0

thermo 1
thermo_style custom step temp press etotal

dump slabtrjeq slab custom 1 slab_only_init.lammpstrj id mol element xu yu zu
dump_modify slabtrjeq append yes element C C H \{elmts\} {elmts} \{elmts\} {elmts} sort id

timestep 1.0
run 10
unfix slab_eq
undump slabtrjeq

displace_atoms cluster01 move v_cluster01_x v_cluster01_y v_cluster01_z units box
displace_atoms cluster02 move v_cluster02_x v_cluster02_y v_cluster02_z units box
displace_atoms cluster03 move v_cluster03_x v_cluster03_y v_cluster03_z units box
displace_atoms cluster04 move v_cluster04_x v_cluster04_y v_cluster04_z units box

# Delete molecules that overlap with the clusters ...
delete_atoms overlap 2.5 slab cluster01 bond yes mol yes
delete_atoms overlap 2.5 slab cluster02 bond yes mol yes
delete_atoms overlap 2.5 slab cluster03 bond yes mol yes
delete_atoms overlap 2.5 slab cluster04 bond yes mol yes

to this:

# ---------------- Evaluate the displacement vector for all systems

Very mysterious. If you send a small demo, we can probably track down the cause.