Constrained molecular dynamics
CP2K can constrain a collective variable (CV) during molecular dynamics with
MOTION/CONSTRAINT/COLLECTIVE. The CV is defined in
FORCE_EVAL/SUBSYS/COLVAR, and COLVAR in the constraint
section selects it by input order. Use
INTERMOLECULAR for a CV whose atoms are
not all in the same molecular object.
A fixed-distance window between two atoms can be set up as follows:
&FORCE_EVAL
...
&SUBSYS
...
&COLVAR
&DISTANCE
ATOMS 1 2
&END DISTANCE
&END COLVAR
&END SUBSYS
&END FORCE_EVAL
&MOTION
&CONSTRAINT
&COLLECTIVE
COLVAR 1
INTERMOLECULAR TRUE
TARGET [angstrom] 2.0
&END COLLECTIVE
&LAGRANGE_MULTIPLIERS ON
FILENAME constraint_force
COMMON_ITERATION_LEVELS 1
&END LAGRANGE_MULTIPLIERS
&END CONSTRAINT
&MD
...
&END MD
&END MOTION
Lagrange-multiplier output
LAGRANGE_MULTIPLIERS writes two records during a velocity-Verlet step:
Shake Lagrangian Multiplierscontains the position-constraint multipliers. These are the values relevant to a configurational constraint force and constrained thermodynamic integration (TI).Rattle Lagrangian Multiplierscontains the velocity-constraint multipliers. They enforce the time derivative of the constraint and are not a second configurational-force sample.
The values are raw CP2K internal quantities; specifying TARGET in another unit does not convert
the printed multipliers. A distance constraint therefore produces a multiplier in hartree/bohr,
while an angular constraint uses the corresponding internal angular unit. The file contains the
multipliers of all active constraints, first intramolecular constraints and then intermolecular
constraints. The ordering within each group is collective-variable, 3-by-3, and 4-by-6 constraints.
Blue-moon ensemble correction
Constrained molecular dynamics generally produces biased statistical distributions. The blue-moon ensemble average is used to correct these biased outputs and retrieve the statistical properties corresponding to unconstrained molecular-dynamics conditions. The printed SHAKE multiplier is not, in general, a complete blue-moon estimator. For one constrained reaction coordinate \(\xi\), the free-energy gradient then follows
where \(k_\mathrm{B}\) is the Boltzmann constant, \(T\) is the temperature, \(\left\langle \cdots \right\rangle_\xi\) denotes the time average of reaction coordinate \(\xi(\mathbf r_1,\dots,\mathbf r_N)\) by MD simulation.
Here \(Z\) defines the scalar mass metric
With the CP2K convention that the constraint force is \(-\lambda\nabla\xi\), the free-energy gradient contains a \(Z^{-1/2}\) reweighting and, for a general coordinate, an additional metric-derivative term (\(G\)). With Cartesian gradient \(\mathbf g=\nabla\xi\), Hessian \(H=\nabla\nabla\xi\), and diagonal Cartesian mass matrix \(M\), this term is
Equivalently, define the local mass-weighted normal derivative \(D_\xi=(M^{-1}\mathbf g/Z)\mathbin{\cdot}\nabla\). Then \(D_\xi\xi=1\) and \(G=\tfrac12 D_\xi\ln Z\). This is not a derivative along the constrained MD trajectory: within a fixed window \(\xi\) is constant, while \(Z\) can vary with the other degrees of freedom. Consequently, dividing changes of \(\ln Z\) between saved frames by changes of \(\xi\) does not provide this derivative. It can instead be evaluated at each frame by differentiating the metric in the specified normal direction; a full Hessian is not mathematically required.
CP2K currently writes \(\lambda\) but does not evaluate or print the complete corrected blue-moon estimator. The required metric terms therefore have to be evaluated during postprocessing for the chosen reaction coordinate. The standalone blue-moon postprocessor evaluates both terms for one fixed collective constraint, including distances, angles, torsions, signed point-plane distances, point-to-bond-center distances and their linear combinations. It requires the matching XYZ trajectory, SHAKE output, actual atomic masses, fixed cell, temperature and exact CV definition. Its README describes supported conventions, data alignment and limitations; it is not a general multiple-constraint estimator. See Komeiji, Chem-Bio Informatics Journal 7, 12 (2007) for the general expression and explicit algorithms for two common coordinates.
For the distance between two atoms,
Here \(Z\) is constant and the metric-derivative term is zero. The reweighting cancels, so the free-energy gradient reduces to \(-\langle\lambda\rangle_\xi\) with the sign convention above.
For the three-atom distance difference
the metric is
and depends on the instantaneous angle. Here \(r_{ij}=|\mathbf r_i-\mathbf r_j|\), \(r_{kj}=|\mathbf r_k-\mathbf r_j|\), \(\boldsymbol\rho_{ij}=(\mathbf r_i-\mathbf r_j)/r_{ij}\) and similarly for \(\boldsymbol\rho_{kj}\). Applying the general Hessian expression above gives
Thus \(G\) vanishes for equal distances or collinear bonds, but not for a general configuration. Both
reweighting and the \(k_\mathrm{B}TG\) term must otherwise be retained. In particular, the \(G=0\)
simplification stated in Eq. (35) of Komeiji does not follow from the general Eq. (8) for unequal,
noncollinear bonds; the postprocessor tests the general expression against independent analytic and
finite-difference evaluations. Its tests also include a recorded CP2K distance-difference trajectory
and a separate canonical-ensemble integration against the derivative of a model partition function.
The former checks the data-processing path, while the latter checks the thermodynamic estimator;
neither is a long-time CP2K sampling-convergence benchmark. Reproduction details and limitations are
given in tools/blue_moon/README.md.
Simplifications for a particular coordinate must not be assumed for an arbitrary COMBINE_COLVAR, coordination number, or multiple simultaneous constraints. Additional constraints, including fixed atoms or rigid molecules, change the accessible degrees of freedom and require an appropriate estimator.
Fixed windows and moving constraints
For equilibrium constrained thermodynamic integration, run independently equilibrated trajectories
at a series of fixed TARGET values, calculate the corrected free-energy gradient in every window,
and integrate it over the reaction coordinate. Check the sampling length, correlation time, window
spacing, integration direction, and unit conversion.
TARGET_GROWTH instead changes TARGET
linearly by TARGET_GROWTH * TIMESTEP at every MD step, optionally stopping at
TARGET_LIMIT. This is a moving-constraint
or slow-growth protocol. A simple example for the performance of CV as the distance between two
atoms:
&FORCE_EVAL
...
&SUBSYS
...
&COLVAR
&DISTANCE
ATOMS 1 2
&END DISTANCE
&END COLVAR
&END SUBSYS
&END FORCE_EVAL
&MOTION
&CONSTRAINT
&COLLECTIVE
COLVAR 1
INTERMOLECULAR TRUE
TARGET [angstrom] 2.0
TARGET_GROWTH [angstrom*fs^-1] 0.0008
TARGET_LIMIT [angstrom] 3.0
&END COLLECTIVE
&LAGRANGE_MULTIPLIERS ON
FILENAME constraint_force
COMMON_ITERATION_LEVELS 1
&END LAGRANGE_MULTIPLIERS
&END CONSTRAINT
&MD
...
&END MD
&END MOTION
CP2K does not integrate the work or turn the resulting trajectory into an equilibrium free-energy profile automatically. A finite pulling rate can cause lag, dissipation, and direction-dependent hysteresis, so such a trajectory must be analysed with a method appropriate to the intended nonequilibrium protocol. In the limit of infinitesimal change of \(\xi\), the irreversible work (\(W\)) describes the energy change between the initial and final state. The resulting work from slow growth (typically the irreversible work) can be related to the free-energy change \((\Delta A)\) via Jarzynski’s identity