Delete_atom command is eating up all my polymer atoms

Dear Lammps User,

So the thing is, In my input deck I have initialized some groups with some number of atoms in it. After just running the the that part in LAMMPS, just to see whether the number of atoms in the groups are same or not, I saw some anomaly between what I executed and what I expected. I am not able to figure out why this is happening.

Below is the snippet of input deck:

#************************************************
# SETTING UP BOUNDARY CONDITIONS
#****************************************************

units		    lj
atom_style	    full
dimension       3
boundary     	p p p
#****************************************************
#SETTING UP POTENTIALS
#****************************************************
pair_style      lj/cut/coul/cut 2.5 9.0
angle_style     harmonic
bond_style      harmonic 
neigh_modify    every 1 delay 5 check yes
special_bonds   coul 0.0 0.0 0.5
dielectric     1.0
#****************************************************
#DEFINING REGION &&&& READING DATA FILE
#****************************************************
read_data    data3.txt extra/atom/types 1 extra/bond/types 1 extra/angle/types 1
region          reg1 cylinder z 25.0 25.0 7.0 -0.5 0.5
delete_atoms region reg1 bond yes
read_data polymer.txt add append offset 1 1 1 0 0
region myreg1 block 0 50 0 50 -100 -10.0 #This region is defined to use for create_atoms 
region myreg2 block 0 50 0 50  10.0  100 #This region is defined to use for create_atoms
create_atoms 2 random 200 341341 myreg1
create_atoms 3 random 200 127569 myreg1
create_atoms 2 random 200 341341 myreg2
create_atoms 3 random 200 127569 myreg2
set type 2 charge 0.5
set type 3 charge -0.5
include	     wallparam.lammps
include         ionsparam.lammps
include         polyparam.lammps
pair_modify  shift yes mix geometric
reset_ids  

#****************************************************
#DEFINING GROUPS
#****************************************************
group gcar type 1 #This is group for graphene sheet which contains 1008 atoms but since 
#I used delete_atoms command so there must atoms less than 1008 atoms
group gpot type 2 #This group for monovalent +ve atoms, say, Potassium. There are 400 atoms here.
group gchl type 3 #This group for monovalent -ve atoms, say, Chloride. There are 400 atoms here
group gpoly type 4 #This group is for polymer chain. There are 20 atoms in here.

#Since I have to fix the graphene sheet in such a manner that it should stay exactly in the middle of the simulation box with 
#its dimension being equal to that of  x,y dimensions of the box.
variable xmax equal bound(gcar,xmax)+1.0
variable xmin equal bound(gcar,xmin)-1.0
variable ymax equal bound(gcar,ymax)+1.0
variable ymin equal bound(gcar,ymin)-1.0

variable xmax1 equal bound(gcar,xmax)-0.5
variable xmin1 equal bound(gcar,xmin)+0.5
variable ymax1 equal bound(gcar,ymax)-0.5
variable ymin1 equal bound(gcar,ymin)+0.5

#This part of the code disects the graphene sheet in 4 parts. I'll explain one and the reasoning for the other 3 is similar.
# say rright, a region which is formed on the right of graphene sheet, it includes only the edge atoms on that side of the graphene
# sheet. The idea is that the 4 edges will have setforce 0 only the mid section will be equilibrated.
region rright block ${xmax1} INF INF INF -0.5 0.5
region rleft block INF ${xmin1} INF INF -0.5 0.5
region rmid block ${xmin1} ${xmax1} ${ymin1} ${ymax1} -0.5 0.5
region rtop block INF INF ${ymin1} INF -0.5 0.5
region rbot block INF INF INF ${ymin1} -0.5 0.5

group gright region rright
group gleft region rleft
group gmid region rmid
group gions union gpot gchl
group gtop region rtop
group gbot region rbot
group gall union gions gpoly

change_box all x final ${xmin} ${xmax}
change_box all y final ${ymin} ${ymax}

the data3.txt is

LAMMPS data file. CGCMM style. atom_style full generated by VMD/TopoTools v1.8 on Wed Sep 28 11:43:09 IST 2022
 1008 atoms
 1467 bonds
 2846 angles
 
 3 atom types
 1 bond types
 1 angle types

-6.04 53.94  xlo xhi
-4.86 54.44  ylo yhi
 -50.000000 50.000000  zlo zhi

# Pair Coeffs
#
# 1  CA

# Bond Coeffs
#
# 1  CA-CA

# Angle Coeffs
#
# 1  CA-CA-CA


 Atoms # full

1 1 1 0.000000 0.000000 0.000000 0.000000 # CA GRA
2 2 1 0.000000 -1.228024 0.709000 0.000000 # CA GRA
3 3 1 0.000000 -1.228024 2.127000 0.000000 # CA GRA
4 4 1 0.000000 0.000000 2.836000 0.000000 # CA GRA
5 5 1 0.000000 2.456048 0.000000 0.000000 # CA GRA
.
.
.
Bonds

1 1 1 2
2 1 1 6
3 1 2 3
4 1 3 4
5 1 4 7
.
.
.
Angles

1 1 2 1 6
2 1 1 2 3
3 1 2 3 4
4 1 3 4 7
5 1 3 4 85

and polymer.txt is:

# LAMMPS data file
20 atoms
19 bonds
18 angles

1 atom types
1 bond types
1 angle types


Atoms

1 1  1 -7.00000 25.00000 25 19.000000 
2 1  1 -7.00000 25.00000 25 18.035000 
3 1  1 -7.00000 25.00000 25 17.070000 
4 1  1 -7.00000 25.00000 25 16.105000 
5 1  1 -7.00000 25.00000 25 15.140000 
.
.
.
Bonds

1 1 1 2 
2 1 2 3 
3 1 3 4 
.
.
.
Angles

1 1 1 2 3 
2 1 2 3 4 
3 1 3 4 5 

What i’m getting after running is this:
when the delete_atom command is commented out

1008 atoms in group gcar
400 atoms in group gpot
400 atoms in group gchl
0 atoms in group gpoly

when the delete_atom command is not commented out, I’m getting this

948 atoms in group gcar
396 atoms in group gpot
393 atoms in group gchl
0 atoms in group gpoly

The group gpoly shouldn’t have 0 atoms. I’m not able to figure out why this is happening. Please help me out.

Thanks

If you write_data right after you read in both data files, I think you’ll find that your polymer chain is type 2, not type 4, because you only offset the atom types by 1 instead of 3. So I would expect that you would never have any atoms in gpoly. That doesn’t explain why you have 20 atoms when you don’t use the delete_atoms command, though.

I see that you’re not using the most recent version of LAMMPS, either, so I would suggest upgrading if changing toff to 3 in the read_data add command doesn’t work.

Also, is your polymer charge-neutral? Unrelated but important for the science.

I shifted the read_data polymer.txt after creating atom type 3 and set toff to 3. It worked. This too is documented well. I should have read it more carefully. When I ran the simulation by commenting out the delete_atom command, I set toff to 3. But when I ran the simulation after including the delete_atom, I mistakenly set toff to 1. I apologise for my naivete.

It is not as of now. I’ll add the same number of counter ions too in the simulation box. I was getting this anomaly in the initial stage itself so I didn’t add it right away.

1 Like