Harmonic models and molecular dynamics simulations of isomorph behavior of Lennard-Jones fluids: Excess entropy and high temperature limiting behavior.

Henchman's approximate harmonic model of liquids is extended to predict the thermodynamic behavior along lines of constant excess entropy ("isomorphs") in the liquid and supercritical fluid regimes of the Lennard-Jones (LJ) potential phase diagram. Simple analytic expressions based on harmonic cell models of fluids are derived for the isomorph lines, one accurate version of which only requires as input parameters the average repulsive and attractive parts of the potential energy per particle at a single reference state point on the isomorph. The new harmonic cell routes for generating the isomorph lines are compared with those predicted by the literature molecular dynamics (MD) methods, the small step MD method giving typically the best agreement over a wide density and temperature range. Four routes to calculate the excess entropy in the MD simulations are compared, which includes employing Henchman's formulation, Widom's particle insertion method, thermodynamic integration, and parameterized LJ equations of state. The thermodynamic integration method proves to be the most computationally efficient. The excess entropy is resolved into contributions from the repulsive and attractive parts of the potential. The repulsive and attractive components of the potential energy, excess Helmholtz free energy, and excess entropy along a fluid isomorph are predicted to vary as ∼T-1/2 in the high temperature limit by an extension of classical inverse power potential perturbation theory statistical mechanics, trends that are confirmed by the MD simulations.

Departures from perfect isomorph behavior in Lennard-Jones fluids and solids.

Isomorphs are lines on a fluid or solid phase diagram along which the microstructure is invariant on affine density scaling of the molecular coordinates. Only inverse power (IP) and hard sphere potential systems are perfectly isomorphic. This work provides new theoretical tools and criteria to determine the extent of deviation from perfect isomorphicity for other pair potentials using the Lennard-Jones (LJ) system as a test case. A simple prescription for predicting isomorphs in the fluid range using the freezing line as a reference is shown to be quite accurate for the LJ system. The shear viscosity and self-diffusion coefficient scale well are calculated using this method, which enables comments on the physical significance of the correlations found previously in the literature to be made. The virial-potential energy fluctuation and the concept of an effective IPL system and exponent, n', are investigated, particularly with reference to the LJ freezing and melting lines. It is shown that the exponent, n', converges to the value 12 at a high temperature as ∼T-1/2, where T is the temperature. Analytic expressions are derived for the density, temperature, and radius derivatives of the radial distribution function along an isomorph that can be used in molecular simulation. The variance of the radial distribution function and radial fluctuation function are shown to be isomorph invariant.

Intrinsic viscuit probability distribution functions for transport coefficients of liquids and solids.

A reformulation of the Green-Kubo expressions for the transport coefficients of liquids in terms of a probability distribution function (PDF) of short trajectory contributions, which were named "viscuits," has been explored in a number of recent publications. The viscuit PDF, P, is asymmetric on the two sides of the distribution. It is shown here using equilibrium 3D and 2D molecular dynamics simulations that the viscuit PDF of a range of simple molecular single component and mixture liquid and solid systems can be expressed in terms of the same intrinsic PDF (P0), which is derived from P with the viscuit normalized by the standard deviation separately on each side of the distribution. P0 is symmetric between the two sides and can be represented for not very small viscuit values by the same gamma distribution formulated in terms of a single disposable parameter. P0 tends to an exponential in the large viscuit wings. Scattergrams of the viscuits and their associated single trajectory correlation functions are shown to distinguish effectively between liquids, solids, and glassy systems. The so-called viscuit square root method for obtaining the transport coefficients is shown to be a useful probe of small and statistically zero self-diffusion coefficients of molecules in the liquid and solid states, respectively. The results of this work suggest that the transport coefficients have a common underlying physical origin, reflecting at a coarse-grained level the traversal statistics of the system through its high-dimensioned potential energy landscape.

Strength of interference screw fixation of meniscus prosthesis matches native meniscus attachments.

PURPOSE: Meniscal surgery is one of the most common orthopaedic surgical interventions. Total meniscus replacements have been proposed as a solution for patients with irreparable meniscal injuries. Reliable fixation is crucial for the success and functionality of such implants. The aim of this study was to characterise an interference screw fixation system developed for a novel fibre-matrix-reinforced synthetic total meniscus replacement in an ovine cadaveric model. METHODS: Textile straps were tested in tension to failure (n = 15) and in cyclic tension (70-220 N) for 1000 cycles (n = 5). The textile strap-interference screw fixation system was tested in 4.5 mm-diameter single anterior and double posterior tunnels in North of England Mule ovine tibias aged > 2 years using titanium alloy (Ti6Al4Va) and polyether-ether-ketone (PEEK) screws (n ≥ 5). Straps were preconditioned, dynamically loaded for 1000 cycles in tension (70-220 N), the fixation slippage under cyclic loading was measured, and then pulled to failure. RESULTS: Strap stiffness was at least 12 times that recorded for human meniscal roots. Strap creep strain at the maximum load (220 N) was 0.005 following 1000 cycles. For all tunnels, pull-out failure resulted from textile strap slippage or bone fracture rather than strap rupture, which demonstrated that the textile strap was comparatively stronger than the interference screw fixation system. Pull-out load (anterior 544 ± 119 N; posterior 889 ± 157 N) was comparable to human meniscal root strength. Fixation slippage was within the acceptable range for anterior cruciate ligament graft reconstruction (anterior 1.9 ± 0.7 mm; posterior 1.9 ± 0.5 mm). CONCLUSION: These findings show that the textile attachment-interference screw fixation system provides reliable fixation for a novel ovine meniscus implant, supporting progression to in vivo testing. This research provides a baseline for future development of novel human meniscus replacements, in relation to attachment design and fixation methods. The data suggest that surgical techniques familiar from ligament reconstruction may be used for the fixation of clinical meniscal prostheses.

Fluid-solid interaction in the rate-dependent failure of brain tissue and biomimicking gels.

Brain tissue is a heterogeneous material, constituted by a soft matrix filled with cerebrospinal fluid. The interactions between, and the complexity of each of these components are responsible for the non-linear rate-dependent behaviour that characterises what is one of the most complex tissue in nature. Here, we investigate the influence of the cutting rate on the fracture properties of brain, through wire cutting experiments. We also present a computational model for the rate-dependent behaviour of fracture propagation in soft materials, which comprises the effects of fluid interaction through a poro-hyperelastic formulation. The method is developed in the framework of finite strain continuum mechanics, implemented in a commercial finite element code, and applied to the case of an edge-crack remotely loaded by a controlled displacement. Experimental and numerical results both show a toughening effect with increasing rates, which is linked to the energy dissipated by the fluid-solid interactions in the region surrounding the crack tip.

Viscuit and the fluctuation theorem investigation of shear viscosity by molecular dynamics simulations: The information and the noise.

The shear viscosity, η, of model liquids and solids is investigated within the framework of the viscuit and Fluctuation Theorem (FT) probability distribution function (PDF) theories, following Heyes et al. [J. Chem. Phys. 152, 194504 (2020)] using equilibrium molecular dynamics (MD) simulations on Lennard-Jones and Weeks-Chandler-Andersen model systems. The viscosity can be obtained in equilibrium MD simulation from the first moment of the viscuit PDF, which is shown for finite simulation lengths to give a less noisy plateau region than the Green-Kubo method. Two other formulas for the shear viscosity in terms of the viscuit and PDF analysis are also derived. A separation of the time-dependent average negative and positive viscuits extrapolated from the noise dominated region to zero time provides another route to η. The third method involves the relative number of positive and negative viscuits and their PDF standard deviations on the two sides for an equilibrium system. For the FT and finite shear rates, accurate analytic expressions for the relative number of positive to negative block average shear stresses is derived assuming a shifted Gaussian PDF, which is shown to agree well with non-equilibrium molecular dynamics simulations. A similar treatment of the positive and negative block average contributions to the viscosity is also shown to match the simulation data very well.

Normal Load and Counter Body Size Influence the Initiation of Microstructural Discontinuities in Copper during Sliding.

Near the interface of two contacting metallic bodies in relative motion, the microstructure changes. This modified microstructure leads to changes in material properties and thereby influences the tribological behavior of the entire contact. Tribological properties such as the friction coefficient and wear rate are controlled by the microstructure, while the elementary mechanisms for microstructural changes are not sufficiently understood. In this paper, the influence of the normal load and the size of the counter body on the initiation of a tribologically induced microstructure in copper after a single sliding pass is revealed. A systematic variation in the normal load and sphere diameter resulted in maximum Hertzian contact pressures between 530 MPa and 1953 MPa. Scanning electron microscopy, focused ion beam, and transmission electron microscopy were used to probe the subsurface deformation. Irrespective of the normal load and the sphere diameter, a sharp line-like feature consisting of dislocations, the so-called dislocation trace line, was identified in the subsurface area at depths between 100 nm and 400 nm. For normal loads below 6.75 N, dislocation features are formed below this line. For higher normal loads, the microstructure evolution directly underneath the surface is mainly confined to the area between the sample surface and the dislocation trace line, which itself is located at increasing depth. Transmission Kikuchi diffraction and transmission electron microscopy demonstrate that the misorientation is predominantly concentrated at the dislocation trace line. The results disclose a material rotation around axes roughly parallel to the transverse direction. This study demonstrates the generality of the trace line phenomena over a wide range of loads and contact pressures and the complexity of subsurface processes under a sliding contact and provides the basis for modeling the early stages in the microstructure evolution.

Modelling the effects of age-related morphological and mechanical skin changes on the stimulation of tactile mechanoreceptors.

Our sense of fine touch deteriorates as we age, a phenomenon typically associated with neurological changes to the skin. However, geometric and material changes to the skin may also play an important role on tactile perception and have not been studied in detail. Here, a finite element model is utilised to assess the extent to which age-related structural changes to the skin influence the tactile stimuli experienced by the mechanoreceptors. A numerical, hyperelastic, four-layered skin model was developed to simulate sliding of the finger against a rigid surface. The strain, deviatoric stress and strain energy density were recorded at the sites of the Merkel and Meissner receptors, whilst parameters of the model were systematically varied to simulate age-related geometric and material skin changes. The simulations comprise changes in skin layer stiffness, flattening of the dermal-epidermal junction and thinning of the dermis. It was found that the stiffness of the skin layers has a substantial effect on the stimulus magnitudes recorded at mechanoreceptors. Additionally, reducing the thickness of the dermis has a substantial effect on the Merkel disc whilst the Meissner corpuscle is particularly affected by flattening of the dermal epidermal junction. In order to represent aged skin, a model comprising a combination of ageing manifestations revealed a decrease in stimulus magnitudes at both mechanoreceptor sites. The result from the combined model differed from the sum of effects of the individually tested ageing manifestations, indicating that the individual effects of ageing cannot be linearly superimposed. Each manifestation of ageing results in a decreased stimulation intensity at the Meissner Corpuscle site, suggesting that ageing reduces the proportion of stimuli meeting the receptor amplitude detection threshold. This model therefore offers an additional biomechanical explanation for tactile perceptive degradation amongst the elderly. Applications of the developed model are in the evaluation of cosmetics products aimed at mitigating the effects of ageing, e.g. through skin hydration and administration of antioxidants, as well as in the design of products with improved tactile sensation, e.g. through the optimisation of materials and surface textures.

Single trajectory transport coefficients and the energy landscape by molecular dynamics simulations.

The Green-Kubo (GK) method is widely used to calculate the transport coefficients of model liquids by Molecular Dynamics (MD) simulation. A reformulation of GK was proposed by Heyes et al. [J. Chem. Phys. 150, 174504 (2019)], which expressed the shear viscosity in terms of a probability distribution function (PDF) of "single trajectory (ST) viscosities," called "viscuits." This approach is extended here to the bulk viscosity, thermal conductivity, and diffusion coefficient. The PDFs of the four STs expressed in terms of their standard deviations (calculated separately for the positive and negative sides) are shown by MD to be statistically the same for the Lennard-Jones fluid. This PDF can be represented well by a sum of exponentials and is independent of system size and state point in the equilibrium fluid regime. The PDF is not well reproduced by a stochastic model. The PDF is statistically the same as that derived from the potential energy, u, and other thermodynamic quantities, indicating that the transport coefficients are determined quantitatively by and follow closely the time evolution of the underlying energy landscape. The PDFs of out-of-equilibrium supercooled high density states are quite different from those of the equilibrium states.

Transport coefficients of the Lennard-Jones fluid close to the freezing line.

Molecular dynamics simulations have been carried out along four Lennard-Jones (LJ) fluid isomorphs close to the freezing line, covering a temperature, T, in the range of 0.8-350 and a number density, ρ, in the range of 1.1-3.0 in LJ units. Analysis of the transport coefficients is via the Green-Kubo time correlation function method. The radial distribution function, percolation threshold connectivity distance, self-diffusion coefficient, and shear viscosity are shown to be invariant along an isomorph to a very good approximation when scaled with Rosenfeld's macroscopic units, although there are some small departures for T ≃ 1 and lower temperatures. The thermal conductivity is shown for the first time also to be isomorph invariant. In contrast, the Einstein and moment-based frequencies, and especially the bulk viscosity, ηb, show poor isomorphic collapse at low T but not surprisingly tend to an "inverse power" potential limiting value in the high T limit. In the case of the bulk viscosity, the significant departures from invariance arise from oscillations in the pressure autocorrelation function at intermediate times, which scale for inverse power potential systems but not for the LJ case, at least in part, as the pressure and bulk elastic moduli are not isomorph invariant.

Shear stress relaxation and diffusion in simple liquids by molecular dynamics simulations: Analytic expressions and paths to viscosity.

The results are reported of an equilibrium molecular dynamics simulation study of the shear viscosity, η, and self-diffusion coefficient, D, of the Lennard-Jones liquid using the Green-Kubo (GK) method. Semiempirical analytic expressions for both GK time correlation functions were fitted to the simulation data and used to derive analytic expressions for the time dependent diffusion coefficient and shear viscosity, and also the correlation function frequency transforms. In the case of the shear viscosity for a state point near the triple point, a sech function was found to fit the correlation function significantly better than a gaussian in the ballistic short time regime. A reformulation of the shear GK formula in terms of a time series of time integrals ("viscuits") and contributions to the viscosity from components based on the initial stress ("visclets") enable the GK expressions to be recast in terms of probability distributions which could be used in coarse grained stochastic models of nanoscale flow. The visclet treatment shows that stress relaxation is statistically independent of the initial stress for equilibrium and metastable liquids, suggesting that shear stress relaxation in liquids is diffusion controlled. By contrast, the velocity autocorrelation function is sensitive to the initial velocity. Weak oscillations and a plateau at intermediate times originate to a greater extent from the high velocity tail of the Maxwell-Boltzmann velocity distribution. Simple approximate analytic expressions for the mean square displacement and the self Van Hove correlation function are also derived.

Instabilities of High Speed Dislocations.

Despite numerous theoretical models and simulation results, a clear physical picture of dislocations traveling at velocities comparable to the speed of sound in the medium remains elusive. Using two complementary atomistic methods to model uniformly moving screw dislocations, lattice dynamics and molecular dynamics, the existence of mechanical instabilities in the system is shown. These instabilities are found at material-dependent velocities far below the speed of sound. We show that these are the onset of an atomistic kinematic generation mechanism, which ultimately results in an avalanche of further dislocations. This homogeneous nucleation mechanism, observed but never fully explained before, is relevant in moderate and high strain rate phenomena including adiabatic shear banding, dynamic fracture, and shock loading. In principle, these mechanical instabilities do not prevent supersonic motion of dislocations.

Incremental viscosity by non-equilibrium molecular dynamics and the Eyring model.

The viscoelastic behavior of sheared fluids is calculated by Non-Equilibrium Molecular Dynamics (NEMD) simulation, and complementary analytic solutions of a time-dependent extension of Eyring's model (EM) for shear thinning are derived. It is argued that an "incremental viscosity," η i , or IV which is the derivative of the steady state stress with respect to the shear rate is a better measure of the physical state of the system than the conventional definition of the shear rate dependent viscosity (i.e., the shear stress divided by the strain rate). The stress relaxation function, C i (t), associated with η i is consistent with Boltzmann's superposition principle and is computed by NEMD and the EM. The IV of the Eyring model is shown to be a special case of the Carreau formula for shear thinning. An analytic solution for the transient time correlation function for the EM is derived. An extension of the EM to allow for significant local shear stress fluctuations on a molecular level, represented by a gaussian distribution, is shown to have the same analytic form as the original EM but with the EM stress replaced by its time and spatial average. Even at high shear rates and on small scales, the probability distribution function is almost gaussian (apart from in the wings) with the peak shifted by the shear. The Eyring formula approximately satisfies the Fluctuation Theorem, which may in part explain its success in representing the shear thinning curves of a wide range of different types of chemical systems.

Slip of Alkanes Confined between Surfactant Monolayers Adsorbed on Solid Surfaces.

The slip and friction behavior of n-hexadecane, confined between organic friction modifier surfactant films adsorbed on hematite surfaces, has been studied using nonequilibrium molecular dynamics simulations. The influence of the surfactant type and coverage, as well as the applied shear rate and pressure, has been investigated. A measurable slip length is only observed for surfactant films with a high surface coverage, which provide smooth interfaces between well-defined surfactant and hexadecane layers. Slip commences above a critical shear rate, beyond which the slip length first increases with increasing shear rate and then asymptotes toward a constant value. The maximum slip length increases significantly with increasing pressure. Systems and conditions which show a larger slip length typically give a lower friction coefficient. Generally, the friction coefficient increases linearly with logarithmic shear rate; however, it shows a much stronger shear rate dependency at low pressure than at high pressure. Relating slip and friction, slip only occurs above a critical shear stress, after which the slip length first increases linearly with increasing shear stress and then asymptotes. This behavior is well-described using previously proposed slip models. This study provides a more detailed understanding of the slip of alkanes on surfactant monolayers. It also suggests that high coverage surfactant films can significantly reduce friction by promoting slip, even when the surfaces are well-separated by a lubricant.

Towards the Irving-Kirkwood limit of the mechanical stress tensor.

The probability density functions (PDFs) of the local measure of pressure as a function of the sampling volume are computed for a model Lennard-Jones (LJ) fluid using the Method of Planes (MOP) and Volume Averaging (VA) techniques. This builds on the study of Heyes, Dini, and Smith [J. Chem. Phys. 145, 104504 (2016)] which only considered the VA method for larger subvolumes. The focus here is typically on much smaller subvolumes than considered previously, which tend to the Irving-Kirkwood limit where the pressure tensor is defined at a point. The PDFs from the MOP and VA routes are compared for cubic subvolumes, V=ℓ3. Using very high grid-resolution and box-counting analysis, we also show that any measurement of pressure in a molecular system will fail to exactly capture the molecular configuration. This suggests that it is impossible to obtain the pressure in the Irving-Kirkwood limit using the commonly employed grid based averaging techniques. More importantly, below ℓ≈3 in LJ reduced units, the PDFs depart from Gaussian statistics, and for ℓ=1.0, a double peaked PDF is observed in the MOP but not VA pressure distributions. This departure from a Gaussian shape means that the average pressure is not the most representative or common value to arise. In addition to contributing to our understanding of local pressure formulas, this work shows a clear lower limit on the validity of simply taking the average value when coarse graining pressure from molecular (and colloidal) systems.

Molecular Dynamics Studies of Overbased Detergents on a Water Surface.

Molecular dynamics (MD) simulations are reported of model overbased detergent nanoparticles on a model water surface which mimic their behavior on a Langmuir trough or large water droplet in engine oil. The simulations predict that the structure of the nanoparticle on a water surface is different to when it is immersed in a bulk hydrophobic solvent. The surfactant tails are partly directed out of the water, while the carbonate core maximizes its extent of contact with the water. Umbrella sampling calculations of the potential of mean force between two particles showed that they are associated with varying degrees with a maximum binding free energy of ca. 10 kBT for the salicylate stabilized particle, ca. 8 kBT for a sulfurized alkyl phenate stabilized particle, and ca. 5 kBT for a sulfonate stabilized particle. The differences in the strength of attraction depend on the proximity of nearest approach and the energy penalty associated with the disruption of the hydration shell of water molecules around the calcium carbonate core when the two particles approach. This is greatest for the sulfonate particle, which partially loses the surfactant ions to the solution, and least for the salicylate, which forms the weakest water "cage". The particles are separated by a water hydration layer, even at the point of closest approach.

On the effect of confined fluid molecular structure on nonequilibrium phase behaviour and friction.

A detailed understanding of the behaviour of confined fluids is critical to a range of industrial applications, for example to control friction in engineering components. In this study, a combination of tribological experiments and confined nonequilibrium molecular dynamics simulations has been used to investigate the effect of base fluid molecular structure on nonequilibrium phase behaviour and friction. An extensive parameter study, including several lubricant and traction fluid molecules subjected to pressures (0.5-2.0 GPa) and strain rates (104-1010 s-1) typical of the elastohydrodynamic lubrication regime, reveals clear relationships between the friction and flow behaviour. Lubricants, which are flexible, broadly linear molecules, give low friction coefficients that increase with strain rate and pressure in both the experiments and the simulations. Conversely, traction fluids, which are based on inflexible cycloaliphatic groups, give high friction coefficients that only weakly depend on strain rate and pressure. The observed differences in friction behaviour can be rationalised through the stronger shear localisation which is observed for the traction fluids in the simulations. Higher pressures lead to more pronounced shear localisation, whilst increased strain rates lead to a widening of the sheared region. The methods utilised in this study have clarified the physical mechanisms of important confined fluid behaviour and show significant potential in both improving the prediction of elastohydrodynamic friction and developing new molecules to control it.

Non-equilibrium phase behavior and friction of confined molecular films under shear: A non-equilibrium molecular dynamics study.

The phase behavior of a confined liquid at high pressure and shear rate, such as is found in elastohydrodynamic lubrication, can influence the traction characteristics in machine operation. Generic aspects of this behavior are investigated here using Non-equilibrium Molecular Dynamics (NEMD) simulations of confined Lennard-Jones (LJ) films under load with a recently proposed wall-driven shearing method without wall atom tethering [C. Gattinoni et al., Phys. Rev. E 90, 043302 (2014)]. The focus is on thick films in which the nonequilibrium phases formed in the confined region impact on the traction properties. The nonequilibrium phase and tribological diagrams are mapped out in detail as a function of load, wall sliding speed, and atomic scale surface roughness, which is shown can have a significant effect. The transition between these phases is typically not sharp as the external conditions are varied. The magnitude of the friction coefficient depends strongly on the nonequilibrium phase adopted by the confined region of molecules, and in general does not follow the classical friction relations between macroscopic bodies, e.g., the frictional force can decrease with increasing load in the Plug-Slip (PS) region of the phase diagram owing to structural changes induced in the confined film. The friction coefficient can be extremely low (∼0.01) in the PS region as a result of incommensurate alignment between a (100) face-centered cubic wall plane and reconstructed (111) layers of the confined region near the wall. It is possible to exploit hysteresis to retain low friction PS states well into the central localization high wall speed region of the phase diagram. Stick-slip behavior due to periodic in-plane melting of layers in the confined region and subsequent annealing is observed at low wall speeds and moderate external loads. At intermediate wall speeds and pressure values (at least) the friction coefficient decreases with increasing well depth of the LJ potential between the wall atoms, but increases when the attractive part of the potential between wall atoms and confined molecules is made larger.

Equilibrium fluctuations of liquid state static properties in a subvolume by molecular dynamics.

System property fluctuations increasingly dominate a physical process as the sampling volume decreases. The purpose of this work is to explore how the fluctuation statistics of various thermodynamic properties depend on the sampling volume, using molecular dynamics (MD) simulations. First an examination of various expressions for calculating the bulk pressure of a bulk liquid is made, which includes a decomposition of the virial expression into two terms, one of which is the Method of Planes (MOP) applied to the faces of the cubic simulation cell. Then an analysis is made of the fluctuations of local density, temperature, pressure, and shear stress as a function of sampling volume (SV). Cubic and spherical shaped SVs were used within a spatially homogeneous LJ liquid at a state point along the melting curve. It is shown that the MD-generated probability distribution functions (PDFs) of all of these properties are to a good approximation Gaussian even for SV containing only a few molecules (∼10), with the variances being inversely proportional to the SV volume, Ω. For small subvolumes the shear stress PDF fits better to a Gaussian than the pressure PDF. A new stochastic sampling technique to implement the volume averaging definition of the pressure tensor is presented, which is employed for cubic, spherical, thin cubic, and spherical shell SV. This method is more efficient for less symmetric SV shapes.

A localized momentum constraint for non-equilibrium molecular dynamics simulations.

A method which controls momentum evolution in a sub-region within a molecular dynamics simulation is derived from Gauss's principle of least constraint. The technique for localization is founded on the equations by Irving and Kirkwood [J. Chem. Phys. 18, 817 (1950)] expressed in a weak form according to the control volume (CV) procedure derived by Smith et al. [Phys. Rev. E. 85, 056705 (2012)]. A term for the advection of molecules appears in the derived constraint and is shown to be essential in order to exactly control the time evolution of momentum in the subvolume. The numerical procedure converges the total momentum in the CV to the target value to within machine precision in an iterative manner. The localized momentum constraint can prescribe essentially arbitrary flow fields in non-equilibrium molecular dynamics simulations. The methodology also forms a rigorous mathematical framework for introducing coupling constraints at the boundary between continuum and discrete systems. This functionality is demonstrated with a boundary-driven flow test case.