Molecular dynamics of open systems: construction of a mean-field particle reservoir

TL;DR

This paper introduces a novel hybrid simulation method for open molecular systems that combines atomistic regions with a mean-field particle reservoir, enabling efficient and accurate modeling of thermodynamic, structural, and dynamic properties.

Contribution

It presents a new approach to simulate open systems by interfacing atomistic regions with a mean-field reservoir of tracers, improving flexibility and efficiency.

Findings

Successfully tested on liquid water properties

Maintains thermodynamic and structural accuracy

Offers conceptual and numerical advantages

Abstract

The simulation of open molecular systems requires explicit or implicit reservoirs of energy and particles. Whereas full atomistic resolution is desired in the region of interest, there is some freedom in the implementation of the reservoirs. Here, we construct a combined, explicit reservoir by interfacing the atomistic region with regions of point-like, non-interacting particles (tracers) embedded in a thermodynamic mean field. The tracer molecules acquire atomistic resolution upon entering the atomistic region and equilibrate with this environment, while atomistic molecules become tracers governed by an effective mean-field potential after crossing the atomistic boundary. The approach is extensively tested on thermodynamic, structural, and dynamic properties of liquid water. Conceptual and numerical advantages of the procedure as well as new perspectives are highlighted and discussed.

Equations

$${\mathbf{F}}_{\alpha \beta }=w(X_{\alpha })w(X_{\beta }){\mathbf{F}}_{\alpha \beta }^{\mathrm{AT}}+[1-w(X_{\alpha })w(X_{\beta })]{\mathbf{F}}_{\alpha \beta }^{\mathrm{CG}}.$$

$$U_{\mathrm{tot}}=U_{\mathrm{tot}}^{\mathrm{AT}}+\sum _{k\in \Delta \cup \mathrm{TR}}\varphi (X_k),$$

Full text

Open systems, particle reservoir, grand canonical molecular dynamics, multiscale simulations

Introduction.

Molecular dynamics (MD) simulations have proven as a powerful means of investigation of molecular systems for half a century now Chandler (2017). The steady increase of computing resources has brought larger system sizes and longer time scales into the scope of direct simulations, which may soon be run in parallel to experiments. Diverse application fields know about situations where the region of interest exchanges particles with an environment, e.g., droplet evaporation Heinen, Vrabec, and Fischer (2016), condensation in porous hosts Höft and Horbach (2015); Peksa et al. (2015), nanoflows Ganti, Liu, and Frenkel (2017), and chemical reactions in living cells Schöneberg et al. (2017); Engblom, Lötstedt, and Meinecke (2018) and catalysts Roa et al. (2018); Matera et al. (2014). Simulations of such open systems face severe limitations as they often depend on the efficient insertion of molecules into dense environmentsDelgado-Buscalioni and Coveney (2003); Delgado-Buscalioni, Kremer, and Praprotnik (2008).In Molecular Dynamics, the insertion of particles into dense simulation boxes has been achieved in various ways: by extending the phase space with an additional variable that controls the insertion and removal of molecules Cagin and Pettit (1991a, b); Ji, Cagin, and Pettit (1992); Weerasinghe and Pettit (1994); Lynch and Pettit (1997), by an extended Hamiltonian with additional fractional particles Palmer and Lo (1994); Lo and Palmer (1995); Shroll and Smith (1999); Kuznetsova and Kvamme (1999); Eslami and Müller-Plathe (2007) or by hybrid techniques where the insertion and removal of particles in the MD simulations are performed with a Monte Carlo approach Boinepalli and Attard (2003); Chempath, Clark, and Snurr (2003); Lsal, Smith, and Kolafa (2005). Similar recent efforts in the direction of grand canonical MD and the construction of a particle reservoir can be found in Refs.25; 26; 27. Here, we present a solution in which, different from the abovementioned schemes, molecules are inserted dynamically by converting point particles from a reservoir into molecules, building on ideas of the Adaptive Resolution Simulation (AdResS) technique Delle Site and Praprotnik (2017); Praprotnik, Delle Site, and Kremer (2005, 2008); Wang, Schütte, and Delle Site (2012). The approach may also be seen as a computational magnifying glass that allows one to model and observe at high resolution only the region of acute interest, thereby focusing computational resources on the essential degrees of freedom Krekeler et al. (2018). We evaluate the concept’s usefulness exemplarily for the simulation of liquid water, based on a number of thermodynamic, structural, and dynamic properties. So far, the reservoir of point particles primarily fulfills the role of a thermodynamic bath that ensures physical consistency in the high-resolution region. As a perspective, the reservoir particles may then act as tracers coupled to hydrodynamic fields in a continuum model, which would enable hybrid multi-scale simulations Delgado-Buscalioni and Coveney (2003); Delgado-Buscalioni, Kremer, and Praprotnik (2008); Petsev, Leal, and Shell (2015); Alekseeva, Winkler, and Sutmann (2016); Hu, Korotkin, and Karabasov (2018). The Adaptive Resolution Simulation (AdResS) technique Praprotnik, Delle Site, and Kremer (2005, 2008); Wang, Schütte, and Delle Site (2012) has been developed in the past years to fulfill the role of such a magnifying glass in molecular dynamics. In AdResS, the space is divided into two large regions: a high resolution (atomistic) and a low resolution (coarse-grained) region. The interface of these regions is characterized by a transition region where interactions, via a space-dependent switching function, smoothly change from atomistic to coarse-grained and vice versa. The method has been applied with success to several challenging systems, including solvation of large molecules Lambeth et al. (2010); Sablić, Praprotnik, and Delgado-Buscalioni (2016); Zavadlav, Podgornik, and Praprotnik (2017); Fiorentini et al. (2017); Agarwal, Clementi, and Delle Site (2017); Shadrack Jabes, Klein, and Delle Site (2018) and ionic liquids Krekeler and Delle Site (2017); Shadrack Jabes and Delle Site (2018); Shadrack Jabes et al. (2018) to cite a few. Recent developments of the method have shown that changing resolution through a space-dependent function slows computational performance and is overall not convenient for the implementation of the algorithm or its transferability from one code to another. For this reason, the design of a simpler interface was proposed that is based on an abrupt change of resolution without a switching function and that was extensively tested and shown to be accurate and efficient Krekeler et al. (2018). In parallel, Kremer and coworkers, within the framework of AdResS utilizing a smooth switching function, proposed to remove the two-body coarse-grained interactions in the coarse-grained region, reducing it to a gas of non-interacting particles Kreis et al. (2015). In the current work, we merge both ideas, creating a system with an abrupt interface between a region of atomistic resolution and a thermodynamic reservoir of non-interacting particles.

Note in this context that AdResS simulations require the action of an additional one-body force at the interface, known as the thermodynamic force. Together with a thermostat, it ensures a thermodynamic equilibrium among the differently-resolved regions such that thermodynamic properties as well as structural and dynamic correlations of the atomistically resolved subdomain match those of a fully atomistic simulation of the entire system. Such a force was derived from first principles of statistical mechanics and thermodynamics Poblete et al. (2010); Fritsch et al. (2012); Wang et al. (2013), and, together with the thermostat, will be the key ingredient in the construction of the reservoir of non-interacting particles below. The combined role of the thermostat and the thermodynamic force is that of creating a mean field in which the non-interacting particles are embedded, so that density and temperature of the tracers fluctuate around their equilibrium values. This requirement will be shown, through a numerical experiment for the case of liquid water at room conditions, to be the essential condition for a proper, grand canonical-like exchange of particles with the atomistic region, as expected by the principles of statistical mechanics of open systems Agarwal et al. (2015); Delle Site (2016); Delle Site and Praprotnik (2017).

AdResS in a nutshell.

The original version of AdResS, illustrated in Fig. 1a, defines the coupling between different regions of resolution through a smooth interpolation of the force between pairs of molecules ($\alpha,\beta$):

[TABLE]

Here, $\boldsymbol{F}_{\alpha\beta}^{\mathrm{AT}}$ is the total force exerted by the atoms of molecule $\beta$ on molecule $\alpha$ due to atomistic interactions; $\boldsymbol{F}_{\alpha\beta}^{\mathrm{CG}}=-\nabla U^{\mathrm{CG}}(\boldsymbol{R}_{\alpha}-\boldsymbol{R}_{\beta})$ is the coarse-grained force between the centers of mass of both molecules, $\boldsymbol{R}_{\alpha}$ and $\boldsymbol{R}_{\beta}$; and the switching function $w(X)$ interpolates smoothly between 1 and 0 in the transition region $\Delta$ with $X=\boldsymbol{n}\cdot\boldsymbol{R}$ denoting the projection of the position $\boldsymbol{R}$ on the interface normal $\boldsymbol{n}$. The global thermodynamic equilibrium across the entire simulation box is assured by the thermostat and by the thermodynamic force $\boldsymbol{F}_{\mathrm{th}}(x)$. The latter is a position-dependent force, $\boldsymbol{F}_{\mathrm{th}}(x)=(M/\rho_{0})\nabla p(x)$, that counteracts the pressure gradient between regions of different resolution Fritsch et al. (2012), which would result from entropy differences due to the mismatch in the number of degrees of freedom. It is calculated iteratively during the equilibration run from the gradient of the mass density $\nabla\rho^{(k)}(x)$ in the $\Delta$ region at the $k$-th iteration step: \boldsymbol{F}_{th}^{(k+1)}(x)=\boldsymbol{F}_{th}^{(k)}(x)-c\bigl{(}M/\rho_{0}^{2}\kappa_{T}\bigr{)}\nabla\rho^{(k)}(x), where $M$ is the mass of a molecule, $\kappa_{T}$ the isothermal compressibility, $c$ a constant to control the convergence rate, and $\rho_{0}$ the equilibrium density Poblete et al. (2010); Fritsch et al. (2012); Wang et al. (2013); Delle Site and Praprotnik (2017). The constant $c$ is tuned by starting from a small value (about 0.0015, i.e. slow convergence) and after a certain amount of iterations, stepwise increased. When variations become larger than the converge rate of the previous step then it is stepwise reduced until final convergence is reached with the criterion: $\frac{max|\rho^{(k)}(x)-\rho_{0}|}{\rho_{0}}\leq 0.02$. After $\boldsymbol{F}_{th}(x)$ has been determined, it remains fixed during the production runs.

In the latest variant of the method, proposed in Ref. 32, the simulation box is divided in only two regions of different resolution, atomistic ($\mathrm{AT}\cup\Delta$) and coarse-grained (CG), with the transition region $\Delta$ using the unchanged atomistic interactions throughout. The coupling is achieved in a direct way (Fig. 1b) by the pair potential $U_{\Delta}=\sum_{\alpha\in\Delta}\sum_{\beta\in\text{CG}}U^{CG}(\boldsymbol{R}_{\alpha}-\boldsymbol{R}_{\beta})$ acting only between pairs of particles at different sides of the $\Delta$/CG interface. The minimum size of the $\Delta$ region is fixed by the cutoff of the intermolecular (non-bonded) interaction. This implies that molecules of the AT region do not have any direct interaction with the molecules of the CG region and as a consequence the corresponding dynamics is determined only by atomistic interactions. The same set up can be used for large molecules, e.g. long polymers, in which case only the non-bonded interactions are subject to the adaptive resolution process according to the position of each polymer segment in the box, while intramolecular interactions remain atomistic in nature regardless of the segment’s position, see e.g. Refs.51; 52. Further, the thermodynamic force and the thermostat act only in the $\Delta$ and CG regions, i.e., the dynamics in the AT region follows the untweaked atomistic Hamiltonian. In $\Delta$, the thermodynamic force is applied together with a capping force to ensure that the sudden change in a molecule’s resolution does not generate unphysically large forces before a local equilibration is established (i.e., avoiding that a CG molecule crosses the CG/AT interface, instantaneously gains atomistic features, and is found in an undue overlap configuration with surrounding atomistic molecules). The evaluation of observables is restricted to the AT region, i.e., molecules in the $\Delta$ region are excluded as this region is contaminated by artifacts of the coupling. This simplified approach was tested with success in the case of both liquid water and ionic liquids Krekeler et al. (2018). The natural question that arises is whether the basic characteristics of molecular adaptivity, used in Ref. 32, can be further reduced, as suggested by the preliminary work of Kreis et al. (2015). In the following, we describe the merging of the ideas of Refs. 32 and 45 and discuss the horizons that such an approach opens.

From interacting coarse-grained particles to non-interacting tracers.

Figure 1c shows the set up for the AdResS simulation with the abrupt interface, a transition region $\Delta$ with atomistic resolution, and the coarse-grained region containing non-interacting particles that we refer to as tracers (TR). Thermodynamic equilibrium is assured by a thermostat, acting on all particles in the $\Delta$ and TR regions, and by the thermodynamic force $\boldsymbol{F}_{th}(x)$. This latter quantity is calculated in the iterative manner as explained above, however, now we lift its restriction to the $\Delta$ region and let it extend into the TR region. In fact we do not need anymore the condition that in the CG region the potential acting on the particles should be, for consistency, the coarse-grained potential only. Let us introduce the corresponding one-body potential: $\phi(x)=-\int_{x_{0}}^{x}\boldsymbol{F}_{th}(x^{\prime})\cdot\boldsymbol{e}_{x}\,\mathrm{d}x^{\prime}$ with $x_{0}$ the position of the AT/$\Delta$ interface. The total potential energy of the system is

[TABLE]

where $U_{\mathrm{tot}}^{\mathrm{AT}}$ is the total potential energy due to atomistic force fields in the AT region and the $k$-sum runs over all molecules in $\Delta$ and all tracers. To avoid numerical instabilities, a capping force is employed in $\Delta$ for molecules which, upon just having entered the atomistic region, experience forces that are orders of magnitude less probable than what is typical in the AT region;equivalent measures are found in Refs. 32 and 45. The capping is achieved by simply truncating each Cartesian component of the total force vector at a prescribed maximum value. The capping occurs only upon insertion of atomistic molecules, i.e., in the vicinity of the border to the TR region, and, despite being a technical artifact, its repercussions on the AT region are not noticeable. In summary, the action of the environment on each tracer depends only on the specific position of the tracer in space, independently from the overall tracers configuration, that is in $\mathrm{TR}\cup\Delta$ region each tracer experiences $\phi(x)$ and the thermostat as an effective mean field.

Test on Liquid Water.

We have tested the proposed scheme in simulations of liquid water at room conditions; technical details can be found at the end. The reference setup consisted of $2.7\text{\times}{10}^{4}$ SPC/E water molecules at ambient conditions within an elongated cuboid box of length $L_{x}=$60\text{,}\mathrm{n}\mathrm{m} and cross-section $\approx($3.7\text{\,}\mathrm{n}\mathrm{m}$)^{2}$. Van der Waals and electrostatic interactions were cut off at a distance of $r_{c}=$1.2\text{\,}\mathrm{n}\mathrm{m}. In the corresponding AdResS setup, the atomistic resolution is kept only in a small region of length $L_{\mathrm{AT}}=$3\text{,}\mathrm{n}\mathrm{m} (along the $x$-axis) expanded by transition regions $\Delta$ at both ends of length $L_{\Delta}=$1\text{\,}\mathrm{n}\mathrm{m}. The remaining volume TR ($L_{\mathrm{TR}}=$55\text{,}\mathrm{n}\mathrm{m}) is filled with an ideal gas of tracer particles, with each tracer representing one water molecule; the mass density of tracers is fixed to its atomistic value, here $\rho_{0}=${10}^{3}\text{\,}\mathrm{k}\mathrm{g}\,\mathrm{m}^{-3}. Employing the ideal gas law, we obtain a pressure in the TR region of $p_{\mathrm{TR}}=(\rho_{0}/M)k_{\mathrm{B}}T\approx$1.400\text{,}\mathrm{b}\mathrm{a}\mathrm{r}. Thus as with many coarse-grained models Fritsch *et al.* ([2012](#bib.bib47)); Kreis *et al.* ([2015](#bib.bib45)); Alekseeva, Winkler, and Sutmann ([2016](#bib.bib34)); Roy, Dietrich, and Höfling ([2016](#bib.bib53)), the pressure is orders of magnitude larger than in the corresponding atomistic model, here $p_{\mathrm{AT}}=$1\text{\,}\mathrm{b}\mathrm{a}\mathrm{r}. The corresponding isothermal compressibility $\kappa_{T}^{\mathrm{(TR)}}=1/p_{\mathrm{TR}}$ exceeds that of water in the AT region by merely a factor of 16. We have intentionally chosen a large reservoir of tracers and a small transition region, $L_{\Delta}\ll L_{\mathrm{TR}}$. The idea behind this choice is to test the computational stability of the algorithm. In fact, in a setup with a smaller TR region and larger AT and $\Delta$ regions, the extended domain with atomistic interactions is likely to more easily stabilize the system towards the wished target, which is fully atomistic. Instead, if for a large TR region and relatively small $\Delta$ and AT region, $\phi(x)$ and the thermostat can equilibrate such a large system and lead to accurate results, then one can be confident that the method is computationally stable. In fact, in the setup employed in Ref. 45, where a relatively small tracer region was coupled to a large atomistic region ($\mathrm{AT}\cup\Delta$) the system was stabilized essentially by the fully atomistic part. In order to remove any doubt about possible artifacts arising from the elongated (tube-like) shape of the simulation box, we have tested also a setup where the AT region is defined by a small sphere embedded in a large cube filled with tracers. The obtained structural properties display the same accuracy of the original tube-like setup.

In the following, we evaluate the new coupling scheme based on the mean density profile and various pair and auto-correlation functions. The mean density profile $\rho(x)$ obtained from AdResS production runs (Fig. 2a) is flat across the entire simulation box and reproduces the equilibrium value $\rho_{0}$ of the atomistic reference simulation within 1% (Fig. 2c), except very close to the $\Delta$/TR interface, where the relative deviation is 2%. These accuracy goals were chosen by us in the equilibration run as convergence criterion for the determination of the thermodynamic force.

The thermodynamic force $\boldsymbol{F}_{th}(x)$, used throughout in the production run, is plotted in Fig. 2b with the corresponding thermodynamic potential $\phi(x)$ shown in Fig. 2d. The latter exhibits a barrier next to the AT region and, towards the TR region, a larger trap of depth $1.1\,k_{\mathrm{B}}T$. In addition, we obtained some fine structure in $\phi(x)$, which is needed to achieve the small roughness of the density profile $\rho(x)$. In the AT region, $\boldsymbol{F}_{th}(x)=0$ and $\phi(x)=0$ by construction, and $\boldsymbol{F}_{th}(x)$ vanishes rapidly in the TR region, within less than the cutoff length $r_{c}$. Hence, the potential $\phi(x)$ rapidly approaches a constant $\phi_{\mathrm{TR}}$ as the distance to the $\Delta$ region increases. Tracer particles that enter the atomistic region gain an energy difference of about $\Delta\phi:=\phi(x_{\Gamma_{TR}})-\phi(x_{\Gamma_{AT}})$, here $x_{\Gamma_{TR}}$ means that the $x$-position of the particle is at the interface, $\Gamma_{TR}$, between the TR region and the $\Delta$ region, and, equivalently, $x_{\Gamma_{AT}}$ means that the $x$-position of the particle is at the interface, $\Gamma_{AT}$, of the AT region with the $\Delta$ region. Interestingly, $\Delta\phi$ compensates exactly the pressure difference between the AT and TR region, $\Delta\phi=(M/\rho_{0})(p_{\mathrm{TR}}-p_{\mathrm{AT}})$. Due to the high pressure of the TR region and making use of the ideal gas law, this reduces to $\Delta\phi\approx(M/\rho_{0})p_{\mathrm{TR}}=k_{\mathrm{B}}T$, that is $\Delta\phi\approx$2.5\text{,}\mathrm{k}\mathrm{J}\mathrm{/}\mathrm{m}\mathrm{o}\mathrm{l} at $T=$300\text{\,}\mathrm{K}. A measure for the orientation of water molecules is given by their electric dipole. The spatial profile of the dipole orientation matches that of a uniform liquid in the whole AT region, which is highly satisfactory (Fig. 3a); merely boundary effects are present towards the $\Delta/\mathrm{TR}$ interface due to the absence of molecular dipoles in the TR region, yet this is outside of the atomistic system of interest. We conclude the action of the thermodynamic force is concentrated in the transition region $\Delta$, while the thermostat in the TR and $\Delta$ regions is sufficient to keep the rest of the large tracer region in equilibrium.

We checked further that the passage of molecules through the transition region is not hindered and that there is a proper exchange between the AT and TR regions. Qualitatively, this is confirmed by coloring all molecules according to their initial region and monitoring their diffusion as time progresses, see supporting figures S1 and S2. As a more stringent test, we analyzed the fluctuations of the particle number in the AT region and compared their probability distribution $p(N)$ with that of the equivalent subregion in a fully atomistic reference simulation, see Fig. 3b. The peak values, $N_{\mathrm{AT}}\approx 1{,}400$, differ slightly by about 1.5% due to marginal differences in the density (Fig. 2c). Yet, after correcting for this, both distributions neatly collapse and approximately resemble a Gaussian distribution (Fig. 3c).

Detailed structural properties in the atomistic region have been verified in terms of the radial distribution functions $g_{\mathrm{OO}}(r)$, $g_{\mathrm{OH}}(r)$, and $g_{\mathrm{HH}}(r)$ of oxygen and hydrogen atoms. The results from the adaptive simulation are virtually indistinguishable from the fully atomistic reference, see Fig. 3d–f; in particular, all maxima and minima of the curves are very well reproduced. Eventually, we have tested the dynamics in the AT region by computing the molecular dipole-dipole autocorrelation functions, $C_{\mu\mu}(t)$, and the velocity-velocity autocorrelation function, $C_{vv}(t)$, as well as the hydrogen bond correlation function $C_{\mathrm{O}-\mathrm{H}}(t)$ (Fig. 3g–i) (for the specific definitions of the aforementioned correlation functions see Agarwal et al. (2015). The agreement with the results from the fully atomistic simulation is rather satisfactory.

Conclusions.

We have constructed a hybrid simulation scheme that couples a small open system of atomistic resolution with a thermodynamic bath of point particles (tracers) such that both subsystems are in thermodynamic equilibrium. The reservoir is characterized by its temperature and mass density, which can be viewed as a thermodynamic mean field the tracers are embedded into while displaying their canonical fluctuations. In the atomistic region, we have shown on the example of liquid water that the scheme reproduces thermodynamic, structural, and dynamic properties with high fidelity. Higher accuracy can be achieved with finer parametrizations of the thermodynamic force and a tighter convergence criterion for the resulting density profile. In practice, the huge TR region used here by purpose can be replaced by a much smaller volume as there are no spatial correlations between the tracers.

The new method provides a stark simplification over previous AdResS setups:

(i) it avoids the smooth force interpolation in the transition region,

(ii) there is no need for double resolution of molecules due to the absence of a coarse-grained potential,

(iii) it uses a computationally trivial reservoir.

Thus, our approach opens a perspective towards efficient non-equilibrium simulations, which would benefit from the easy insertion of reservoir particles. Further, the simple structure of the coupling suggests extensions that satisfy additional constraints, such as momentum conservation, with the ultimate goal of coupling the atomistic region to fluid-mechanical continuum fields.

Simulation details

The atomistic reference system consists of 27,648 SPC/E water molecules in a cuboid box of dimensions 60.00\text{\,}\mathrm{n}\mathrm{m}$\times($3.72\text{\,}\mathrm{n}\mathrm{m}$)^{2}$ at temperature $T=$300\text{\,}\mathrm{K}, mass density $\rho_{0}=${10}^{3}\text{,}\mathrm{k}\mathrm{g},\mathrm{m}^{-3}, and pressure $1\text{\,}\mathrm{b}\mathrm{a}\mathrm{r}$. Simulations were performed with GROMACS Abraham *et al.* ([2015](#bib.bib54)) 5.1.0 in the NVE ensemble with a time step of $2\text{\,}\mathrm{f}\mathrm{s}$. The van der Waals and electrostatic interactions were treated with the “switch” method and the reaction-field method, respectively Agarwal *et al.* ([2015](#bib.bib49)), and the cutoff radius for each was set to $r_{c}=$1.2\text{\,}\mathrm{n}\mathrm{m} as in our previous study Krekeler et al. (2018). After equilibration, observables were recorded from one production run over $10\text{\,}\mathrm{n}\mathrm{s}$.

The corresponding AdResS simulations were run with a modified and extended GROMACS code, with the following essential differences to the reference setup: the elongated box is divided along the long $x$-axis in one AT region ($L_{\mathrm{AT}}=$3.0\text{,}\mathrm{n}\mathrm{m}), two $\Delta$ regions with $L_{\Delta}=$1\text{\,}\mathrm{n}\mathrm{m}, and two large TR regions connected through periodic boundaries ($L_{\mathrm{TR}}=$55\text{,}\mathrm{n}\mathrm{m}$$ in total). Thus, on average only 2,300 water molecules are treated atomistically ($\mathrm{AT}\cup\Delta$ region). For tracers entering the $\Delta$ region, atomistic forces are capped at a threshold of $2.000\text{\,}\mathrm{k}\mathrm{J}\,\mathrm{m}\mathrm{o}\mathrm{l}^{-1}\,\mathrm{n}\mathrm{m}^{-1}$. The stochastic dynamics integrator (Langevin thermostat) was applied in the $\Delta\cup\mathrm{TR}$ region with a relaxation time constant of $0.05\text{\,}\mathrm{p}\mathrm{s}$ and a time step of $2\text{\,}\mathrm{f}\mathrm{s}$, while in the AT region the thermostat part was switched off. Here we have employed the simplest option of a single thermostat acting on $\Delta\cup\mathrm{TR}$, however the option of thermostating the two regions separately may be technically more efficient (see e.g. Ref.55) and it will be explored in future work. The thermodynamic force $\boldsymbol{F}_{th}(x)$ points along the $x$-axis by symmetry and calculated from cubic splines parametrized on a uniform grid of spacing $0.16\text{\,}\mathrm{n}\mathrm{m}$. It is obtained iteratively during the equilibration procedure, which is interrupted every $1\text{\,}\mathrm{n}\mathrm{s}$ to readjust $\boldsymbol{F}_{th}(x)$; after 30 iterations, the desired flatness of the density profile $\rho(x)$ of 1% was achieved. The AdResS production run was then performed with the converged thermodynamic force for a duration of $10\text{\,}\mathrm{n}\mathrm{s}$.

Acknowledgements.

We are grateful to Matej Praprotnik for useful discussions. This research has been funded by Deutsche Forschungsgemeinschaft (DFG) through the grant CRC 1114: “Scaling Cascades in Complex Systems”, project C01. This work has also received funding from the European Union’s Horizon 2020 research and innovation program under the grant agreement No. 676531 (project E-CAM).The simulations were performed with the HPC resources provided by the North-German Supercomputing Alliance (HLRN), grant no. bec00171

Supporting material