Molecular dynamics (MD) is a useful technique for theoretical investigation of molecular systems such as biomolecular systems and other macromolecular systems.  A primary limitation in the application of MD to the study of complex processes involving macromolecules, e.g. biomolecules, is the small time step size of conventional MD.  Whereas the latter is typically measured in femtoseconds, some dynamical processes of interest happen in microseconds and longer time scales.  Bridging the time scale gap between simulations and the phenomena of interest has been an area of active research for more than a decade.  A variety of techniques have been introduced in order to increase the time step in molecular dynamics simulations. 
One common approach is to constrain bond lengths using either the SHAKE (Ryckaert et al.,  1977) or RATTLE (Anderson, 1983) algorithms.  Although application of these methods allows for a modest (˜ a factor of 2) increase in the time step, time-dependent quantities may be affected.  Additionally, the constraint methods have not been shown to work well for bond angle degrees of freedom when applied to the case of macromolecules (van Gunsteren et al.,  1982). 
An additional complication is that biological and some other complex systems are multi-scale in nature.  For example, the dynamics of proteins contain motions over different time scales, from atomic vibrations in the order of femtoseconds to collective motions at millisecond scales.  FIG.  1 depicts the dynamics of molecules such as protein molecules, to illustrate the variation in time scales. 
Traditional time stepping integrators (e.g. Verlet) commonly used in molecular dynamics (MD) are not able to address this time scale problem.  A typical time-step for these methods is 1 femtosecond.  This makes atomistic simulation of biomolecules computationally extremely expensive.  Multi-scale numerical methods, in which the presence of fast scales does not affect the time integration of slow scales, are urgently needed for efficient simulation of large biomolecular systems.  Such approaches, in theory, can essentially enhance performance of molecular simulation since the most computationally expensive long-range electrostatic interactions contribute to the dynamics on relatively long time scales (compared with internal vibrations) and thus ideally do not need to be calculated frequently.  Also, such approaches enhance the data locality which makes them better suited for implementation on parallel computers than traditional MD schemes. 
One approach is to separate the dynamics into fast, uninteresting modes, and slow, functionally relevant modes and perform MD in the reduced space.  Among the most popular approaches to find reduced dynamical spaces for biomolecules are normal mode analysis (NMA) (Levitt et al.,  1985) and principal component analysis (PCA) (Balsera et al.,  1996).  Those methods have been combined with various computational schemes (e.g. LIN (Zhang et al.,  1993), ACM (Zhang et al.,  2003), LLMOD (Kolossvary et al.,  2001), SMD (Space et al.,  1993), NML (Sweet et al.,  2008)) to yield simulation techniques which, in fact, have not succeeded in either providing the desired accuracy or in achieving substantial computational speed-up (Sweet et al.,  2008) for atomistic simulation of macromolecules. 
Other work has attempted to build multiple time step (MTS) integrators for MD that allow for time steps of differing lengths according to how rapidly a given type of interaction is evolving in time.  The prototypical algorithm is the Verlet-l/r-RESPA/Impulse integrator (Grubmüller et al.,  1991; Humphreys et al.,  1994), which splits the forces into fast (short-range) and slow (long-range) components and evaluates the former more frequently than the latter.  The ratio between frequencies of evaluation of the long-range forces (outer step-size) and short-range forces (inner step-size) measures the gain in simulation time and will be further referred to as “the step-size ratio”. 
FIG.  2 is a diagramatic explanation of the MTS idea, which splits the forces in a system into bonded “fast” forces and long range non bonded “slow” forces (which tend to be non-linear), evaluating the slow forces less frequently.  For this, multiple timestepping integrators are required to solve modified ODEs (Ordinary Differential Equations). 
For biomolecular applications, the computational complexity of the fast short-range force evaluations scales linearly in the number of atoms in the system, N, while it scales quadratically in N for the slow long-range force evaluations.  Furthermore, while the short-range fast forces are easy to compute in parallel, long-range slow forces require global data communication and hence are more difficult to parallelize efficiently.  Therefore, in theory, MTS methods can dramatically speed up MD simulations by reducing the number of expensive slow force evaluations. 
Although in theory MTS methods can dramatically speed up MD simulations by reducing the number of expensive slow force evaluations, in practice, however, impulse MTS methods such as the popular Verlet-l/r-RESPA suffer from severe resonance instabilities that limit practical performance gain (Ma et al.,  2003; Izaguirre et al.,  2001).  For solvated biomolecular systems, for example, the stability limit means that the “the step-size ratio” becomes equal to ˜4.  Performance of impulse MTS methods was recently improved in the Langevin stabilized MTS methods (Izaguirre et al.,  2001) by reducing resonance induced instabilities through the introduction of mollified MTS methods by Izaguirre et al.,  1999 and weak coupling to a stochastic heat bath (Langevin dynamics) (Izaguirre et al.,  2001) to weaken non-linear instabilities.  This allowed an increase of the step-size ratio to 6-12 for solvated biomolecular systems (Izaguirre et al.,  1999, 2001). 
In addition to the limitations which are specific to each listed method, all those methods have a common drawback—they do not exactly sample from the target temperature even if the simulations are stable and are subject to a thermostat (Pastor et al.,  1988; Bond and Leimkuhler, 2007).  This error can be controlled with a loss of computational efficiency by increasing the frequency of slow force updates. 
It is desirable to provide a method and apparatus for simulation which overcome or at least mitigate some of the disadvantages of the prior art.