There are many situations where you want to efficiently sample so-called "rare events"—such as conformational changes of molecules, bond formation/cleavage, or molecular desorption—that are difficult to reach simply by running conventional molecular dynamics (MD) simulations.
Enhanced sampling methods like umbrella sampling, metadynamics, and Steered MD (SMD) are useful for sampling these rare events. The most critical factor in these techniques is the selection of Collective Variables (CVs) or reaction coordinates. Your choice of CV changes the accuracy of the free energy landscape, sampling efficiency, and computational cost.
This article explains commonly used reaction coordinates alongside concrete simulation examples and how to write the corresponding PLUMED scripts. If you have not installed PLUMED yet, please also refer to an article entitled "How to install PLUMED"
Reaction Coordinates Using Single Atom Coordinates
Example 1: Pulling a Metal Atom from a Copper Slab Surface into a Solvent
To evaluate metal corrosion, catalyst degradation, or the elution of ions in solids, it is crucial to quantitatively determine the activation free energy when a metal atom desorbs from the slab surface into the solvent. Such a simulation can be executed by using the Z-axis coordinate of a single target atom as the reaction coordinate.
How to write the PLUMED script:
plumed_setting = [
f"UNITS LENGTH=A ENERGY=eV",
# Get coordinates of the target (PLUMED start atom index from :1)
f"pos: POSITION ATOM=1 NOPBC",
# Steered MD
"restraint: MOVINGRESTRAINT ARG=pos.z "
"STEP0=0 AT0=10.0 KAPPA0=2.5 "
"STEP1=200000 AT1=20.0 KAPPA1=2.5",
# Output
f"PRINT STRIDE=1000 ARG=pos.z,restraint.bias,restraint.force2 FILE=COLVAR",
f"FLUSH STRIDE=1000"
]- The
POSITIONcommand calculates the components of the position of the specified atom. The x, y, and z components can be extracted for the defined variable (pos) aspos.x,pos.y, andpos.z.
Important Considerations in Simulations
Note 1: Fix the lower layers of the slab
When pulling an atom from the surface, it is strongly recommended to fix the 1st and 2nd layers of the slab using FixAtoms or a similar method. If they are not fixed, the entire system will translate, and you will not be able to pull out just the specific atom.
Note 2: Periodic Boundary Conditions
When handling specific coordinates in PLUMED, you must be careful about the system explosion due to periodic boundary wrapping. PLUMED internally manages the periodic boundary box as [-L/2, L/2]. Therefore, coordinates that engines like ASE handle in the range of [0, L] will have their [L/2, L] portion wrapped to [-L/2, 0] inside PLUMED.
For example, suppose you pull an atom to around Z = 25.0 Å in a system with a Z-direction box size of 50.0 Å. The moment the atom slightly exceeds 25.0 Å, PLUMED internally determines that it has warped instantly to -25.0 Å. Because the spring potential for pulling is set around 25.0 Å, but the atom's coordinate is mapped to -25.0 Å, a massive, sudden force is applied, causing the system to explode.
To solve this, specify the NOPBC option when acquiring coordinates. This ensures that the unwrapped coordinates output by the MD engine are used exactly as they are, preventing discontinuous forces caused by unintended coordinate wrapping. (For details, refer to the PLUMED Manual).
Reaction Coordinates Using Center of Mass Coordinates
Pulling a single atom can cause a molecule to deform or decompose unnaturally. To avoid such unexpected deformations, it is appropriate to use the molecule's Center of Mass (COM) as the reaction coordinate.
Example 2: Pulling a Benzene Molecule Adsorbed on a Pt Slab into a Solvent
By using the relative distance between the center of mass of the slab and the molecule as the reaction coordinate, you can gradually move the adsorbed molecule away from the slab into the solvent without breaking it.
How to write the PLUMED script:
plumed_setting = [f"UNITS LENGTH=A ENERGY=eV",
# Specify the center of mass positions for the slab and the adsorbed molecule (Z-coordinate of COM)
"gslab: GROUP ATOMS=1-240",
"gmol: GROUP ATOMS=241-252",
"com1: COM ATOMS=gslab NOPBC", # Slab
"com2: COM ATOMS=gmol NOPBC", # Mol
"pcom1: POSITION ATOM=com1 NOPBC", # Slab
"pcom2: POSITION ATOM=com2 NOPBC", # Mol
# Calculate the difference in the Z-coordinates of the centers of mass
"distz: MATHEVAL ARG=pcom1.z,pcom2.z FUNC=y-x PERIODIC=NO",
# Steered MD settings
"restraint: MOVINGRESTRAINT ARG=distz "
"STEP0=0 AT0=10.0 KAPPA0=2.5 "
"STEP1=200000 AT1=20.0 KAPPA1=2.5",
# Output settings
f"PRINT STRIDE=1000 ARG=distz,restraint.bias,restraint.force2 FILE=COLVAR",
"FLUSH STRIDE=1000"
]The
GROUPcommand defines groups of atoms, allowing you to refer to a specific list of atoms with a single label.The
COMcommand calculates the center of mass of the specified atom group. The calculated center of mass is saved as a virtual atom.The
POSITIONcommand get the coordination information of the virtual atom created byCOM.The
MATHEVALcommand allows you to calculate CVs using arbitrary mathematical formulas. Here, the relative Z-coordinate is defined byy-x(=pcom2.z - pcom1.z).
When calculating COM, if the molecule spans across a periodic boundary, the center of mass position may become inaccurate (e.g., if atoms are split with one side at -24 Å and the other at +24 Å). In such cases, append NOPBC when calculating the center of mass.
Example 3: Movement of a Solute Molecule Between Water and Oil Phases
If you pull the molecule itself directly, the entire system might move due to interactions between the molecule and the surrounding solvent. In such cases, using the distance between the centers of mass of two domains as the CV can prevent the entire system from translating due to the bias potential. Additionally, in cases where the relative displacement of two molecular domains is the physically meaningful coordinate, using the COM distance between the two domains is an effective approach.
Below is an example of achieving a simulation that moves a solute molecule into an oil layer by using the distance between the center of mass of water and the center of mass of the solute molecule as the reaction coordinate.
How to write the PLUMED script:
plumed_setting = [
"UNITS LENGTH=A ENERGY=eV",
# Calculate the center of mass for solute and solvent
"mol: COM ATOMS=1-18 NOPBC",
"wat: COM ATOMS=19-3618 NOPBC",
# Calculate the distance between the centers of mass
"dist: DISTANCE ATOMS=mol,wat COMPONENTS",
# Steered MD settings
"restraint: MOVINGRESTRAINT ARG=dist.z "
"STEP0=0 AT0=5.4 KAPPA0=5.0 "
"STEP1=200000 AT1=30.0 KAPPA1=5.0",
# Output
"PRINT STRIDE=100 ARG=dist.x,dist.y,dist.z,restraint.bias,restraint.force2,restraint.work FILE=COLVAR_SMD",
"FLUSH STRIDE=1000"
]The
COMcommand calculates the center of mass of a specified atom group and labels it as a virtual atom. Here, the centers of mass for the solute molecule and the solvent water molecules are calculated.For
ATOMSin theDISTANCEcommand, you can directly specify not only actual atomic indices but also the virtual labels (mol,wat) defined byCOM.
Specific Interatomic Distances: Elongation and Dissociation
When the distance between two specific atoms directly characterizes a phenomenon—such as the dissociation of a covalent bond, the elongation of a polymer chain, or the collapse of a specific hydrogen bond network—using the distance between those two atoms as the reaction coordinate is optimal. This is also one of the simplest reaction coordinates.
Example 4: Unfolding of Decaalanine
Here, we look at an example of unfolding the α-helix structure of decaalanine (a chain of 10 alanines) by pulling the terminal C atoms.
How to write the PLUMED script:
plumed_setting = [f"UNITS LENGTH=A ENERGY=eV",
# Define the distance between atoms
"dist: DISTANCE ATOMS=9,99",
# Umbrella potential settings
f"restraint: RESTRAINT ARG=dist KAPPA=0.2 AT=15.0",
# Steered MD settings
"restraint: MOVINGRESTRAINT ARG=dist "
"STEP0=0 AT0=15.0 KAPPA0=0.2 "
"STEP1=200000 AT1=32.0 KAPPA1=0.2",
# Output settings
f"PRINT STRIDE=500 ARG=dist,restraint.bias,restraint.force2 FILE=COLVAR",
"FLUSH STRIDE=1000"
]- By specifying two atomic indices in
ATOMSwithin theDISTANCEcommand, you can use the distance between those atoms as the reaction coordinate.
Dihedral Angles: Conformational Exploration and Isomerization
Example 5: Conformational Exploration of the Backbone Dihedral Angles φ/ψ of Alanine Dipeptide
For analyzing the folding of protein backbones or rotational isomers of organic molecules, using dihedral angles rather than distances allows you to efficiently induce conformational changes without unnaturally distorting bond lengths or bond angles.
How to write the PLUMED script:
plumed_setting = [
"UNITS LENGTH=A ENERGY=eV",
# Definition of dihedral angles
"phi: TORSION ATOMS=5,7,9,15", # phi: (C_prev, N, CA, C)
"psi: TORSION ATOMS=7,9,15,17", # psi: (N, CA, C, N_next)
# Well-tempered Metadynamics
(
f"METAD "
f"ARG=phi,psi "
f"SIGMA=0.35,0.35 "
f"HEIGHT=0.01036 "
f"PACE=1000 "
f"BIASFACTOR=10.0 "
f"TEMP=300.0 "
f"GRID_MIN=-pi,-pi "
f"GRID_MAX=pi,pi "
f"GRID_BIN=360,360 "
f"LABEL=metad "
f"FILE={out_dir}/HILLS"
),
# Output
f"PRINT ARG=phi,psi,metad.bias STRIDE=1000 FILE={out_dir}/COLVAR",
f"FLUSH STRIDE=10000",
]- You can calculate torsion angles by specifying the indices of four atoms with the
TORSIONcommand.
Enzymatic Reactions: Bond Formation and Cleavage
In phenomena like enzyme-catalyzed reactions, simultaneous changes in two distances—the "forming bond" and the "breaking bond"—are often observed. If you use only a single distance as the reaction coordinate in this scenario, you will overlook the change in the other distance and fail to accurately represent the transition state. A widely used method here is to set the reaction coordinate as the difference between the two distances:
RC = dbreak - dnuc
Reactant State: Shorter and longer dnuc → RC < 0
Transition State: Both distances are roughly equal → RC ≒ 0
Product State: Longer and shoter dnuc → RC > 0
Therefore, as the reaction progresses, the reaction coordinate smoothly changes from negative to positive, allowing a single variable to effectively represent the entire reaction path, including the transition state.
How to write the PLUMED script:
plumed_setting = [f"UNITS LENGTH=A ENERGY=eV",
# Calculate the distances between the involved atoms
"dnuc: DISTANCE ATOMS=1,2", # Forming bond [A] - [B]
"dbreak: DISTANCE ATOMS=3,4", # Cleaving bond [C] - [D]
# Define the reaction coordinate
"rc: MATHEVAL ARG=dnuc,dbreak FUNC=y-x PERIODIC=NO",
# Steered MD settings
"restraint: MOVINGRESTRAINT ARG=rc "
"STEP0=0 AT0=-2.0 KAPPA0=2.2 "
"STEP1=200000 AT1=2.0 KAPPA1=2.2",
# Output settings
f"PRINT STRIDE=50 ARG=rc,restraint.bias,restraint.force2 FILE=COLVAR",
]
Conclusion
In structural sampling and reaction acceleration methods using PLUMED, the accuracy of the resulting free energy landscape and the sampling efficiency change depending on "which variables are selected as CVs."
Setting reaction coordinates can be difficult because it depends heavily on the system. However, it is best to start by trying the basic CVs introduced here (distance, center of mass, angle, differences) and observe the system's behavior (such as the presence of hysteresis or the influence of hidden degrees of freedom). If the desired conformational change does not occur, it is recommended to design more complex CVs (e.g., coordination number, contact map, RMSD, radius of gyration) or combine multiple CVs to expand into a two-dimensional space.