Feature request: Ten Wolde criterion for identification of crystal clusters

Dear Steve and collaborators,

I was very excited when lammps introduced the orientorder/atom compute a few months ago. That motivated me to add the ‘twconnections/atom’ compute, which would allow one to identify a solid-like particle in a system with a large thermal fluctuation.

The compute would work like this:

1- Compute a normalized Ybar_lm complex (2l+1)-dimensional vector for each particle i (please refer to the orientorder/atom entry on the manual)

2- Compute the scalar product between i and its nearest neighbor particles j

3- If the scalar product is larger than a threshold, j is connected to i. The number of connections on each particle i is returned as output.

So, the compute would have the cutoff radius, the order of the harmonics and the threshold as parameters.

Is there anyone working on this implementation? If not, I would be glad to contribute. I believe this should be fairly simple.

Best regards,

Luis Goncalves

P.S.: Ten Wolde’s work reference is

P. R. ten Wolde et al., J. Chem. Phys, 104, 9932 (1996)

Luis - sounds interesting. Aidan may wish to comment
as he wrote the compute orientorder/atom.

There is a compute coord/atom which uses
a simple distance criterion for counting the coord #/atom.

Also there is a compute cluster/atom which uses
a simple distance criteria for identifying clusters.

It should not be hard to add options to either or
both of those computes (or possibly just derive new computes from them)
that uses your criteria. I.e. Instead of using
the Ri for each atom, use the per-atom value from compute orientorder/atom,
and compare the I,J values to a threshold.

We’d be happy to add such options/computes to LAMMPS.



This sounds like a nice addition to LAMMPS. I think the best thing would be to construct a new compute based on cluster/atom, ideally using C++ inheritance. Call it say compute cluster/tenwolde/atom (class ComputeClusterTenWoldeAtom, compute_cluster_tenwolde_atom.h/cpp). You could then replace the distance criterion with ten Wolde’s criterion, with a threshold on the dot product. Note that by using cluster analysis, you don’t need to place a threshold on the number of connections, because, liquid atoms will not form large clusters if the dot product threshold is set high enough. That deviates from ten Wolde, and so maybe another name is appropriate, say compute cluster/orientorder/atom.

You can jsut cut and paste the low level compute orientorder/atom functions into you new class. There are ways to avoid duplicating those functions:

  1. Do multiple inheritance, which I have never attempted
  2. Make the low-level functions static and put in a boop.h file, see KSPACE/kissfft.h


Hi, Aidan!

Actually, it would be very important to get the ‘number of connections’. When a new system is simulated, we run a bulk liquid and a bulk crystal at the same temperature and obtain histograms of connections as a function of the threshold. Then, the threshold is set based on a clear, unambiguous identification of the phases via the histogram pattern.

So I think a new compute based on the coord/atom using the functions of orientorder/atom would be simpler and less hard-coded. A cluster analysis could be performed on a second stage (e.g., via group and cluster/atom).

I also never did multiple inheritance, but I will look it up. :wink:


Thank you for your help!


I see. If you can unambiguously classify atoms in that way and then do distance-based cluster analysis, then that simplifies the problem. There may be systems where that does not work as well as a direct cluster analysis based on the dot product.

Here is a better idea I think.

One compute can call another. See compute heat/flux for
how it takes 3 other computes as inputs and invokes
them internally.

So I suggest changing the syntax of compute coord/atom to:

compue ID group coord/atom cstyle args …


cstyle = cutoff Rc
cstyle = orient/order cID

cID is the ID of the compute orient/order command.

Then coord/atom can invoke that cID and use
the per-atom values it generates to tally the
coord #. Note that you will likely need the orient/atom
values for ghost atoms, so compute coord/atom can
invoke comm->forward_comm_compute() and provide pack/unpack
functions for storing the ghost atom values. See
how compute cluster/atom does this, for example.

A similar extension could be made to compute cluster/atom
to allow it to use a compute’s values (orient/atom)
as input for determining clusters.

Note that if you do this for clusters, you can
use the cluster IDs as chunk IDs in the compute chunk/atom
and fix ave/chunk command. This means you
could calculate properties of each cluster of atoms
connctted via their orient/order params. E.g. KE
or potEng of each cluster.

This extension would also allow future options
for coord/atom and cluster/atom to use other per-atom
params as criteria. E.g. centro/atom or CNA/atom.


Yes, that would be good too. coord/atom would have a special option to count ‘connected neighbors’ based on the scalar product. The only problem is that right now orientorder/atom returns just the Ql and for the ten Wolde criterion the whole vector Ybar_lm is used. So there would be just changes in two computes

1- orientorder/atom returns a normalized Ybar_lm vector in addition to the Ql scalar

2- coord/atom receives a new option ‘orientorder/atom cID threshold’ that counts connected atoms based on the scalar products between neighbors. ‘cutoff’ and ‘degree’ comes from the orientorder/atom compute.

That is even simpler. Did I get your idea right, Steve?


That’s a good idea, because it allows for more code reuse. The main things that are still missing are:

  1. A modification to compute orientorder/atom that would separately store the full vector of 2*l+1 components of each order parameter, instead of just summing them to get Q_l. This would not be too hard to add as an option.
  2. The dot product of the vectors for atom i and j. Since that is not a per-atom quantity, it can not reside in orientorder/atom. The easiest thing would be to do it in coord/atom, as an alternative to computing Rij.



I did not this message from you on the mail list. It must be held up somewhere, probably by sourceforge. Anyway what you say matches my own response, right down to the numbering 1 and 2. It sounds like you are well able to take this on yourself. I can give you some helpful advice, but I can’t provide much time.


That is even simpler. Did I get your idea right, Steve?

yes, I think so. As Aidan says, you can
add options to compute orient/atom to
have it calculate the various new per-atom values
that are needed. Compute coord/atom can
extract them, communicate them for ghost atoms,
and do the pairwise dot-products (some with ghost atoms)
in place of the Rij criterion, as Aidan said.

Ditto for compute cluster/atom and what
it uses as a clustering criterion.


OK, I’ll do it.

If I find any problems, I’ll contact Aidan for support.