Skip to content

Action: VORONOI_COORDINATION

Module colvar
Description Usage
Calculate a scalar reduction of smooth Voronoi coordination defects. used in 0 tutorialsused in 0 eggs
output value type
the selected reduction of the smooth occupancy defects scalar

Details and examples

Calculate a scalar reduction of smooth Voronoi coordination defects.

Reactive processes such as proton transfer are difficult to describe with a fixed molecular identity because the atom that carries the proton can change. This Action assigns every atom in ASSIGNED continuously to the atoms in CENTERS and reduces the resulting coordination defects to one scalar. The same assignment is shared by VORONOI_DISTANCE and VORONOI_POSITION. Together the three Actions describe the amount, separation, and location of coordination defects without hard-coding water, glycine, a catalyst, or an atom-list position.

The construction follows the descriptors introduced for acid-base equilibria \cite Grifoni2019AcidBase and condensed-phase tautomerism \cite Grifoni2020Tautomeric. Applications to solvated glycine, interfacial water ions, electric-field effects, and electrocatalytic nitrogen reduction are discussed in \cite Zhang2024Glycine, \cite Zhang2025Interfaces, \cite Zhang2025ElectricField, and \cite Zhang2026NRR.

An illustrated bilingual guide with downloadable teaching structures, application examples, performance data, and additional references is available on the Reactive Soft-Voronoi CV project page.

VORONOI_COORDINATION, VORONOI_DISTANCE, and VORONOI_POSITION are part of PLUMED's default colvar module and have no external library dependency. In a PLUMED installation that contains these Actions, use them directly in the input; no LOAD line or optional module is required. Bias examples using OPES_METAD additionally require the opes module.

Soft assignment and coordination defects

For a center \f and an assigned atom \f, let \f be their minimum-image distance. The assignment weight is

\f[ w_{ij}=\frac{\exp(-\kappa d_{ij})} {\sum_k\exp(-\kappa d_{kj})}. \f]

By default, \f is evaluated with the minimum-image convention. NOPBC instead uses the direct coordinate difference, so the input coordinates must already follow a consistent unwrapped image convention.

The denominator contains all CENTERS for the same assigned atom, so \f. KAPPA is positive and has inverse units of the current PLUMED length unit. Increasing KAPPA sharpens the assignment toward the nearest center; reducing it spreads an assigned atom over more centers. The implementation uses a shifted softmax, which improves numerical stability without changing the mathematical value.

For monitoring or structural diagnosis, a larger KAPPA can provide sharper state labels once the relevant basins are already known. For biased sampling, a smaller KAPPA usually gives smoother assignment and force changes when an assigned atom switches between nearby centers, although a value that is too small can blur distinct chemical basins. Scan KAPPA on representative reactant, transition, product, and host-switching frames before applying a bias. Values such as 5 for smoother sampling or 100 for sharper diagnosis are application examples, not transferable defaults.

The smooth occupancy and coordination defect of center \f are

\f[ n_i=\sum_j w_{ij}, \qquad q_i=n_i-\nu_i . \f]

REFERENCE supplies \f. A single value is broadcast to all CENTERS; otherwise provide exactly one value per center in the same order as CENTERS. For example, REFERENCE=2 with two water O centers is the vector (2,2), not a one-center calculation. Fractional entries are valid when the model deliberately shares one reference occupancy over symmetry-related centers; they should not be introduced by averaging chemically nonequivalent sites. The identities

\f[ \sum_i n_i=N_{\mathrm{assigned}}, \qquad \sum_i q_i=N_{\mathrm{assigned}}-\sum_i\nu_i \f]

provide useful checks on a new chemical mapping. The Action does not infer elements, molecules, water, or a special reactive site from atom order. CENTERS and ASSIGNED must be nonempty, internally unique, and disjoint.

VORONOI_COORDINATION returns

\f[ Q_p=\sum_{i\in S} a_i q_i^p , \f]

where SELECT defines \f and COEFFICIENTS supplies \f. SELECT defaults to all CENTERS and COEFFICIENTS defaults to one. SIGN restricts the sum to \f or \f when POSITIVE or NEGATIVE is selected.

  • POWER=1 preserves the signed defect. With sign filtering it is non-differentiable exactly at \f.
  • POWER=2 measures defect activity without cancellation. With sign filtering its value and first derivative are continuous at \f, but its second derivative has a cusp there.
  • COEFFICIENTS can distinguish chemically different selected centers or reproduce a published scalar. Coefficients do not change the assignment.

Analytical coordinate and box derivatives are provided. The derivative with respect to a defect is \f for POWER=1 and \f for POWER=2 inside the selected sign branch, and zero outside it. These derivatives are propagated through every assignment weight. A CV intended for biasing should avoid a POWER=1 sign boundary, a neighbor-list membership change, or another non-smooth surface discussed below.

The geometric defects are not formal electronic charges. Their physical meaning comes from the chosen atom sets and reference occupancies and should always be checked for the system of interest.

Translating chemistry into keywords

Build the input from chemistry rather than from atom-list positions:

  1. Put atoms that can receive an assigned atom in CENTERS. For proton transfer these are commonly O and N atoms.
  2. Put only the transferable atoms in ASSIGNED. Hydrogen atoms that cannot participate in the process need not be included.
  3. Give every center its neutral or intended occupancy in REFERENCE. Typical examples are 2 for a water O and a model-dependent value for a reactive O or N.
  4. Start with the exact full-pair calculation. Inspect representative neutral, product, transition, and multi-defect configurations before choosing SELECT, POWER, SIGN, or COEFFICIENTS.
  5. Add VORONOI_DISTANCE only when separation is needed, and VORONOI_POSITION only when a fixed spatial frame is physically meaningful.

REFERENCE follows CENTERS order. COEFFICIENTS follows SELECT order. If the centers are reordered, reorder the corresponding numeric vector as well. Atom-valued SELECT, GROUP1, and GROUP2 lists use absolute atom numbers and do not rely on a center being last, first, or one of a fixed number of species.

Exact and neighbor-list calculations

Without NLIST, every CENTER-ASSIGNED pair is included and the finite-system definition above is exact. NLIST truncates the candidate centers and renormalizes the weights over the retained candidates, so it is an approximation rather than an algebraically exact acceleration. Every ASSIGNED atom must retain at least one CENTER or the calculation stops with an error.

Omitting NLIST therefore has no hidden cutoff. When NLIST is enabled, NL_CUTOFF is an absolute CENTER-ASSIGNED distance in the active PLUMED length units, not a bond cutoff or a water-specific constant. If \f is the nearest-center distance and \f is a trial cutoff, the screening relation

\f[ \exp[-\kappa(R-d_{\min})] \le \epsilon \f]

can provide an initial estimate for making a single omitted score small. Smaller KAPPA generally requires a larger cutoff. This estimate does not replace convergence because omitted scores, normalization errors, and force derivatives can accumulate over many centers.

NL_CUTOFF and NL_STRIDE are required whenever NLIST is present. NL_SKIN is optional, defaults to zero, and can be used only when NL_STRIDE is greater than one.

Before using NLIST in production:

  1. evaluate representative configurations with the exact full-pair form;
  2. increase NL_CUTOFF until values and forces agree within the required tolerance;
  3. start with NL_STRIDE=1, then increase it only after checking consecutive MD steps. With NL_STRIDE greater than one, NL_SKIN can add a Verlet buffer. The list is rebuilt early if any requested atom moves by more than half the skin or if the periodic box changes, while the evaluated pairs are still filtered at the true NL_CUTOFF on every step.

For replica exchange, every exchange step must also be a neighbor-list update step. Choose the update schedule accordingly; normally NL_STRIDE should divide the exchange stride.

For a fixed box, the half-skin displacement check ensures that a pair cannot cross NL_CUTOFF before it was present in the buffered candidate list. A larger skin reduces rebuilds but retains more candidate pairs, so both the skin and stride should be benchmarked. On PLUMED 2.11 and later, sufficiently large periodic candidate lists use the built-in link-cell broad phase; smaller lists retain the threaded pair scan because it is faster there. Runtime plugins built against older PLUMED versions keep the compatible pair scan, but can still amortize its rebuild cost with NL_SKIN.

A cutoff copied from another system is not a convergence test. Pair-list changes can introduce small discontinuities because the retained weights are renormalized. Prefer exact mode for derivative validation and for small or moderate systems. Use NLIST only after a value-and-force convergence scan demonstrates a useful speed/accuracy tradeoff for the target system.

CPU parallelism and practical scaling

The exact calculation evaluates \f pairs per step. For water with every O in CENTERS and every H in ASSIGNED, this is \f. Exact mode therefore remains quadratic even when it is parallel: additional CPU workers reduce elapsed time but do not change the asymptotic cost.

The pair, normalization, occupancy, and derivative loops use OpenMP when it is available. Set the PLUMED and OpenMP thread counts to the CPU cores allocated to each molecular-dynamics rank, for example:

export PLUMED_NUM_THREADS=4
export OMP_NUM_THREADS=4
plumed driver --plumed plumed.dat --ixyz trajectory.xyz --box 3.0,3.0,3.0

Benchmark 1, 2, 4, and 8 threads on the same frames because small candidate lists can spend more time entering parallel regions than doing pair work, and large lists can become limited by memory bandwidth. Avoid CPU oversubscription. In a GPU molecular-dynamics run these Actions still execute on CPUs, so a GPU allocation alone does not accelerate them. MPI also partitions pairs unless SERIAL is present; SERIAL is a debugging mode that repeats the work on every rank, not a performance option.

For large systems, a converged NLIST is the only option here that reduces the number of retained assignment pairs. NL_CUTOFF controls truncation accuracy, whereas NL_STRIDE controls how often membership is rebuilt. Start with NL_STRIDE=1, converge NL_CUTOFF against exact values and derivatives, and only then test a larger stride against the maximum atomic displacement between updates. Report CV time per call as well as whole-simulation throughput so that GPU force-model time and output time are not mistaken for CV cost.

Worked example 1: water autoionization

For water, oxygen atoms can be used as CENTERS, hydrogen atoms as ASSIGNED, and the neutral reference occupancy is two. In a configuration containing one hydronium and one hydroxide, the corresponding defects approach +1 and -1.

Click on the labels of the actions for more information on what each action computes
tested on2.11
UNITSThis command sets the internal units for the code. More details LENGTHthe units of lengths=A
WaterO: GROUPDefine a group of atoms so that a particular list of atoms can be referenced with a single label in definitions of CVs or virtual atoms. More details ATOMSthe numerical indexes for the set of atoms in the group=1-4
WaterH: GROUPDefine a group of atoms so that a particular list of atoms can be referenced with a single label in definitions of CVs or virtual atoms. More details ATOMSthe numerical indexes for the set of atoms in the group=5-12

ionization: VORONOI_COORDINATIONCalculate a scalar reduction of smooth Voronoi coordination defects. This action has hidden defaults. More details CENTERSAtoms that receive the smooth assignment=WaterO ASSIGNEDAtoms that are distributed over CENTERS=WaterH KAPPAPositive soft-assignment sharpness in inverse PLUMED length units=5 REFERENCEOne intended occupancy broadcast to all CENTERS, or one value per CENTER in CENTERS order=2 POWER Power of the selected occupancy defects; supported values are 1 and 2=2
positive_amount: VORONOI_COORDINATIONCalculate a scalar reduction of smooth Voronoi coordination defects. More details CENTERSAtoms that receive the smooth assignment=WaterO ASSIGNEDAtoms that are distributed over CENTERS=WaterH KAPPAPositive soft-assignment sharpness in inverse PLUMED length units=5 REFERENCEOne intended occupancy broadcast to all CENTERS, or one value per CENTER in CENTERS order=2 POWER Power of the selected occupancy defects; supported values are 1 and 2=2 SIGN Use ALL, POSITIVE, or NEGATIVE defects=POSITIVE
negative_amount: VORONOI_COORDINATIONCalculate a scalar reduction of smooth Voronoi coordination defects. More details CENTERSAtoms that receive the smooth assignment=WaterO ASSIGNEDAtoms that are distributed over CENTERS=WaterH KAPPAPositive soft-assignment sharpness in inverse PLUMED length units=5 REFERENCEOne intended occupancy broadcast to all CENTERS, or one value per CENTER in CENTERS order=2 POWER Power of the selected occupancy defects; supported values are 1 and 2=2 SIGN Use ALL, POSITIVE, or NEGATIVE defects=NEGATIVE
positive_signed: VORONOI_COORDINATIONCalculate a scalar reduction of smooth Voronoi coordination defects. More details CENTERSAtoms that receive the smooth assignment=WaterO ASSIGNEDAtoms that are distributed over CENTERS=WaterH KAPPAPositive soft-assignment sharpness in inverse PLUMED length units=5 REFERENCEOne intended occupancy broadcast to all CENTERS, or one value per CENTER in CENTERS order=2 POWER Power of the selected occupancy defects; supported values are 1 and 2=1 SIGN Use ALL, POSITIVE, or NEGATIVE defects=POSITIVE
negative_signed: VORONOI_COORDINATIONCalculate a scalar reduction of smooth Voronoi coordination defects. More details CENTERSAtoms that receive the smooth assignment=WaterO ASSIGNEDAtoms that are distributed over CENTERS=WaterH KAPPAPositive soft-assignment sharpness in inverse PLUMED length units=5 REFERENCEOne intended occupancy broadcast to all CENTERS, or one value per CENTER in CENTERS order=2 POWER Power of the selected occupancy defects; supported values are 1 and 2=1 SIGN Use ALL, POSITIVE, or NEGATIVE defects=NEGATIVE

PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=ionization,positive_amount,negative_amount,positive_signed,negative_signed FILEthe name of the file on which to output these quantities=COLVAR

For an isolated ion pair, ionization approaches 2, the two squared branches approach 1, and the signed branches approach +1 and -1. In a neutral frame all five values approach zero. Soft values between these limits are expected during proton transfer.

Worked example 2: exact-to-NLIST acceleration

The numeric cutoff, skin, and stride below are only input examples, not transferable recommendations. Use the three stages in order: establish an exact reference, converge NL_CUTOFF while rebuilding every step, and only then test whether a displacement-safe skin and a longer stride improve performance without changing values or forces.

Click on the labels of the actions for more information on what each action computes
tested on2.11
UNITSThis command sets the internal units for the code. More details LENGTHthe units of lengths=A
WaterO: GROUPDefine a group of atoms so that a particular list of atoms can be referenced with a single label in definitions of CVs or virtual atoms. More details ATOMSthe numerical indexes for the set of atoms in the group=1-4
WaterH: GROUPDefine a group of atoms so that a particular list of atoms can be referenced with a single label in definitions of CVs or virtual atoms. More details ATOMSthe numerical indexes for the set of atoms in the group=5-12

# 1. Exact full-pair reference
exact: VORONOI_COORDINATIONCalculate a scalar reduction of smooth Voronoi coordination defects. This action has hidden defaults. More details CENTERSAtoms that receive the smooth assignment=WaterO ASSIGNEDAtoms that are distributed over CENTERS=WaterH KAPPAPositive soft-assignment sharpness in inverse PLUMED length units=5 REFERENCEOne intended occupancy broadcast to all CENTERS, or one value per CENTER in CENTERS order=2 POWER Power of the selected occupancy defects; supported values are 1 and 2=2

# 2. Converge NL_CUTOFF with a rebuild every step
trial: VORONOI_COORDINATIONCalculate a scalar reduction of smooth Voronoi coordination defects. This action has hidden defaults. More details CENTERSAtoms that receive the smooth assignment=WaterO ASSIGNEDAtoms that are distributed over CENTERS=WaterH KAPPAPositive soft-assignment sharpness in inverse PLUMED length units=5 REFERENCEOne intended occupancy broadcast to all CENTERS, or one value per CENTER in CENTERS order=2 POWER Power of the selected occupancy defects; supported values are 1 and 2=2 NLIST Use an approximate neighbor-list truncation of the assignment candidates NL_CUTOFFCandidate cutoff in PLUMED length units; every ASSIGNED atom must retain at least one CENTER=8.0 NL_STRIDENumber of steps between neighbor-list updates=1

# 3. Only after stage 2, amortize rebuilds with a tested skin and stride
fast: VORONOI_COORDINATIONCalculate a scalar reduction of smooth Voronoi coordination defects. This action has hidden defaults. More details CENTERSAtoms that receive the smooth assignment=WaterO ASSIGNEDAtoms that are distributed over CENTERS=WaterH KAPPAPositive soft-assignment sharpness in inverse PLUMED length units=5 REFERENCEOne intended occupancy broadcast to all CENTERS, or one value per CENTER in CENTERS order=2 POWER Power of the selected occupancy defects; supported values are 1 and 2=2 NLIST Use an approximate neighbor-list truncation of the assignment candidates NL_CUTOFFCandidate cutoff in PLUMED length units; every ASSIGNED atom must retain at least one CENTER=8.0 NL_SKINNonnegative Verlet buffer in PLUMED length units; requires NL_STRIDE greater than one=1.0 NL_STRIDENumber of steps between neighbor-list updates=10

PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=exact,trial,fast FILEthe name of the file on which to output these quantities=COLVAR
DUMPDERIVATIVESDump the derivatives with respect to the input parameters for scalar values (generally CVs, functions or biases). More details ARGthe labels of the values whose derivatives should be output=exact,trial,fast FILEthe name of the file on which to output the derivatives=DERIVATIVES STRIDE the frequency with which the derivatives should be output=1

Worked example 3: applying a bias

First monitor the unbiased CV and check its scale. The following compact input then biases the total ionization activity with OPES. PACE, BARRIER, TEMP, and all production settings are system-dependent and must be justified for the simulation being run.

Click on the labels of the actions for more information on what each action computes
tested on2.11
UNITSThis command sets the internal units for the code. More details LENGTHthe units of lengths=A ENERGYthe units of energy=kj/mol
WaterO: GROUPDefine a group of atoms so that a particular list of atoms can be referenced with a single label in definitions of CVs or virtual atoms. More details ATOMSthe numerical indexes for the set of atoms in the group=1-4
WaterH: GROUPDefine a group of atoms so that a particular list of atoms can be referenced with a single label in definitions of CVs or virtual atoms. More details ATOMSthe numerical indexes for the set of atoms in the group=5-12
ionization: VORONOI_COORDINATIONCalculate a scalar reduction of smooth Voronoi coordination defects. This action has hidden defaults. More details CENTERSAtoms that receive the smooth assignment=WaterO ASSIGNEDAtoms that are distributed over CENTERS=WaterH KAPPAPositive soft-assignment sharpness in inverse PLUMED length units=5 REFERENCEOne intended occupancy broadcast to all CENTERS, or one value per CENTER in CENTERS order=2 POWER Power of the selected occupancy defects; supported values are 1 and 2=2
opes: OPES_METADOn-the-fly probability enhanced sampling with metadynamics-like target distribution. This action has hidden defaults. More details ARGthe labels of the scalars on which the bias will act=ionization PACEthe frequency for kernel deposition=500 BARRIERthe free energy barrier to be overcome=40 TEMP temperature=300
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=ionization,opes.bias FILEthe name of the file on which to output these quantities=COLVAR STRIDE the frequency with which the quantities of interest should be output=10

Passing a one-frame plumed driver test only establishes parsing and local evaluation. Before production biasing, verify analytical derivatives, energy/force units, restart behavior, the unbiased CV distribution, and a short molecular-dynamics force-path run.

Validation checklist

For every new chemical system:

  1. Check CENTERS, ASSIGNED, REFERENCE, SELECT, and COEFFICIENTS against a labeled structure rather than a topology-order assumption.
  2. Confirm the occupancy and defect conservation identities on neutral and reactive frames.
  3. Compare analytical and numerical derivatives away from sign cusps, zero-distance configurations, and periodic branch cuts.
  4. Translate the whole system by a lattice vector and verify the same value.
  5. Reorder CENTERS and ASSIGNED, reorder numeric vectors consistently, and verify the same value and forces.
  6. If NLIST is requested, converge both values and forces against exact mode.
  7. Compare serial and intended MPI/OpenMP execution on the same frames.
  8. Run a short, fixed-seed molecular-dynamics smoke test before enhanced sampling, then inspect finite COLVAR, bias, force, and restart output.

Troubleshooting and production cautions

  • Action VORONOI_COORDINATION is not known: use a PLUMED installation that contains these Actions and verify that the simulation is loading the same PLUMED kernel as plumed driver.
  • REFERENCE must contain ... values: provide one value for broadcast or one value per CENTER, in CENTERS order.
  • every atom in SELECT must also be present in CENTERS: SELECT is a subset, not an independent chemical group.
  • an ASSIGNED atom has no CENTER inside the candidate list: NL_CUTOFF is too small for that configuration or the atom sets are wrong. The calculation intentionally stops instead of returning a biased normalization.
  • Large changes after enabling NLIST indicate an unconverged truncation, not a harmless implementation detail.
  • Very large KAPPA makes the assignment nearly discrete and can produce sharp force changes near equidistant centers. Choose it from physical and numerical validation, not only from endpoint values.
  • POWER=1 sign filtering, ABSOLUTE positions, periodic position branch cuts, and neighbor-list membership changes are not globally smooth bias coordinates. Choose a smooth path or constrain the sampled region when these features are unavoidable.

Input

The atoms that serve as the input for this action are specified using one or more of the keywords in the following table.

Keyword Type Description
CENTERS atoms Atoms that receive the smooth assignment
ASSIGNED atoms Atoms that are distributed over CENTERS
SELECT atoms Subset of CENTERS included in the scalar reduction; the default is all centers

Full list of keywords

The following table describes the keywords and options that can be used with this action

Keyword Type Default Description
CENTERS input none Atoms that receive the smooth assignment
ASSIGNED input none Atoms that are distributed over CENTERS
SELECTThis keyword do not have examples input none Subset of CENTERS included in the scalar reduction; the default is all centers
KAPPA compulsory none Positive soft-assignment sharpness in inverse PLUMED length units
REFERENCE compulsory none One intended occupancy broadcast to all CENTERS, or one value per CENTER in CENTERS order
POWER compulsory 1 Power of the selected occupancy defects; supported values are 1 and 2
SIGN compulsory ALL Use ALL, POSITIVE, or NEGATIVE defects
NUMERICAL_DERIVATIVESThis keyword do not have examples optional false calculate the derivatives for these quantities numerically
NOPBCThis keyword do not have examples optional false ignore the periodic boundary conditions when calculating distances
SERIALThis keyword do not have examples optional false Perform the calculation redundantly on each rank for debugging
NLIST optional false Use an approximate neighbor-list truncation of the assignment candidates
NL_CUTOFF optional not used Candidate cutoff in PLUMED length units; every ASSIGNED atom must retain at least one CENTER
NL_STRIDE optional not used Number of steps between neighbor-list updates
NL_SKIN optional not used Nonnegative Verlet buffer in PLUMED length units; requires NL_STRIDE greater than one
COEFFICIENTSThis keyword do not have examples optional not used One coefficient, or one value per atom in SELECT