cna/atoms outputs incorrect structures for polymer simulation

Hello all,

I am having trouble using the cna/atoms compute command in a model polymer simulation. My simulation runs all the way through, but my dump file reads that cna finds all 5s. I know this result is incorrect because I can visualize the structures on ovito and run the same cna analysis there, and find a large fraction of particles have crystallized into a bbc structure.

In my input file, I run the command “compute cna all cna/atom 1.5” where 1.5 corresponds to the first minimum in the pair distribution function. I have extended the cutoff radius of my potential to 3.0 so that it is more than twice this radius and have modified my custom dump command to include “c_cna”. I wonder if this problem has something to do with the the atom_style being “bond”? Below is my input file.

Best,
Robert

#----------- Initial Setup -----------------------#

units lj

atom_style bond

boundary p p s # Film boundary conditions

I wonder if this problem has something to do with the the atom_style being “bond”?

yes, I think so. If you have bonds defined and special_bonds lj 0.0 for
the 1-2 bonds, that means 2 I,J atoms that are bonded do not appear
in the neighbor list. So compute cna/atom is not including them
in its analysis.

If you do a rerun command on the dump files of the bonded system
and do not define bonds in your rerun script, then a compute cna/atom
in the rerun script should do the right thing.

You could also change 0.0 to something like 1.0e-6 and see if
that works.

Axel - I can’t recall if there is another way to force compute cna/atom
(or similar computes) to use a neighbor list w/out the special bond
restrictions removed? I don’t think so …

Steve

>> I wonder if this problem has something to do with the the atom_style being "bond"?

yes, I think so. If you have bonds defined and special_bonds lj 0.0 for
the 1-2 bonds, that means 2 I,J atoms that are bonded do not appear
in the neighbor list. So compute cna/atom is not including them
in its analysis.

If you do a rerun command on the dump files of the bonded system
and do not define bonds in your rerun script, then a compute cna/atom
in the rerun script should do the right thing.

You could also change 0.0 to something like 1.0e-6 and see if
that works.

Axel - I can't recall if there is another way to force compute cna/atom
(or similar computes) to use a neighbor list w/out the special bond
restrictions removed? I don't think so ...

i agree. it might be helpful, if we put a test and warning message
into those computes to recommend replacing 0 with 1.0e-10 or so for
special_bonds when using any of those computes computing local
structure or clusters using a neighbor list.
we should file this as an issue, so we won't forget.

axel.