CONTACT_MATRIX does not account for per-atom weights (MFILTER_MORE output)

[russojohn]

Dear PLUMED developers,

I am trying to construct a CV for ice nucleation using PLUMED, and I am running into an issue with CONTACT_MATRIX that I do not fully understand.

I have successfully compiled GROMACS with PLUMED and verified that the basic multicolvars work correctly. In particular, I can compute the Steinhardt order parameter and identify “ice-like” molecules:

O: GROUP ATOMS=1-828:4
q6: Q6 SPECIES=O SWITCH={RATIONAL D_0=0.30 R_0=0.02 D_MAX=0.35}
lq6: LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=0.30 R_0=0.02 D_MAX=0.35}
ice: MFILTER_MORE DATA=lq6 SWITCH={GAUSSIAN D_0=2.00 R_0=0.01}

I have checked that ice behaves as expected and effectively returns a per-particle quantity that is 0 or 1 depending on whether the molecule is ice-like.

The issue arises when trying to build a contact matrix between ice-like molecules. What I would like is a matrix of the form:

M_{ij} = s(r_{ij}) , \text{ice}_i , \text{ice}_j

i.e., the usual distance-based contact matrix multiplied by the crystalline indicator of the two atoms.

However, when I do:

mat: CONTACT_MATRIX GROUP=ice SWITCH={CUBIC D_0=0.3 D_MAX=0.35}

I always obtain a matrix that depends only on distances, identical to what I get using the full oxygen group. It seems that CONTACT_MATRIX ignores the values of ice and only uses the underlying atom list.

From the documentation, it appears that CONTACT_MATRIX defines contacts purely based on distances, but I was wondering whether there is a way (or recommended approach) to incorporate per-atom weights such as ice_i * ice_j into the matrix elements.

So my questions are:

1. Is there a way in PLUMED to construct a contact matrix weighted by a per-atom quantity (e.g., ice_i * ice_j)?
2. If this is not directly supported, what is the recommended way to build a CV corresponding to the size of the largest crystalline cluster?

Any guidance or suggestions would be greatly appreciated.

Thanks in advance for your help!

John