Custom Flow in Molecular Dynamics

Driving an inertial many-body system out of equilibrium generates complex dynamics due to memory effects and the intricate relationships between the external driving force, internal forces, and transport effects. Understanding the underlying physics is challenging and often requires carrying out case-by-case analysis. To systematically study the interplay between all types of forces that contribute to the dynamics, a method to generate prescribed flow patterns could be of great help. We develop a custom flow method to numerically construct the external force field required to obtain the desired time evolution of an inertial many-body system, as prescribed by its one-body current and density profiles. We validate the custom flow method in a Newtonian system of purely repulsive particles by creating a slow motion dynamics of an out-of-equilibrium process and by prescribing the full time evolution between two distinct equilibrium states. The method can also be used with thermostat algorithms to control the temperature.


I. INTRODUCTION
The precise application of a space-and time-resolved external force field can be used to drive a many-body system out-of-equilibrium in a controlled way. Analysing the response of a system to an external field is a primary method to calculate transport coefficients [1] such as shear [2,3] and bulk [4,5] viscosities, and the thermal conductivity [6]. Imposed pressure gradients, patterned substrates, capillary forces, electromagnetic fields, and centrifugal forces are examples of external fields that can be used in lab-on-a-chip devices [7] for the control of microflows [8]. However, the dynamics are complex due to far from trivial relationships between the external driving, the interparticle interactions, and transport effects. It is therefore difficult to predict the time evolution of a many-body system under the influence of an imposed external field, and case-by-case analyses are often required [9][10][11].
We consider here the inverse problem: to impose the desired dynamics and then find the corresponding external field. Such inversion is a valuable tool even at the level of individual particles. It allows for example the independent and simultaneous motion of several particles (that differ in either the shape [12] or the magnetic properties [13]) in arbitrary directions using a single external field. We focus here on many-body systems. We aim at specifying the dynamics of a many-body system, as given by the time evolution of the one-body density and the one-body current fields, and then find, with computer simulations, the corresponding space-and time-resolved external field.
From a fundamental view point, such inversion can be used to time-reverse the dynamics of a many-body system [14] and might offer new insights on irreversible processes [15,16]. The inversion can be also used as an * Matthias.Schmidt@uni-bayreuth.de † delasheras.daniel@gmail.com; www.danieldelasheras.com alternative to Gauss's principle of least constraint [17] in order to impose constraints on a dynamical system. From an applied view point, controlling the time evolution of the system instead of the external force acting on it can also be useful for the calculation of transport coefficients and relaxation times, especially in nanochannels where deviations from the Navier-Stokes formalism and from bulk behaviour are expected [18][19][20][21][22]. The study of memory effects [23,24], the design of lab-on-a-chip devices [25,26], and the determination of the slip length at the nanoscale [27] can also benefit from such inverse methodology. From a theoretical perspective, the existence in equilibrium of a unique mapping between the density distribution and the conservative external force forms the basis of quantum [28,29] and classical [30] density functional theory. In time-dependent quantum mechanical systems, the Runge-Gross theorem [31] ensures the existence of a unique mapping between the density distribution and a time-dependent external potential. A classical analogue of the Runge-Gross theorem was proposed by Chan and Finken [32]. The existence of a unique mapping between the kinematic fields and the external force field plays a central role in power functional theory, an exact variational principle for nonequilibrium classical many-body overdamped Brownian [33] and Hamiltonian systems [34] as well as for many-body quantum systems [35].
The rapid increase in computational power has made possible the development of numerical inverse methods that implement these unique mappings in equilibrium for both classical [36,37] and quantum [38,39] systems. It is hence possible to prescribe an equilibrium density distribution and find the corresponding external potential that generates the density using e.g. Monte Carlo simulations in the case of an equilibrium classical system. We have also developed a custom flow method for time-dependent overdamped Brownian systems [37]. The method is a valuable tool to generate specific flow and density patterns in a completely controlled way [14]. Custom flow is based on the exact one-body force balance equation that, in overdamped Brownian systems, relates the fric-tion (against the solvent) force field, the internal force field, the external force field, and the thermal diffusion. Using Brownian dynamics simulations, custom flow finds the external force required to generate the desired (imposed) time evolution of both the one-body density and the one-body current distributions. We prescribe elements that enter into the force balance equation, namely the density and the current distributions, and use an iterative scheme to find the generating external force field.
The force sampling method [40] uses a closely related idea: by sampling the one-body internal force field and using the force balance equation it is possible to obtain the one-body density distribution of an equilibrium system. The density distribution obtained using the forces acting on the particles [40][41][42][43][44] are more accurate than those obtained via the traditional counting of particles at space points.
Here, we present a custom flow method for classical many-body systems following Newtonian Dynamics. The method is motivated by the exact one-body force balance equation, Sec. II A, and it constructs iteratively the external force field that is required to generate the desired (target) time-evolution of both the density and the current distributions, Sec. II B. The method constitutes the solution of a complex inverse problem in statistical physics and implements numerically the map between the kinematic fields (density and current) and the external force field. Custom flow can be used with both conservative and non-conservative forces as well as with thermostats. We validate the method in a model system, Sec. III, of purely repulsive particles using several test cases, Sec. IV, including one with the Bussi-Donadio-Parrinello thermostat [45].

II. THEORY
A. One-body force balance equation Consider a classical system with N identical and mutually interacting particles following Newtonian dynamics. The equations of motion of the i−th particle are where m is the mass of the particle, r i denotes its position, p i = mv i is the momentum of the particle with v i its velocity, and f i is the total force acting on the particle, which in general consists of an imposed time-dependent external contribution, f ext (r i , t), and an internal contribution, −∇ i u(r N ). Here ∇ i is the partial derivative with respect to r i and u(r N ) is the interparticle potential energy with, r N = {r 1 . . . r N } the complete set of particle positions.
In molecular dynamics (MD) simulations the equations of motion, (1) and (2), are integrated in time. The observables of interest can be obtained as space-and timeresolved one-body fields. For example, the one-body density distribution is given by where δ(r) is the three dimensional Dirac delta distribution, the sum runs over all particles N , and r is the position vector. The brackets · denote a statistical average, which out-of-equilibrium is done at each time t over different realizations of the initial conditions (that is, the positions and the velocities of the particles at the initial time t 0 = 0). Differentiation of eq. (4) with respect to time yields the one-body continuity equatioṅ where the overdot indicates a time derivative and the one-body current is defined as Differentiating eq. (6) and using equations (1), (2), and (3) results in the exact one-body force balance equation See e.g. Ref. [34] for a more detailed derivation of eq. (7).
Here, the one-body internal force field is f int (r, t) = F int (r, t)/ρ(r, t), where F int is the internal force density field given by The last term in eq. (7) describes the transport effects that arise due to the one-body description. This transport term involves the one-body kinetic stress tensor where v i v i is a dyadic product such that τ is of second rank.

B. Custom flow in inertial systems
We present here an iterative method to construct the external force that generates a prescribed time evolution of the one-body fields ρ and J. The most general form of the iteration scheme reads +γ∇ ln ρ(r, t) ρ (k) (r, t) .
Here k denotes the iteration index and α, β, and γ are free non-negative prefactors that in general can carry spatial and temporal dependencies. The fields ρ and J are the prescribed (target) fields, andJ is also known since it follows directly from the prescribed J via partial time derivative.
The procedure to find f ext (r, t) iteratively using equation (10) is the following: (i) Run an MD simulation, evolving all the initial microstates from t 0 to t 0 + ∆t and sampling the space-resolved one body fields ρ (k) (r, t 0 + ∆t), J (k) (r, t 0 +∆t), andJ (k) (r, t 0 +∆t); (ii) use the sampled fields at iteration k to construct the external force for the next iteration k + 1 according to eq. (10), and (iii) iterate until the process converges and the external force at time t 0 +∆t is found. Convergence is achieved once the sampled fields ρ (k) , J (k) , andJ (k) are the same (within the desired numerical accuracy) as their target counterparts ρ, J, andJ. Next, advance the time from t 0 +∆t to t 0 + 2∆t and repeat the previous steps until the external force at t 0 + 2∆t is found. The process is repeated for the complete time evolution which is discretized in time steps ∆t. The idea behind Eq. (10) is simple but very useful: at each time, the external force at iteration k + 1 is that at the previous iteration k plus terms that (i) correct the deviations in the sampled fields with respect to the target fields and (ii) vanish if the target and sampled fields are identical. For example, if the current at a given position is higher (lower) than the desired one, the external force at that position decreases (increases) in the next iteration. Other correction terms are possible provided that they change the external force at each iteration in the right direction. For example, the third term on the right hand side of eq. (10) can be replaced by something like γ∇(ρ − ρ (k) ). The precise form of eq. (10) is motivated by the exact force balance equation, as we show in the Appendix A.
The non-negative prefactors α, β, and γ control how much the external force changes in one iteration. Using the exact one-body force balance equation (7) we obtain suitable expressions for them (see the Appendix A for a detailed calculation): Here k B is the Boltzmann constant and T 0 denotes the temperature of the initial state. Recall that the above expressions for α, β, and γ are not unique. The prefactors only fix the amount of change between two iterations. The method can in principle be also implemented by simply using non-negative constant prefactors. Also, not all three prefactors need to be present. Actually, having only α or only β is sufficient to find the external force iteratively. In cases for which the time dependent density distribution ρ determines the full dynamical evolution of the system it is also possible to work only with the coefficient γ. Such cases occur only if the current is free of both rotational and harmonic terms such that ρ alone fully determines J via eq. (5). In our particular implementation of Eq. (10) we iterate using only the target and sample currents. Hence, we set α to the value in Eq. (11) and set both β and γ to zero. Then, the iterative custom flow method we use here reads This iteration scheme is repeated at every ∆t. That is, we use eq. (14) to iteratively find f ext (r, t 0 + ∆t). We then advance time to t 0 + 2∆t and use eq. (14) to find f ext (r, t 0 + 2∆t). The process repeats until the complete time evolution is found.
The same algorithm can be used to find a suitable collection of initial microstates at t 0 such as e.g. microstates from an equilibrium system with a prescribed one-body density distribution. To this end we can start with an homogeneous equilibrium system and use custom flow to find microstates of another equilibrium system with the desired density distribution that serves as our initial state at t 0 . Alternatively, such initial set of microstates can be also found using the inversion between the external field and the density distribution for equilibrium systems described in Ref. [37].
At each time we initialize the iterative process (k = 0) using the external force which follows by making both f int and ∇ · τ zero everywhere in eq. (7). Another possible initialization is to include the contributions of the internal force and the kinetic stress tensor at the previous sampling time step ∆t on the right hand side of eq. (15). In principle the time step ∆t can be as small as the integration time step dt of the MD simulation. In practice, however, we use a larger time step (∆t/dt = 10) that does not compromise the accuracy of the calculation but still reduces the computational effort. In the MD integration algorithm of the equations of motion, we keep the external force constant between two consecutive time steps, t and t + ∆t. Interpolating the external force between two consecutive time steps can be problematic in cases where f ext changes drastically (such as e.g. if an external force is switched on at a specific time).
The iteration scheme (14) is particularly well suited since it is general and it requires to sample only J (k) in order to find the external force for the next iteration k+1. Sampling ofJ, that can be done by individually sampling the terms on the right hand side of eq. (7), is therefore not required. The proposed version of custom flow in eq. (14) does not require knowledge of any additional contribution or modification to the one-body force balance equation that might arise if e.g. a thermostat algorithm acts on the many-body level. As we demonstrate below our method works also with thermostats.

III. MODEL AND SIMULATION DETAILS
We implement the method in a system of N = 50 particles that interact via the purely repulsive Weeks-Chandler-Anderson interparticle-interaction potential [46]: Here, r is the distance between two particles, σ is the length scale, is the energy scale, and r c /σ = 2 1/6 is the cutoff radius which is set at the minimum of the 12-6 Lennard-Jones potential. We use τ = mσ 2 / as the time scale. The particles are located in a periodic threedimensional simulation box with lengths L x , L y , and L z . The origin of coordinates is located at the center of the simulation box. The system is inhomogeneous only in thê x-direction, which is discretized with bins of size 0.05 σ. The system remains homogeneous and without average flow in the other two directions,ŷ andẑ.
The many-body equations of motion, eqs. (1) and (2), are integrated in time using the velocity Verlet algorithm [47] with an integration time step dt/τ = 10 −4 . The particle positions are initialized randomly with the condition that the particles do not interact with each other (all interparticle distances are larger than r c ). The particle velocities are initialized according to a Maxwell-Boltzmann distribution. Hence, each velocity component is generated from a Gaussian distribution with zero mean and standard deviation 2k B T /m with T absolute temperature, which needs to be prescribed. The center of mass is set initially at rest. We then let the system equilibrate for 1 τ with no external force applied.
For the custom flow method we use a time step of ∆t/dt = 10 and average at each time over 2 · 10 6 trajectories (different realizations of the initial positions and velocities of the particles at t = 0) to obtain accurate results. Since the time step ∆t is small and the prefactor α has been carefully selected, only three iterations are required at each time for the method to converge to the desired external force. The one-body fields are sampled at every ∆t and used according to eq. (14) to find the external force for the next iteration.

IV. RESULTS
We illustrate the validity of the method with three examples. In Sec. IV A we measure the time evolution of the one-body fields ρ and J in a system subject to a spatially inhomogeneous external force which is switched on at t = 0 and then kept constant in time. We next construct the time-dependent external force required to slow down the observed dynamics by an arbitrary factor. In Sec. IV B, we incorporate a thermostat to demonstrate its easy implementation within custom flow. Finally, in Sec. IV C, we prescribe the full time evolution of the onebody density and the one-body current and then find the corresponding external force.

A. Slow motion dynamics
As a first example, we use custom flow to modify the time scale of a dynamical process. The particles are located in a box with dimensions L x /σ = 4, L y /σ = 8, and L z /σ = 10. We start with a homogeneous system at equilibrium at t = 0 and initial temperature k B T / = 0.5. Then, we switch on the following external potential with V 0 / = 1. Hence, the corresponding external force, which is constant in time for t > 0, acts only in thexdirection We let the system evolve for a total time of 10 τ , which is long enough to reach proximity to a new equilibrium state with an inhomogeneous density profile. The time evolution of the one-body density and current profiles are shown in Fig. 1a) and 1b), respectively. The external force field, shown in Fig. 1c), accelerates the particles toward the minima of the external potential, located at x/σ = ±1. Two density peaks grow from the homogeneous density distribution at t = 0, reaching the largest amplitude at t/τ ≈ 0.65. At short times, the shape of the one-body current resembles that of the external force (sine wave) and it increases in amplitude until it reaches its largest value at t/τ ≈ 0.3. Then, the amplitude of the current decreases until the current starts to flip sign, which occurs when the density peaks reach the largest amplitude (t/τ ≈ 0.65). Next, the density peaks decrease in amplitude and get broader since only some particles have enough momenta to overcome the external potential barrier. Most particles, however, cannot overcome the potential barrier. Instead, they climb partially the barrier (contributing to the broadening of the density peaks). Then, once the kinetic energies of the particles have been transformed into external energy, the particles start to move backwards towards the minima of the potential. This backward motion leads to an increase of the density peaks and the process repeats again in time. Now, however, the process is less intense since both the energy stored in the current and the external energy have been partially dissipated due to e.g. interparticle collisions and converted into thermal energy.
The described time evolution repeats in time creating a damped oscillatory behaviour. Eventually, the system reaches an equilibrium state at t/τ ≈ 9.5 with vanishing one-body current (within our numerical accuracy). Recall that this highly nontrivial time evolution of the one-body density and current profiles is produced by a simple external force, see eq. (18), that is switched on at t = 0 and then kept constant in time. A video showing the time evolution of the one-body fields is provided in the Supplemental Material [48]. Next, we use custom flow to find the external force required to reproduce this complex dynamics but in slow motion, i.e. slowed down by an arbitrarily chosen factor. By changing the time scale of the process we expect the external force of the slow motion system to be time dependent, which is indeed the case. We want to scale the time by a factor a. Hence, in the new system the density profile ρ a at time t is the same as the density profile in the original system at time at. That is, ρ a (r, t) = ρ(r, at). The time derivative in the continuity equation (5) implies that the current in the new system is also scaled by a factor a. That is, J a (r, t) = aJ(r, at). Therefore, scaling the time leads to a factor a in front of the current that needs to be considered when prescribing the target fields. Similarly, the time derivative of the current gets an additional scaling factorJ a (r, t) = a 2J (r, at).
In Figs. 1d) to 1f) we show the slow motion dynamics with scaling factor a = 0.5. Hence, the slow motion runs for 20 τ , i.e. twice the original total time. The sampled one-body density and current profiles coincide with their target. The time evolution of the one-body density, Fig. 1d), is the same as in the original system, Fig. 1a), but it proceeds only at half speed. Similarly the evolution of the one-body current, Fig. 1e) is two times slower and the amplitude is half of the original target current profiles, Fig. 1d). The time dependent external force that generates the slow motion (found with custom flow) is shown in Fig. 1f). At short times the shape of the external force resembles that in the original system, Fig. 1c), but its amplitude is reduced by a factor four. This was expected since at t = 0 the system is in equilibrium under no external force. That is, at t = 0 both ∇ · τ and F int are homogeneous and cancel each other in eq. (7). At short times, onlyJ a contributes to the external force in the force balance equation and hence the external force has the same shape as in the original system but it is rescaled by a factor a 2 sinceJ a = a 2J as discussed above.
Although the maximum amplitude of the density peaks occurs at 1.3 τ , the amplitude of the external force continues to grow and its shape deviates from a sinusoidal wave. The maxima and the minima of the external force are shifted towards the location of the density peaks at x/σ = ±1. While the amplitude of the density profile decreases, the extrema of the external force shift towards the minima of the one-body density at x = 0 and x/σ = ±2.
When the slow motion system reaches the equilibrium state, the shape of the external force resembles that in the original system but, interestingly, the amplitude is slightly smaller in the slow motion system. In slow motion less energy is dissipated due to the reduced value of the one-body current. Hence, also the temperature is slightly different. For the original time evolution the final temperature after equilibrium is reached is k B T / = 0.77 and in slow motion it is k B T / = 0.68 [49]. This temperature difference is responsible for the different amplitudes of the external force. Note that in equilibrium the transport term reduces to Hence, according to eq. (7) it is clear that two equilibrium systems with the same density distribution but at different temperatures are generated by external forces with different amplitudes.
Robustness of external force. Custom flow generates a noisy external force, see Fig. 1f). The iterative scheme, eq. (14), minimizes the error in the current profile since it is designed to converge if target and sampled currents coincide. As a result, the statistical fluctuations in J are directly translated into the external force field. Such fluctuations arise due to the finite number of realizations (recall that at each time we average over 2 · 10 6 trajectories that have evolved from different realizations of the initial particle positions and velocities). One could then think that the external force is tailored to the set of initial conditions and cannot be used in other circumstances. We demonstrate in the following that this is not the case and that the external force is indeed quite robust and independent of the small details. To this end we first smooth the external force in Fig. 1f) using a Fourier transform and eliminating the high frequency modes (only the lowest 15 modes are retained). The smoothed external force is depicted in Fig.  1i). Next, using the smoothed external force and also a different set of initial conditions (with the same number of microstates, i.e. 2 · 10 6 ) we sample the time evolution of the density and the current one-body fields, shown in Figs. 1g) and 1h), respectively.
Both the density and the current profiles obtained with the original external force and with the smoothed external force and a different set of initial conditions are similar and reproduce accurately the target profiles. Of course small differences occur, compare e.g. the current profiles at t/τ = 1 in Figs. 1e) and 1h). Hence, we conclude custom flow is suitable to reproduce the prescribed dynamics using other sets of initial conditions provided that enough trajectories are used. Furthermore, we want to stress that for the cases considered here the noise in the external force is not relevant to generate the target fields ρ and J. We expect that other filters that keep the structure and eliminate the noise can also be used to smooth the external force. Sampled (solid-lines) density (d) and current (e) profiles of a slow motion system with dynamics slowed down by a factor a = 1/2 (note that the indicated times are twice those in panels (a) and (b)). The dashed lines indicate the target profiles. (f) External force, color-coded as in (d), obtained with custom flow, that produces the slow motion dynamics. Sampled density (g) and current profiles (h) for a system under the influence of the external force field (i) which is a smoothed version of that shown in (f). The set of initial microstates used to obtain the averaged fields in the second and in the third column are not the same. A video showing the time evolution is provided in the Supplemental Material [48].

B. Thermostats
It is often the case that MD simulations are performed at constant temperature. In the following we show custom flow is also valid with algorithms to control the temperature (thermostats). Several types of thermostats can be implemented in MD [50]. In general, a thermostat acts on the many-body level by rescaling the particle velocities and modifying therefore the equations of motion and the integration algorithm. Hence, thermostats generate new terms in the one-body force balance equation (7). However, custom flow is designed such that sampling of these terms is not required to advance the iterative process. Custom flow uses only the external force at the previous iteration and the sampled current field, cf. eq. (14). The external force constructed with custom flow changes, of course, if a thermostat is used but the implementation of the method remains unchanged.
To illustrate the use of thermostats, we implement the well known Bussi-Donadio-Parrinello (BDP) thermostat [45], an extension of the Berendsen thermostat [51], that stochastically ensures a thermalized distribution of the kinetic energy.
Out of equilibrium, the kinetic energy has a contribution due to the net flow of the system which is unrelated to the temperature [50]. To control the temperature considering only the velocity fluctua-tions around the mean velocity [52], we can rescale the particle velocities using only the thermal kinetic energy: where it is important to note that v(r, t) is the spaceand time-dependent velocity profile (and not the center of mass velocity). The implementation of eq. (21) is particularly simple within custom flow since the velocity profile is known in advance.
To demonstrate that custom flow can be used with thermostats, we find the external force that generates the same time evolution of ρ and J as that in Figs. 1a) and 1b) but using the BDP thermostat. We set the time constant (required in the algorithm to control the temperature) to five times the integration time step dt of the simulation. We show in Fig. 2 results from both using the total kinetic energy and only the thermal energy, eq. (21), to rescale the velocities.
The sampled one-body density and one-body current are displayed in Figs. 2a) and 2b), respectively. The results from both thermostats are on top of each other and are indistinguishable from the case without a thermostat, Fig. 1. In Fig. 2c) we show the external force constructed with custom flow and smoothed using the same procedure as discussed above. We have verified that the smoothed external force also reproduces the target time evolution. The inclusion of the thermostat has a strong influence on the external force, which is now time dependent, in contrast to the original (constant energy) dynamics for which the external force is constant in time, eq. (18). Furthermore, the external force required to produce the desired dynamics depends on the thermostat. This is clearly visible at e.g. t/τ = 0.35 at which the current profile has its largest amplitude. At that time the flow kinetic energy has also its largest value and hence the two versions of the thermostat rescale the velocities differently. Differences are also noticeable in both the internal force field, see Fig. 2e), and especially in the transport term ∇ · τ , Fig. 2f), since τ is directly related to the total kinetic energy cf. eq. (9).
Custom flow can help to understand how different thermostat algorithms modify the physical properties of a system. In Fig. 2c) we plot the external force at t/τ = 0.35 calculated by using the sampled f int and ∇ · τ into the force balance equation (A1), which is exact only without thermostats. Using the thermal kinetic energy in the BDP thermostat results in an approximated external force, via eq. (A1), that is almost on top of the actual force generated with custom flow. The external force obtained with the original BDP thermostat (that uses the total kinetic energy) via the force balance equation (A1) shows a clear deviation from the force obtained in custom flow. We therefore conclude that using the total kinetic energy in the BDP thermostat induces a nontrivial contribution to the force balance equation that alters the flow. In contrast, using the thermal energy, only the velocity fluctuations are rescaled and the flow is left unchanged.
As expected, the time derivative of the currentJ, Fig. 2d), is very small for t = 0.35 τ since J reaches at that time the maximum amplitude. Also for t/τ ≥ 9.9 the system is very close to equilibrium andJ vanishes within the numerical accuracy.
Finally, we show in Figs. 2g), 2h), and 2i) the time evolution of the total kinetic energy, the thermal energy, and the internal potential energy, respectively. Shown are the original dynamics (no thermostat) and both the BDP thermostat and the modified version that uses only the thermal kinetic energy. The total kinetic energy is constant in time for the BDP thermostat as it should be by construction. In the modified version, the thermal kinetic energy, Fig. 2g), is constant in time but the total kinetic energy varies with time since the flow kinetic energy is kept unchanged. The original time evolution is clearly not at constant temperature since both the kinetic and the thermal kinetic energy vary substantially over time. The total internal potential energy is for neither thermostat constant in time, Fig. 2i). For a short period of time around t/τ = 0.35 there is a significant difference between both versions of the BDP thermostat due to the large amplitude of the one-body current.
Custom flow can be used as a new tool to analyse the quality and the physical consequences of the inclusion of thermostats in the dynamics of many-body systems. We have shown here that separating the flow and thermal kinetic energies, especially in systems with large magnitude of the one-body current, is advisable.

C. Tailoring inhomogeneous density profiles
In the previous examples we obtained ρ and J in a simulation for a fixed external force and used modified versions of them as target fields. In this last example, we show there is freedom to prescribe the fields provided that they represent a physical system. For example, the target ρ and J must obey the continuity equation.
We set a simulation box with dimensions L x /σ = 10, L y /σ = 5 and L z /σ = 5, and prescribe the one-body density with average density ρ 0 = N/(L x L y L z ) = 0.2 σ −3 , and constants Aσ 3 = 0.05 (maximum amplitude of the density inhomogeneity) and T 0 /τ = 0.5. At t = 0 the density is homogeneous, see eq. (22), and the system is in equilibrium. The density profile evolves according to eq. (22) for 0 < t < T 0 (two peaks grow from the initial homogeneous state). At t ≥ T 0 the one-body current is set to zero everywhere and therefore the inhomogeneous density profile remains stationary. The target current J follows from eq. (22) and the space integral of the continuity equation (5): with J x the x-component of J and C a constant that we set such that the total integral of the current vanishes dxJ x = 0. That is, for convenience we choose to not have motion of the center of mass. We calculate the target current analytically using eq. (23). Note that in our effective one-dimensional system with periodic boundary conditions the time evolution of ρ determines the current J up to a constant only. However, in higher dimensions the continuity equation alone is not enough to determine the current from the time-evolution of the density profile since any divergence-free field can be added to the current without altering ρ.
Using custom flow to construct the external force that generates the time evolution prescribed in eq. (22) yields the results shown in Fig. 3 (a movie is also included in the Supplemental Material [48]). Panels a) and b) of Fig. 3 show for different times the sampled one-body density and the one-body current, respectively. Both fields are in perfect agreement with their respective target fields, cf. eqs. (22) and (23). The largest amplitude in J occurs at t/τ = 0.25 which is also the time of the largest change in the density profile. The constructed external force, shown in Fig. 3c), is highly non-trivial and it is closely related to the behaviour ofJ, shown in panel Fig. 3d). Initially, the external force accelerates the particles towards the maxima of the one-body density. Then, around t/τ = 0.25 the external force flips sign and decelerates, therefore, the particles. At t/τ = 0.5 there is a jump in the time evolution of the external force due to the im- posed vanishingJ (compare the profiles at t/τ = 0.499 and 0.502). Custom flow finds the correct external force despite this drastic change in time. AfterJ vanishes, both ρ and J do not change anymore. Interestingly, the external force continues to evolve in time. This can be explained by memory effects occurring in both the internal force field f int , Fig. 3e), and the transport term ∇ · τ , Fig. 3f). Even thoughJ vanishes at t/τ > 0.5, the external force still needs to vary in time in order to cancel the time evolution of the internal and the transport terms. Custom flow is therefore a valuable tool to study memory effects [23,24,[53][54][55] since it allows to isolate memory contributions in the force balance equation from the time evolution of the density and the current fields.

V. DISCUSSION AND CONCLUSIONS
We have presented a numerical iterative scheme, see eqs. (10) and (14), to construct the external force required to achieve a given (prescribed) time evolution of a Newtonian many body-system, as specified by the onebody density and the one-body current. We have previously shown how in overdamped Brownian dynamics the exact one-body force balance equation can be directly used to construct a reliable custom flow method [37]. The external force is generated as the sum of different contributions that are sampled in the simulation. Here, we have followed a different and more general approach. We construct iteratively the external force by adding at each iteration terms that correct the external force in the right direction and that vanish when target and sample fields coincide. Although we have restricted ourselves here to inertial molecular dynamics, the method is general and can be also used in e.g. overdamped Brownian dynamics and Langevin dynamics. There the corresponding force balance equation can be used to recalculate suitable expressions for the prefactors α, β, and γ since they might be different.
The more general iterative scheme, eq. (10), modifies the external force based on three types of target-sampled differences that occur in the gradient of the density, in the current, and in the time derivative of the current. This is analogous to the three different fundamental equations in classical mechanics: the d'Alembert's principle based on particle displacements, the Jourdain's principle based on variations of particle velocities, and the Gibbs-Appell-Gauss's principle based on variations of particle accelerations. It is therefore not surprising that only one prefactor α, β, or γ in eq. (10) need to be present (using only γ is restricted to curl-free target current profiles as discussed in Sec. II B).
Since the calculation of the external force requires to sample only the current, but not the individual contributions to the force balance equation, it is straightforward to use custom flow together with a thermostat. Applying custom flow to systems with the same target fields but different thermostat algorithms provides new insight into the different mechanisms to control the tempera-ture since it is possible to precisely analyse how the individual contributions to the force balance equation are affected by different thermostats. The choice of thermostat can heavily influence the results [56,57] and custom flow might help to make an educated selection on how to control the temperature out-of-equilibrium, which is a delicate issue [58].
Custom flow can be used to prescribe target fields such that at least one contribution to the force balance equation vanishes. For example, the current vanishes after a certain time in the example of Sec. IV C. This can facilitate the study of memory effects and the structure of memory kernels, a topic of current interest [23,24,[53][54][55].
The external forces constructed with custom flow are in general noisy since they are tailored to the finite set of initial microstates used during the iterative process. Nevertheless, we have shown that a smoothed version of the external force, constructed by filtering out the high frequency terms, also produces the target dynamics within the numerical accuracy. We note however that we have stayed away from instabilities that might be the source of convergence issues.
Although they are only model situations, the examples considered here are demanding; the target fields vary substantially over distances comparable to the particle size and we have designed a case in Sec. IV C for which the resulting external force is discontinuous in time. Custom flow has in all cases found the external force field that produces the target fields within the numerical accuracy. Nevertheless, convergence issues can occur in e.g. strongly driven systems, flows with rapid spatial variations, and near the onset of mechanical and fluid instabilities. Integer arithmetic [16,59] and using small values for the prefactors α, β, and γ might help to mitigate some of the problems that might appear.
The existence of a unique mapping between the density distribution and a time-dependent external potential is at the core of time-dependent density functional theory [32]. Such mapping is not completely general but restricted to the occurrence of gradient-like forces only. Similarly, in the widely spread dynamical density functional theory [30,60], the internal force field is drastically approximated as a functional of the density distribution only. No functional dependence on the flow occurs. These limitations are solved in the formally exact power functional theory [33,34] that considers a functional dependence on all kinematic fields. Such dependence is required to properly describe e.g. shear migration [61], phase coexistence of active particles [62], and laning formation in binary mixtures [63]. Power functional theory relies on a mapping between the external force and both the density and the current distributions. Such mapping is indispensable to e.g. describe systems in which the current field contains non-gradient contributions (i.e. rotational and harmonic contributions) since the continuity equation links only the divergence of the current and the time evolution of the density profile. It is therefore perfectly possible to construct families of systems that e.g. share the same time evolution of the density profile but have different current profiles [14] and are therefore generated by different external forces. Custom flow provides the numerical evidence of the existence of the unique mapping between the external force and the kinematic fields.
Custom flow has proven to be an excellent tool to develop approximated power functionals in overdamped Brownian systems [14] and we expect it to be also of great help to develop approximate power functionals in Newtonian systems. To study large scale systems it can be useful to extend custom flow to adaptive resolution techniques for multiscale molecular dynamics simulations [64][65][66][67].

ACKNOWLEDGMENTS
We thank Sophie Hermann for a critical reading of the manuscript. This work is supported by the German Research Foundation (DFG) via project number 447925252. the first term of the right hand side of eq. (A3). Note that this is necessarily the case if the iterative process converges. Comparing eqs. (A4) and (10) yields β = m ρ . To find an expression for α we approximateJ(r, t) anḋ J (k) (r, t) in eq. (A4) bẏ J(r, t) = J(r, t) − J(r, t − ∆t) ∆t , (A5) where we have used that the sampled J (k) and the target J coincide at time t − ∆t. This is again necessarily the case if the process converges. Although a central time difference would be more precise, we use here the backward time difference since J (k) (r, t + ∆t) is unknown at time t. Inserting eqs. (A5) and (A6) into eq. (A4) and comparing the result to eq. (10) results in α = m ρ∆t .
Finally to find a suitable expression for γ we use the equilibrium expression for the transport term (19) to roughly approximate the second term in eq. (A4) by From which we obtain γ = k B T by comparison to eq. (10). We note that eq (A7) is a crude approximation but we are only interested in a suitable expression for the prefactor γ. The precise values of the prefactors is not critical for the method to converge (provided they are small enough such that the iterative process is stable). However, having suitable expressions is relevant to achieve a fast convergence since the prefactors control the amount of change of the external force from iteration to iteration.