problem deleting bonds and atoms

LAMMPS LAMMPS Mailing List Mirror
1

Hi. I’m simulating a bead-sping polymer network and am trying to create a spherical void in the middle of the sample by deleting both bonds and atoms in a region ‘def’. The script

units lj
atom_style bond
read_data data.chain
special_bonds lj/coul 1 1 1
neighbor 0.4 bin
neigh_modify every 1 delay 1
communicate single vel yes
bond_style quartic
bond_coeff 1 1434.3 -0.7589 0.0 1.5 67.2234
dump mydump all custom 100000 chain_un*.xyz id type xu yu zu
timestep 0.01
fix 1 all npt temp 1.0 1.0 1.0 iso 0.0 0.0 1.0
thermo 100
thermo_style custom step temp press pxx pyy pzz lx ly lz
pair_style lj/cut 2.5
pair_coeff * * 1.0 1.0 2.5
pair_modify shift yes
region middle sphere 32.9266 32.9266 32.9266 19 side in units box
group def region middle
delete_bonds def atom 1 remove
delete_atoms group def
restart 10000 robsrestart
run 10000

produces the following crash:

Reading data file …
orthogonal box = (0.516376 0.516376 0.516376) to (65.8532 65.8532 65.8532)
1 by 1 by 1 MPI processor grid
248500 atoms
248500 velocities
248892 bonds
Finding 1-2 1-3 1-4 neighbors …
4 = max # of 1-2 neighbors
4 = max # of 1-3 neighbors
8 = max # of 1-4 neighbors
12 = max # of special neighbors
Finding 1-2 1-3 1-4 neighbors …
4 = max # of 1-2 neighbors
4 = max # of special neighbors
25598 atoms in group def
System init for delete_bonds …
Deleting bonds …
224422 total bonds, 224422 turned on, 0 turned off
Deleted 25598 atoms, new total = 222902
Setting up run …
ERROR on proc 0: Bond atom missing in box size check (domain.cpp:560)

I’m using the latest LAMMPS. Have tried various permutations of the above commands, to no avail. Any ideas as to how to delete the bonds and atoms without getting the ‘bond atom missing’ error?

Thanks,
Rob

2

Two thoughts. The output of delete bonds doesn't look
like any bonds were deleted, since none were turned off.

Second, this from the delete_bonds doc page:

For all styles, an interaction is only turned off (or on) if all the atoms involved are in the specified group. For style multi this is the >only criterion applied - all types of bonds, angles, dihedrals, impropers in the group turned off.

This will lead to the error you are getting, since a bond between
a def and non-def atom will not be deleted. When you delete
the def atoms, you thus have a dangling bond. So you want
to use the "any" keyword with delete_bonds. But you should
first figure out why no bonds are being flagged for removal.

Steve

3

Hi Rob
Sending this again. I forgot to CC the list, and I left out some details

  I wish I could say that my pet program "moltemplate" would
automatically solve your problem, but it won't, at least not
completely. (Moltemplate can deletes atoms by index location in an
array, not by spatial location.)

  You may have to write a short script to identify a list of atom ids
which fall in the region that you want to delete. Then go through the
Bonds, Angles, and Dihedrals sections, to throw away lines containing
references to any of these deleted atoms. (Python has associative
containers like "sets" which make this an easy script to write.)

   One advantage to this, is that way you can be sure that bonds,
angles, dihedrals are deleted whenever any of the atoms lie in the
deleted region. However, the "delete_bonds" command is not so
thorough. From the docs:

"For all styles, an interaction is only turned off (or on) if all the
atoms involved are in the specified group."
http://lammps.sandia.gov/doc/delete_bonds.html

I doubt this helps much, but I attached several scripts. One of them,
"extract_lammps_data.py", breaks a lammps data file into pieces:

extract_lammps_data.py Atoms < DATAFILE > Atoms.txt
extract_lammps_data.py Bonds < DATAFILE > Bonds.txt
extract_lammps_data.py Angles < DATAFILE > Angles.txt
extract_lammps_data.py Dihedrals < DATAFILE > Dihedrals.txt
extract_lammps_data.py header < DATAFILE > header.txt

(although I suppose you could just use a text editor)

To delete atoms, bonds, and angles, you could try setting their type
to 0 or -1. Otherwise if that does not work, you could just omit them
and then renumber the IDs of all the atoms, bonds, angles, and
dihedrals. One way to do this is using one of the other (poorly
named) scripts I attached.

renumber_DATA_first_column.py < BondsGaps.txt > BondsContiguous.txt
renumber_DATA_first_column.py < AnglesGaps.txt > AnglesContiguous.txt

Cheers

Andrew

extract_lammps_data.py (4.25 KB)

renumber_DATA_first_column.py (1.83 KB)

remove_duplicate_atoms.py (1.25 KB)

remove_duplicates_nbody.py (1.4 KB)

4

Hi, Steve. Actually it was definitely deleting bonds - the LAMMPS output shows there are 248892 bonds initially and 224492 after the delete_bonds command, i.e. 24400 bonds were deleted. And later stress-strain tests had confirmed that the bonds were deleted, i.e. the material gets softer.

That aside, I tried the modification you suggested, and further introduced a smaller spherical region for delete_atoms, to suppress the issue of having bonded atoms across the boundary. I changed the relevant part of the script to

sphere for deleting bonds

region middle sphere 32.9266 32.9266 32.9266 19 side in units box

slightly smaller sphere for deleting atoms

region middle2 sphere 32.9266 32.9266 32.9266 18 side in units box
group def region middle
group def2 region middle2

delete all bonds of type 1

delete_bonds def bond 1 any remove

delete atoms

delete_atoms group def2

And now it runs, but something is still awry. I write a restart file at step 1000, and when this occurs, I get a “Bond/angle/dihedral extent > half of periodic box length (domain.cpp:590)” error, i.e. the output is

LAMMPS (6 Dec 2012)
Scanning data file …
4 = max bonds/atom
Reading data file …
orthogonal box = (0.516376 0.516376 0.516376) to (65.8532 65.8532 65.8532)
1 by 1 by 1 MPI processor grid
248500 atoms
248500 velocities
248892 bonds
Finding 1-2 1-3 1-4 neighbors …
4 = max # of 1-2 neighbors
4 = max # of 1-3 neighbors
8 = max # of 1-4 neighbors
12 = max # of special neighbors
Finding 1-2 1-3 1-4 neighbors …
4 = max # of 1-2 neighbors
4 = max # of special neighbors
25598 atoms in group def
21779 atoms in group def2
System init for delete_bonds …
Deleting bonds …
224250 total bonds, 222257 turned on, 1993 turned off
Deleted 21779 atoms, new total = 226721
Setting up run …
Memory usage per processor = 145.344 Mbytes
Step Temp Press Pxx Pyy Pzz Lx Ly Lz
0 0.99970098 0.22851151 0.22108468 0.25434036 0.21010949 65.336839 65.336839 65.336839
1000 1.0035593 -0.027422399 -0.03216926 -0.024519467 -0.025578471 65.284217 65.284217 65.284217
ERROR: Bond/angle/dihedral extent > half of periodic box length (domain.cpp:590)

If I don’t delete the atoms, no such error occurs. So, it’s getting closer, but ???

Thanks,
Rob
Thanks,
Rob