# Adaptive Resolution Molecular Dynamics Technique: Down to the Essential

###### Abstract

We investigate the role of the thermodynamic (TD) force, as an essential and sufficient technical ingredient for an efficient and accurate adaptive resolution algorithm. Such a force applied in the coupling region of an adaptive resolution Molecular Dynamics (MD) set-up, assures thermodynamic equilibrium between atomistically resolved and coarse-grained regions, allowing the proper exchange of molecules. We numerically prove that indeed for systems as relevant as liquid water and 1,3-dimethylimidazolium chloride ionic liquid, the combined action of the TD force and thermostat allows for computationally efficient and numerically accurate simulations, beyond the current capabilities of adaptive resolution set-ups, which employ switching functions in the coupling region.

## I introduction

Recent developments in the field of multiscale MD methods have brought the idea of adaptive molecular resolution for simulation of liquids to a high standard level, both conceptually and technically jcp-adress (); annurev (); ensing1 (); truhlar (); csany (). In particular, among all methods, the Adaptive Resolution Simulation (AdResS) jcp-adress (); annurev (); pre1-adress () has been intensely developed and applied to a large number of systems in different fields (see recent Refs. matej-shear (); matej-dna (); kurt-nucl-ac (); anim-pccp (); shad-adts () as an example of variety). The root model of the adaptive resolution technique (AdResS) is based on partitioning the simulation box in three regions, namely: a region at atomistic resolution (AT region), a region at coarse-grained resolution (CG region) and, as an interface between AT and CG, a hybrid region ( region) (see Fig. 1).

In the region, the force between molecules and is defined by a space-dependent interpolation of the AT and CG forces: ; is the interpolating function, smoothly passing from the value 1 in AT to the value 0 in CG. Although based on empirical arguments, such an algorithm has been proven to be computationally robust for several highly challenging applications (see e.g. wat1 (); wat2 (); full1 (); prlado ()). Further developments of the root model have removed most of the empirical character and led AdResS through two different but complementary routes: (a) Grand Canonical MD for the AT region (GC-AdResS) and Open Boundary MD (OBMD) for the AT region and coupling to continuum prl2012 (); prx (); njp (); obmd (); obmd-DNA (); (b) Global Hamiltonian for the entire simulation box (H-AdResS) prlraff1 (); prlraff2 (). The former consists of re-framing the force-based technique within a solid conceptual background of Grand Canonical ensemble for AT, with mathematically sound criteria of validity (see in particular Refs.prx (); njp ()). The latter instead introduces a technical variation, that is the space-dependent coupling of the potential instead of forces: ; is equivalent to but acting on , the center of mass of molecules and . Common to both derivations is the addition of a one-body external thermodynamic term that, while conceptually derived from basic principles of statistical mechanics, in technical terms allows for a density balance between the various regions. In GC-AdResS (and AdResS), this term is named a TD force and is calculated self-consistently in the equilibration step before the production run: ; is the mass of the molecule, a (converge-driven, but well defined) tunable constant, is the molecular density as a function of the position in the region at the -th iteration prl2012 (). Conceptually, its application corresponds to equalizing the chemical potential of the various resolutions to the reference atomistic chemical potential jcp-simon (); jctc-han (); prx (); jcp-mu (). Instead in H-AdResS, the density balancing operation is made in terms of a balance of an artificial global free energy (because a global Hamiltonian with a physical meaning cannot exist pre2007 ()), in an adiabatic interpretation of the switching from the AT to the CG potential in the region (Kirkwood Hamiltonian)prlraff1 (); jcppep (). Independently from the specific coupling, both routes go under the guiding principles for systems with open boundaries physrep (). Recent applications of GC-AdResS for extremely challenging systems as imidazolium-based ionic liquids, i.e., a high density liquid mixture of molecular anions and cations, have highlighted the capability of the TD force in establishing equilibrium in the system and produce results of high standard krek (); shadrack (); krekpshad (). As a consequence, a natural question concerns the possibility that, contrary to the pioneering idea of AdResS, the TD force would be actually a more fundamental technical characteristics for the adaptivity than the smooth interpolation of the region. Additional arguments to support the concept that the TD force, as formulated above, has a fundamental role in equilibrating and stabilizing a system, are provided also in recent work of other groups (see, e.g., Ref. hmmm ()).

In this work, we investigate the possibility that the TD force represents the essential technical ingredient in the molecular resolution adaptation process and show that indeed it is able to assure proper equilibrium and stability even in case of a direct AT/CG interface without any space-dependent smoothing (see Fig. 1). In such a case, standard AdResS becomes automatically H-AdResS with and vice versa. Despite that such a scheme implies the existence of a global Hamiltonian, our interest will not be related to this aspect. As reported in physrep (), although in MD Hamiltonian-based algorithms represent a traditional source of genuine rigor, this principle does not necessary hold anymore for systems with open boundaries that exchange matter with an environment. Actually, for AdResS, at a conceptual level the existence of a global Hamiltonian does not imply a mathematical superiority compared to a force-based approach; on the contrary, as anticipated above, a global Hamiltonian implies a unphysical global (micro-)canonical ensemble which may lead to a misinterpretation of results (see also pre2007 ()). Moreover, in actual adaptive resolution simulations of realistic systems, energy dissipation is, for several reasons, inevitable even if one has a well defined global Hamiltonian. For example, it is mandatory to apply a repulsive force capping that prevents atoms overlapping in the region even for molecules of medium size in the standard interpolation scheme. Such a force capping induces dissipation and requires the action of a thermostat for achieving the desired chemico-physical equilibrium jan (). Furthermore, an equilibrating thermostat is mandatory for performing simulations (even full-blown atomistic ones) of any system with a finite cut-off for the interactions, e.g., the reaction field methodneumann () for electrostatics used in AdResS. It must be also underlined that the hybrid or transition region where the two resolutions are coupled is an artificial region by construction, in any AdResS version. For this reason, one is free to operate at technical level in such regions within some generic/macroscopic constraints (see discussion below) to make the algorithm more efficient. The key point of such an approach is that despite the artificial character of the region, physical properties of the AT region do not depend on the pragmatic computational (empirical) parametrization of the transition region. In fact in Ref.njp () AdResS was mapped onto the Bergman-Lebowitz model of open system with the only constraint that the and CG region are, at any time, at the same molecular density and temperature of the thermodynamic point chosen for the AT region. In this sense, (GC-)AdResS for the atomistic region, considered as an open region embedded in a generic thermodynamic environment, is not an empirical approach but has got a real physical meaning at any time. In this perspective, here we will focus on the technical validity of the method in reproducing, in the AT region, results of a reference NVT full-blown atomistic simulation, that is, as underlined above, we will consider the AT region as a system with open boundaries (Grand Ensemble) njp (); physrep (). A major technical advantage of this formulation concerns the simplicity of implementation and transferability of the algorithm in any MD code. Its implementation consists only of (i) the partitioning of the simulation box in two regions, (ii) the definition of the CG molecules in the neighboring list of the AT molecules at the border (and vice versa), (iii) the addition of the external iterative loop for the convergence of the TD force, and (iv) if necessary, finding a proper molecular orientation avoiding molecular overlaps when entering the region from the CG side (see below). A side consequence of such simplicity is the removal of internal loops required by the switching function, (or in H-AdResS), in locating a specific pair of molecules at each time step. Such a removal leads to a sizable gain in computational performances, as it was predicted by a previous computational analysis of AdResS amdahl (). There is also one relevant methodological advantage, namely that the atomistic radial distribution functions (RDFs) are the same as in a full atomistic simulation not only in the AT region but also in the region. This condition automatically satisfies the necessary condition derived in prx () which assures, at mathematical/rigorous level, that the probability distribution function in the AT region is the same of a full atomistic simulation up to a third order (i.e. up to a 3-particle distribution function); this property in MD implies a very high accuracy (within ) in the statistical determination of physical properties in the AT region. Next, we will report the formal structure of the direct coupling and discuss its numerical implementation. We will then show numerical results for liquid water at ambient conditions; in addition, with the current method, we will reproduce some relevant data of a recent published work regarding the solvation and the free energy of aggregation of two micelles in water. The original results were obtained with the previous approach at a higher computational cost compared to the computational effort required by the current approach. Finally, we will show numerical results for 1,3-dimethylimidazolium chloride ionic liquid in the extreme case of uncharged coarse-grained models for the anions and cations.

## Ii Methodology: From smooth coupling to a direct interface

Following Fig. 1, the simulation box is divided in an AT and a CG region. Let us define , the cut-off distance of molecular interactions. In the AT region, all molecules interact via the atomistic potential: , with the atomistic potential between all the atoms of molecule and molecule . Similarly in the CG region, all molecules interact via a coarse-grained potential: , with the coarse-grained potential between the centers of mass of molecule and molecule . In the AT region at the interface, that is in the region of dimension along the axis, named region in 1, AT molecules interact with CG molecules via the coarse-grained potential, i.e. the coupling AT-CG interaction. Since CG molecules do not have atomistic degrees of freedom there is no other way to interact with AT molecules: , with the center of mass of the atomistic molecule and the center of mass of the coarse-grained molecule . Effectively, the region acts similarly to the transition region of standard AdResS. Thus, atomistic molecules closer to the GC region experience the CG character of the interaction more than AT molecules located at larger distance from the CG region; that is the passage from one resolution to the other is not artificially smoothed via as in standard AdResS, but it is implicitly gradual (i.e. function of , not strictly abrupt). Similarly to standard AdResS and following the same numerical procedure, in , we define the TD force for the density balance. As already explained, one of the aims of the switching function in combination with the repulsive force capping and the local thermostat in the standard AdResS is to avoid fatal large forces due to potential overlaps (or unphysical short distances) between atoms of neighboring molecules at the CG/ interface. In the new implementation, the overlaps must be avoided in another way, e.g., by fixing the center of mass of a problematic molecule and running a few MD steps with a repulsive force capping to find an energetically permissible orientation. In a Monte-Carlo implementation, for example, this would be automatically taken care of by the rejection criterion cameron (). In the present work, we adopt the simplest approach, that is we avoid molecular overlaps by simply capping the interactions at close distances. These are all the technical ingredients required for accurate and efficient GC-AdResS simulations of liquids. Conceptually (but not technically), this should be equivalent to the standard approach with the Heaviside step function used as the switching function supplemented with the molecular orientation adjustment.

## Iii Results

### iii.1 Liquid water at room condition

We consider systems of liquid water at room conditions of different sizes, as treated in previous work with AdResS and whose technical details are reported in the Appendix. Here, we report the results for a system of molecules in total. Our aim is to show the accuracy of the scheme in predicting basic properties in the AT region. The necessary condition that any adaptive resolution technique must satisfy is that the particle number density is uniform across the simulation box. Such a goal in AdResS is achieved through the combined action of the TD force and a thermostat so that the thermodynamic state point corresponds to the desired state of reference.

Fig. 2 depicts the convergence of the density to the target value (within a accuracy) as the calculation of the TD force progresses. Once the TD force is determined (see lower panel of Fig.2) the density is maintained uniform for the production run. Compared to the standard AdResS, the new scheme requires, as expected, some more iteration steps, however, their cost is negligible. Next, we must prove that structural properties are accurately reproduced.

Fig. 3 presents the various RDFs calculated in the AT region, results show a highly satisfactory agreement with reference data. Moreover, RDFs represent second order (two-body) approximation of the -particle probability distribution. Hence, the accuracy shown in these results assures at least a second order accuracy in the statistical calculation of physical observable in the AT region (see also prx (); njp (); krek (); shadrack ()). The natural question arising is whether such a high accuracy is an artifact of the scheme, which, by having a direct interface, may actually create a barrier to the proper exchange of molecules from one region to the other. In fact, in such a case, the initial configuration (starting from an equilibrated full atomistic setup of the whole simulation box) remains essentially the same in the AT and CG region separately.

Fig. 4 shows that samples of particles located in the AT or CG or regions properly diffuse into the other regions as the simulation progresses. In addition, Fig.5 reports the probability distribution, , in the AT region and show that it is satisfactorily close (within ) to the target atomistic simulation reference.

Next, a water molecule has a dipole whose orientation in the liquid does not have any preferential direction. However, coarse-grained molecules do not have a dipole and thus one may envision a sharp artificial order at the interface, typical of polar molecules in contact with a non-polar substrate. In turn, this ordering may have strong influence on the liquid structure of the AT region.

Fig. 6 shows the absence of a preferential alignment of water dipoles along a direction in the AT region, where in fact a reference full atomistic simulation is perfectly reproduced. Instead, as expected, at the interface between and CG regions the molecular dipole tend to align along a preferential direction as it happens to liquid water confined by a non-polar surface. Finally, the important role of the thermostat may induce doubts about the utility of the method for the calculation of time-dependent statistical properties since the dynamics in the AT region would be inevitably perturbed by its action on the corresponding molecules. However, this problem was already treated in GC-AdResS by allowing the action of the thermostat in the CG and regions only, so that the dynamics in the AT region is the natural dynamics of a Grand Canonical system njp (); jcppi (). We have implemented the same principle in the current scheme and confirm its capability to reproduce results of the full atomistic simulation of reference regarding quantities that characterize the essential features of the liquid. Specifically, the setup with a thermostat acting only in the coarse-grained region shows, as for the setup with a global thermostat reported before, its capability to reproduce (i) the detailed atomistic local structure of the liquid (radial distribution functions, (Fig. 7)), (ii) the expected exchange of particles in time from one region to another (Fig.8), (iii) the expected particle number distribution in the atomistic region (Fig.9), (iv) the absence in the AT region of artificial orientational effects due to the abrupt change of molecular structure between and the CG region. In the latter case, to offer an alternative view, we reproduce the tetrahedral order parameter, , as a function of the distance (Fig.10). The latter quantity, for a single molecule taken as a reference, is defined as: , with label the four nearest Oxygen neighbors to the central Oxygen atom of the -th molecule, is the angle formed between the vectors r and r. Finally is calculated as the average over all the water molecule. In addition, we report the dipole-dipole time autocorrelation function, , calculated in the AT region (where there is no thermostat acting) and compared it with the same quantity calculated in the equivalent subregion of the full atomistic (NVE) simulation. The dipole auto correlation function is defined as:

(1) |

where calculates the correlation between the electric dipole moment of molecule at time 0 and . Results show a good agreement with the full atomistic simulation of reference and thus it shows the capability of the method to properly reproduce dynamic quantities as well.

Finally, Fig. 12 shows the performance of the current implementation compared to the old implementation for SPC/E water in the GROMACS code gromacs ().

This is a particularly challenging situation since GROMACS is the most efficient MD code with a specific optimization of the full atomistic simulation of this water model. The enhancement in the performance is sizable (see also the discussion about ionic liquids in the section below).

### iii.2 Two micelles in Water

In a recent work by some of us shad-adts (), GC-AdResS has been employed to determine the free energy of aggregation of two micelles in water. This study represents a relevant subject involving a relatively large system and a large number of calculations, thus computational efficiency plays a major role. Fig. 13 reports a schematic illustration of the two different computational approaches compared here: the current version (top) and the older version of the code (bottom), basic technical details of the simulation can be found in the Appendix, further details can be found in Ref.shad-adts ().

Fig.14 shows the (symmetrized) average density of water along the center of mass-center of mass direction of the two micelles for a distance of 6.0 nm.

Such a distance was found to correspond to the radius of the smallest atomistic region which can reproduce results of a full atomistic simulation of reference. Results clearly show that the new approach can reproduce the previous AdResS results and the full atomistic data of reference very well. The computational speed up in such case is equal to a factor of 2.4 w.r.t. the full atomistic simulation to be compared with a 1.4 of the old code. Furthermore, in Fig.15 we plot three relevant points of the free energy of aggregation calculated with the full atomistic approach, with the standard AdResS and with the method proposed in this work. Once again the agreement is highly satisfactory, that is, the most relevant quantity of the study is reproduced with high atomistic accuracy and at the same time at a sizable lower computational price.

### iii.3 1,3-dimethylimidazolium chloride Ionic Liquid

We have chosen 1,3-dimethylimidazolium chloride ionic liquid as a test system of the new code because this is a typical prototype of ionic liquid and it is highly challenging for MD (in particular for AdResS) given the relevance of the long range electrostatic interactions. The setup of the simulation is the same as that of Ref.krek () (see Fig.16), where the previous version of the code was used.

Further technical details of the simulation are reported in the Appendix. In particular, we consider an uncharged coarse-grained model, this latter is highly challenging for the abrupt interface employed by the proposed approach. Similarly to the case of liquid water, we calculate the quantities that characterize the essential features of the liquid and results are reported in Figs. 17 (see corresponding TD force in Figs.18), 19, 20, 21, 22. Similarly to the case of liquid water, also in this case the results show that the method has the same accuracy of the previous version of AdResS and an improvement in the computational performance w.r.t. the previous implementation has been found; specifically a speeding factor of about 1.4-1.5 w.r.t. the previous version of the code.

## Iv Conclusions

In conclusion, we have proven that the action of the TD force in adaptive resolution schemes is technically sufficient for developing a highly efficient and easy-to-implement/transferable algorithm. Moreover, it assures higher accuracy in the AT region. In future perspective, such a scheme resolves also the main concern of a recently proposed adaptive resolution scheme for molecules with electrons cpcel (); adtsluigi (). In fact, in that scheme, the only technical concern is the absence of a switching function for a numerical smooth coupling between the quantum and the classical region. The use of a switching function would make the quantum treatment artificial and uncontrollably unphysical. Thus, one can consider the scheme proposed in this work as an alternative path for the construction of an efficient and transferable numerical recipe for multiscale simulation of last generation.

## V Appendix: Technical Details

### v.1 Liquid Water

As a test system, we simulated a box of water with 6912 SPC/E water molecules and compared the results of NVT runs with full atomistic MD simulations and also with results from the standard GC-AdResS simulations. We used a modified GROMACS 5.1.0 gromacs () code where the box size is a=15.00 nm, b= 3.72071 nm, c= 3.72071 nm. In the AdResS set-up, the extension of the atomistic region is in total 4.4 nm. The hybrid region was set to 1.3 nm (a total of 2.6 nm), that is slightly larger than the cutoff length of 0.9 nm of interactions. The electrostatic interactions are treated by the standard reaction field method. We use, as in the standard GC-AdResS implementation, the stochastic dynamics integrator (Langevin thermostat). The simulation temperature is 300 K and the timestep is 2 fs. The thermodynamic force was calculated via 50 ps NVT simulations. Once the thermodynamic force has converged, we run a 20ns production run. For the case, where a local thermostat was employed, in order to speed up calculations, we kept the same set-up of the full NVT simulation but we considered a smaller atomistic region (1.6 nm in total). The value of the constant is equal to 0.07 in all simulations presented in this work, for any system. Our experience suggests that for liquids made of small molecules using this value of avoids large jumps in the density as the iterative process evolves and thus leads to a smooth convergence within 25-30 iterations at the worst. It is also possible to choose different values and a systematic study of the behavior of the TD as a function of for different systems and sizes is under development. In any case, being an iterative procedure, does not play a key conceptual role. Regarding the simulations used to determine the computational performance of the method w.r.t. to the previous implementations, they were conducted with the modified GROMACS 5.1.0 gromacs () code with systems of 6912, 13824, 20736, 27648, 34560, 41472, 48384 SPC/E water molecules. The total atomistic region was set to 2.2 nm, and the hybrid region to 1.3 nm (i.e. a total 2.6 nm) while the CG region was systematically increased and run for 5000 steps with a timestep of 2 fs. The coarse-grained model of water was derived as in Ref.wat1 () by a straightforward Inverse Boltzmann Iterative procedure (IBI) ibi (). For NVE simulation, we use the reaction field method for calculating the electrostatic interactions in the system,with dielectric constant, , as this tends to give good energy conservation and has been reported in the gromacs manual. For treating the Vander Waals interactions, we use the “switch” method. The cutoff radius for interactions is nm.

### v.2 Micelles in Water

We used the same setup as described in Ref. shad-adts (). The relevant details here are the following: the system consists of a box of nm in linear dimension, containing 190000 water molecules. The atomistic region is a sphere centered in the middle of the box of radius nm containing a hybrid region of a spherical tackiness of nm. For the free energy calculations of the two micelles at the different (chosen) points of aggregation we have run simulations with the above mentioned modified GROMACS 5.1.0 gromacs () version and compared them with the results of the previous study. All technical details for the reproduction of the results can be found in the technical section of Ref.shad-adts ().

### v.3 1,3-dimethylimidazolium chloride ionic liquid

We used the same setup for the ionic liquids simulations reported in Ref. krek (). The only difference is that our new system contained 2000 ion pairs of 1,3-Dimethyl-imidazolium chloride instead of the published 1000 ion pairs. We duplicated the 1000 ion pairs in x direction. We run a 2ns NpT simulation, with T=400K and a Parrinello-Rahman barostat par_rah (). Then we equilibrated the configuration via NVT for 10 ns at T=400K. The resulting box size was a=23.98304 nm and b=c=3.99718 nm. For the GC-AdResS simulations, we used the neutral CG potential already described in Ref. krek (). In the coarse-grained region, the ion pairs are modeled as neutral spherical beads, specifically, one sphere for the cation and one sphere with an effective excluded volume, for the anion. The interaction parameters, as for liquid water, are developed using a straightforward Inverse Boltzmann Iterative procedure (IBI) ibi (). Such coarse-grained model have already been tested in several application of AdResS krek (); shadrack (); krekpshad (). data for analysis are obtained in a production run of 10 ns.

###### Acknowledgements.

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). M. P. acknowledges financial support through the grant P1-0002 from the Slovenian Research Agency.## References

- (1) M. Praprotnik, L. Delle Site, and K. Kremer, J. Chem. Phys. 123, 224106 (2005)
- (2) M. Praprotnik, L. Delle Site, and K. Kremer, Annu. Rev. Phys. Chem. 59, 545 (2008)
- (3) B. Ensing, S. O. Nielsen, P. B. Moore, M. L. Klein, and M. Parrinello, J. Chem. Theory Comput. 3, 1100 (2007)
- (4) A. Heyden and D. G. Truhlar, J. Chem. Theory Comput. 4, 217 (2008)
- (5) L. Mones, A. Jones, A. W. Götz, T. Laino, R. C. Walker, B. Leimkuhler, G. Csanyi and N. Bernstein, J. Comput. Chem., 36, 633, (2015)
- (6) M. Praprotnik, L. Delle Site, and K. Kremer, Phys. Rev. E 73, 066701 (2006)
- (7) J. Sablic, M. Praprotnik, and R. Delgado-Buscalioni, Soft Matter 13, 4971 (2017)
- (8) J. Zavadlav, R. Podgornik, and M. Praprotnik, Sci. Rep. 7, 4775 (2017)
- (9) R. Fiorentini, K. Kremer, R. Potestio, and A. C. Fogarty, J. Chem. Phys. 146, 244113 (2017)
- (10) A. Agarwal, C. Clementi, and L. Delle Site, Phys. Chem. Chem. Phys. 19, 13030 (2017)
- (11) B. S. Jabes, R. Klein and L. Delle Site, Adv.Th.Sim. 1, 1800025 (2018)
- (12) M. Praprotnik, S. Matysiak, L. Delle Site, K. Kremer, and C. Clementi, J. Phys.: Condens. Matter 19, 292201 (2007)
- (13) S.Matysiak, C.Clementi, M.Praprotnik, K.Kremer and L.Delle Site, J. Chem. Phys. 128, 024503 (2008)
- (14) B. P. Lambeth, C. Junghans, K. Kremer, C. Clementi, and L. Delle Site, J. Chem. Phys. 133, 221101, (2010)
- (15) A. Poma and L. Delle Site, Phys. Rev. Lett. 104, 250201 (2010)
- (16) S. Fritsch, S. Poblete, C. Junghans, G. Ciccotti, L. Delle Site, and K. Kremer, Phys. Rev. Lett. 108, 170602 (2012)
- (17) H. Wang, C. Hartmann, C. Schütte, and L. Delle Site, Phys. Rev. X 3, 011018 (2013)
- (18) A. Agarwal, J. Zhu, C. Hartmann, H. Wang, and L. Delle Site, New J. Phys. 17, 083042 (2015)
- (19) R. Delgado-Buscalioni, J. Sablic and M. Praprotnik, Eur. Phys. J. Spec. Top. 224, 2331 (2015)
- (20) J. Zavadlav, J. Sablic, R. Podgornik and M. Praprotnik, Biophys J. 114, 2352 (2018)
- (21) R. Potestio, S. Fritsch, P. Espanol, R. Delgado-Buscalioni, K. Kremer, R. Everaers, and D. Donadio, Phys. Rev. Lett. 110, 108301 (2013)
- (22) R. Potestio, P. Espanol, R. Delgado-Buscalioni, R. Everaers, K. Kremer, and D. Donadio, Phys. Rev. Lett. 111, 060601 (2013)
- (23) S. Poblete, M. Praprotnik, K. Kremer, and L. Delle Site, J. Chem. Phys. 132, 114101 (2010)
- (24) H. Wang, C. Schütte, and L. Delle Site, J. Chem. Theory Comput. 8, 2878 (2012)
- (25) A. Agarwal, H. Wang, C. Schütte, and L. Delle Site, J. Chem. Phys. 141, 034102 (2014)
- (26) L. Delle Site, Phys.Rev.E, 76, 047701 (2007)
- (27) P. Espanol, R. Delgado-Buscalioni, R. Everaers, R. Potestio, D. Donadio and K. Kremer, J. Chem. Phys. 142, 064115 (2015)
- (28) L. Delle Site and M. Praprotnik, Phys.Rep. 693, 1-56 (2017)
- (29) C. Krekeler and L. Delle Site, Phys. Chem. Chem. Phys. 19, 4701 (2017)
- (30) B. S. Jabes, C. Krekeler, R. Klein and L. Delle Site, J. Chem. Phys. 148, 193804 (2018)
- (31) B. S. Jabes and C. Krekeler, Computation, 6, 21 (2018)
- (32) D. de las Heras and M.Schmidt, https://arxiv.org/abs/1711.00006
- (33) J. H. Peters, R. Klein and L. Delle Site, Phys.Rev.E 94, 023309 (2016)
- (34) M. Neumann, J. Chem. Phys. 82, 5663, (1985)
- (35) C. Junghans, A. Agarwal and L. Delle Site, Comp. Phys. Comm. 215, 20 (2017)
- (36) C.F. Abrams, J.Chem.Phys. 123, 234101 (2005)
- (37) A. Agarwal and L. Delle Site, J. Chem. Phys. 143, 094102 (2015); ibid Comp. Phys. Comm. 206, 26 (2016)
- (38) M. J. Abraham, T. Murtola, R. Schulz, S. Pall, J. C. Smith, B. Hess, E. Lindahl, Software X, 1-2, 19, (2015)
- (39) L. Delle Site, Comp. Phys. Comm. 222, 94 (2018)
- (40) L.Delle Site, Adv.Th.Sim. (2018), in press
- (41) D.Reith, M.Pütz and F Müller-Plathe, J. Comp. Chem. 24, 1624 (2003)
- (42) F.Dommert, K.Wendler, R.Berger, L.Delle Site and C.Holm, ChemPhysChem 13, 1625 (2012)
- (43) H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, J. Chem. Phys. 81, 3684 (1984)
- (44) M. Parrinello, A. Rahman, J. Appl. Phys. 52, 7182 (1981)
- (45) C.Hardacre, J.D.Holbrey, S.E.J.McMath, D.T.Bowron, and A.K.Soper, J. Chem. Phys. 118, 273 (2003)