Accelerating molecular dynamics simulations by a hybrid molecular dynamics-continuum mechanical approach

ABSTRACT In this contribution, molecular dynamics (MD) simulations in combination with continuum mechanical (CM) approaches are performed to investigate particle movements under uniaxial deformations of an amorphous polymer at particle resolution. A coarse-grained (CG) model of atactic polystyrene is used as an exemplary model system. We propose a hybrid molecular dynamics-continuum mechanical (MD-CM) approach to simulate the deformation behavior of the polymer. As a reference, purely molecular dynamics systems are used. The methods are compared with regard to the local displacement of the particles and the global stress-strain behavior of the overall system. The good reproducibility of the system’s mechanical behavior with the hybrid molecular dynamics-continuum mechanical method is shown. Furthermore, it is demonstrated that CPU time can be significantly reduced with the hybrid calculation model.


Multiscale Simulations of Polymers
Polymers offer great potential since their properties can be tailored to specific applications. A major challenge is the numerical investigation of polymers particularly in the case of inhomogeneities as they appear, e.g., in polymer composites or polymer fracture problems. From an engineering point of view, the material response under mechanical loadings is of specific interest, which requires a holistic consideration to capture the large range of resolutions involved. In order to appropriately describe the overall macroscopic material behavior, mechanisms at atomistic scales have to be taken into account, but should be confined only to those regions that are of particular relevance in order to constrain the computational effort. In this regard, socalled domain decomposition techniques are well-suited to run multiscale simulations at a reasonable computational cost, where only distinguished parts have to be treated at atomistic resolution, cf. e.g. Ref. [1] The remaining regions, which should be as large as possible, use coarse-scale techniques as, e.g., continuum mechanics at a significantly less computational cost.
In the simplest case, the regions to be considered at fine and coarse resolution are defined a priori, i.e., are based on expectations about the loadings they will be exposed to in the subsequent deformation procedure. This, however, requires knowledge about how deformation evolves throughout the entire system and which parts will be higher or lower strained. In the case of fracture simulations, a complex strain state occurs, which renders the a priori definition impossible. Figure 1 shows one of our first multiscale studies of polymer fracture, [2] in which we choose a particlebased resolution in the vicinity of the crack tip and a much coarser continuum-based treatment elsewhere. For these simulations, we use our domain decomposition "Capriccio" method, whose original version has been introduced and further developed by our group in recent years, cf. Refs. [3][4][5][6][7] The Capriccio method couples a continuum, described by the finite element method (FEM), with a particle-based domain relying on molecular dynamics (MD) and employing stochastic boundary conditions. [8] Aside from polymer fracture as referred to here, it has shown its capabilities in the context of polymer nanocomposites. [9][10][11][12] As evident from fracture investigations in general, the highest loads appear near the crack tip, such that they will move throughout the domain according to the crack propagation. In this case, it is desirable to introduce adaptivity and to select computationally expensive finescale regions based on the effective load they are exposed to. In this concept, the coarse-scale description can be kept wherever only low strains occur, such that significant parts of the deformation history could be described at the macroscale without any actual need of expensive treatments at high resolution. This, in turn, requires the replacement of the coarse-scale description "on-the-fly" by its fine-scale counterpart in distinguished parts. Aside from fracture problems, comparable requirements may occur in many other cases, e.g., nanocomposites or complex loading situations.

A Novel Hybrid Molecular Dynamics-continuum Mechanical Method
Thus, we propose a technique that partly substitutes MD simulations by continuum mechanical (CM) considerations, which we thus label as the "hybrid MD-CM" method. By doing so, we allow to literally "jump" from one deformation state to another whenever expensive MD simulations are not required in full detail. This approach might be beneficial for a wide range of MDbased deformation tests studying the mechanical properties of polymers and polymer nanocomposites: Bulk polymer systems in such studies require a large number of particles, which leads to high computational cost limiting the possible investigations. However, many effects of specific interest as, e.g., fracture and inelasticity, only occur at large strains, such that the elastic range typically comprising only small strains is not of great interest and it is desirable to reduce the numerical effort in this regime. In this regard, our hybrid MD-CM method enables MD simulations accelerated by a factor of two in the elastic range and hence frees capacities for a more in-depth analysis of large deformations.
In the present contribution, we first demonstrate that our novel hybrid MD-CM method is, in contrast to competing approaches, well-suited to accelerate MD simulations of polymers. Second, we show that our technique does not require a priori knowledge about the material behavior in terms of continuum mechanics, i.e., that we do not need a sophisticated continuum mechanical characterization of the MD sample in advance as we have conducted in other contexts. [13][14][15] Third, we substantiate that the hybrid MD-CM method performs well also in load cases that have not been subject to the calibration procedure, such that universal usage within specific load ranges is ensured.

Competing Approaches
Over the past few years, several concepts have been developed to establish a link between continuum mechanical deformation and atom motions in MD. None of them, however, is well-suited for the multiscale simulation of polymers that we are aiming at. To give a short overview of competing approaches, we briefly introduce three strategies as representatives.
Möller et al. [16] present their FE2AT approach using finite element calculations to provide data for atomistic simulations in terms of appropriate initial and boundary conditions. The methodology comprises three basic steps, namely, i) creating a finite element model of the considered atomistic system, ii) obtaining the equilibrium configuration according to the loading conditions by performing an anisotropically elastic finite element simulation, and finally iii) deriving the atomic positions corresponding to the deformed FE configuration in the elastic regime. The last step assigns each atom from the atomistic representation to its corresponding finite element in the FE mesh and predicts the atom's displacement by interpolating the FE displacements of the associated element nodes. This method, though, has been specifically designed for crystalline systems and thus cannot straightforwardly be extended to polymers: As we show later, assigning atom positions in a polymer based on a continuum mechanical displacement prediction leads to large deviations in terms of the mechanical stress.
A further promising approach combining continuum deformations and molecular trajectories can be found in Ref. [17] which is based on the so-called coarse-grained Parrinello-Rahman (CG-PR) method. The CG-PR method extends the Cauchy-Born rule such that hierarchical multiscale simulations of amorphous materials can be performed, which capture inelasticity. [18] This approach relies on representative unit cells (r-cells), which comprise many atoms and thus represent individual MD systems. Each FE element is assigned to its individual MD system and therefore requires the calculation of an r-cell at each quadrature point. This implies a massive number of atomistic simulations to be conducted in parallel, which results in expensive computations [17] and hence counteracts our aim to reduce numerical cost particularly for deformation states that do not require full atomistic treatment.
A third procedure to bypass computational limitations of standard MD simulations is based on an approach by Maeda and Takeuchi [19] and labeled as athermal, quasistatic (AQS) simulation. [20] It adjusts the particle positions affinely according to the applied deformation and subsequently minimizes the potential energy of the particle system using a conjugate gradient algorithm. To mimic appropriate deformations, the temperature needs to be low, and the deformation rate small. [20] Typically, strain steps in the order of 10 À 4 are used, [20,21] which is by two orders of magnitude lower than those to be employed in our approach and thus would also compromise our aim of reducing the computational effort.

Outline
In this contribution, we introduce our hybrid MD-CM method and demonstrate how aspects of continuum mechanics can accelerate MD simulations by computing a reduced number of MD time steps. We begin in Section Model and simulation methodology with the system preparation of a thermoplastic sample polymer and give an overview of the general deformation framework. Section Uniaxial deformation simulations compares the system behavior under uniaxial tension at fine (pure MD) and coarse (pure CM) scales and introduces the hybrid MD-CM method. Section Simulations discusses and summarizes our simulation results and validates the hybrid MD-CM method.

System Preparation
In general, simulating polymers at atomistic resolution results in high computational cost, restricting the system size to be examined tremendously. For the purpose of our investigations, we used a coarse-grained model of atactic polystyrene (PS) from the study by Qian et al. [22] This allows for a sufficient number of polymer chains of appropriate length to be utilized with significantly less computational effort. Here, substituting a single monomer at its center of mass by a so-called superatom enables the system size to be reduced by a factor of 16 compared to the corresponding fully atomistic model. We use two types of superatoms representing the R and S enantiomers. Among other aspects, the choice of the CG model can be justified by considering the vastly separated length and time scales inherent to polymer systems: In this case, the structure and dynamics at coarser length scales, which approach the chain size, remain mostly unaffected either by the local monomer structure or by high-frequency motions of individual atoms. [23] To describe the interactions between the superatoms in terms of inter-and intramolecular nonbonded and bonded interactions in the coarse-grained model, we use the force field developed by Ghanbari et al. [24] and Qian et al. [22] derived at a temperature of 590 K and a pressure of 101.3 kPa. The force fields were generated by using iterative Boltzmann inversion to convert atomistic radial distribution functions (RDF) into coarse-grained profiles. [24,25] With these force fields, bond interactions up to the second neighbor atoms are considered besides the inter-and intramolecular nonbonded interactions. The torsional component of the potential is negligibly low and has therefore been omitted. [26] The molecular dynamics calculations are performed with the classical molecular dynamics code LAMMPS [27] using message passing with a spatial decomposition of the simulation domain by a 2 x 4 x 5 MPI (message passing interface) processor grid.
The main steps of molecular dynamics simulation are performed according to Ref. [28] , namely, system preparation, minimization, equilibration, and production. The model system consists of 300 PS chains with a length of 200 superatoms per chain. The initial chain structure is created by spatially distributing the polymer chains by using a self-avoiding random walk algorithm based on that implemented by Ref. [24] The cubic simulation box has an initial edge length of 80 nm, so the initial density is low enough to avoid overlapping. To avoid computational instabilities, the energy of the system is minimized in a first step by iteratively adjusting atom coordinates. Therefore, the Polak-Ribière [29] version of the conjugated gradient algorithm is used. In a first step, the initial structure was used to perform an equilibration run to relax high energies and to bring the system to the appropriate state for the production run. For all the following steps, a time step size of Δt ¼ 5 fs and periodic boundary conditions (PBC) are applied. Initially, the systems are equilibrated for 10 ns under the NPT ensemble using a Nosé-Hoover barostat [30] and thermostat [31] at a temperature of 100 K and at atmospheric pressure of p ¼ p atm followed by an annealing run up to a temperature of 590 K for 50 ns in the NPT ensemble to ensure entanglement of the chains. In a next step, the system is cooled down to the desired temperature of 100 K, which is sufficiently far below the glass transition temperature of T g ¼ 170 K in the coarsegrained resolution. [32] Due to the considerably high applied constant cooling rate of 50 K ns -1 , the resulting systems have to be considered as supercooled melts, rather than real solid PS below glass transition temperature. [33] In our considerations, the structure is eventually equilibrated for another 50 ns at a temperature of 100 K using an NVT ensemble. More details of the applied temperature profile and the resulting mass density throughout the equilibration process can be found in the supplementary information.
To ensure statistical reliability, five systems are prepared and averages as well as standard deviations are calculated in the following evaluations.

System Characterization
An important step during the system preparation is to check whether the equilibrated systems are suitable and appropriate for further production runs. To this end, we use the quantities given in the supplementary information, which we classify as thermodynamic quantities, local structure metrics, global structure metrics, and dynamic quantities. [8]

Kinematics of Affine Deformations
In our hybrid approach, we partly use affine deformations of the particle system in the loading direction. Since these are based on continuous field formulations typically employed in continuum mechanics, we label them as continuum mechanical ones. By doing so, we separate them from the discrete and highly non-affine particle motions in purely molecular descriptions of polymer systems. Also note that, as we show below, the hybrid MD-CM method does not require a continuum mechanical constitutive description, but later, possibly inelastic, extensions employing such formulations are envisioned to further accelerate the MD simulations.
According to Figure 2(a) with the undeformed configuration Ω 0 and the deformed one Ω t , the displacement of particle I at time t following an affine deformation of the simulation box reads with R I being the initial particle position vector and r CM I its current one. Following the same notation, the displacement of the same particle obtained from a classical particle-based consideration is Hence, the deviation between the molecular dynamics and the continuum mechanical displacement of particle I follows the below equation:

Molecular Dynamics Setup
Uniaxial deformation is applied by deforming the periodic simulation box in the loading direction (xdirection) by scaling the box length via with L x � 208Å being the initial length of the simulation box, , x being the length after deformation, and λ x being the stretch in the tensile direction, which is linked to the associated strain ε xx by λ x ¼ ε xx þ 1 (cf. Figure 2(b)). In y-and z-directions, the simulation box is allowed to deform freely according to the barostat. In each load step, which comprises a single MD time step, the particles first move according to the affine deformation of the MD simulation box before they adjust their positions according to the subsequent MD run. This procedure naturally requires very small load step sizes in order to enable sufficient equilibration in each step. The temperature of the system and the pressure for the lateral simulation cell surfaces are controlled by a Nosé-Hoover thermostat and barostat with coupling times of 100Δt and 1000Δt, respectively. The global system data are sampled every 10 4 time steps, and a separate trajectory file is written every 10 4 time steps. The production runs are performed for a maximum strain of ε max xx ¼ 8% and different constant strain rates, The total number of time steps is chosen as n ¼ 3:2 � 10 6 with the strain rate predefined by adjusting the time step size Δt as in Ref. [13] This enables a coarser temporal discretization in the case of small strain rates with only moderate particle displacements and has been proven in our previous work to not affect the results substantially.
We are aware of the fact that the strain rates used in the following are significantly higher than those that can be achieved by experiments. Nevertheless, limited by the time scale in MD, suitable results are to be expected. [32] The virial stress tensor σ, which is identified as the continuum mechanical Cauchy stress tensor, is calculated, based on the global pressure tensor p, [34] the atmospheric pressure p atm , and the unit tensor 1.

Molecular Dynamics Reference System
In total, five reference systems are generated according to the procedure described in Section Molecular dynamics set-up. In the MD setup, the particle positions r MD I t ð Þ follow from the MD equations of motion. The overall strain is applied with a constant strain rate governed by (5).

Affine Deformation of the Particle Domain
In the special case of uniaxial tension as considered here, the continuous displacement field reads According to the chosen setup, U 0 i ¼ 0 holds and a maximum strain of ε max xx ¼ 8% is applied. Substituting the coordinate X i of position vector X in the initial configuration by the associated quantity R I i of R I in (7) eventually renders the coordinates of the particle positions, In the tensile direction, we apply ε xx , which is governed by the associated deformation of the MD reference system. In lateral directions, however, we need to specify ε yy and ε zz based on Poisson's ratio, which introduces the necessity of a constitutive description. Based on Ref. [13] we utilize isotropy such that In Section Comparison of particle displacements, we investigate the sensitivity with respect to the choice of ν. Finally, it should be remarked that the present approach is a special case of the Cauchy-Born rule, [35] which reads for multilattice crystalline structures. Here, F is the deformation gradient, � R I is the mean position of atom I, � s are the shift vectors, and ω t ð Þ is the displacement describing time-dependent microscopic fluctuations due to thermal vibrations. In our purely affine deformation approach, we only adjust the particle positions based on the uniform continuum deformation within the MD simulation box. To define a thermodynamic state to the system, the thermostat rescales the velocities of the atoms, which represent the term ω t ð Þ and ensure the target temperature of 100 K on average. Furthermore, � s ¼ 0 holds in our study. As we will discuss in Section Simulations, the Cauchy-Born rule is not suitable to capture the behavior of amorphous polymer systems as considered here since its basic assumption of atoms following a quasiuniform deformation [1] is violated in the present case. According to our findings in Section Simulations, the application of the Cauchy-Born rule to amorphous polymers leads to unphysically high stresses, which occur even in the case of homogeneous uniform deformations. This is in line with the findings of Shan and Nackenhorst, [36] who concluded that the pronounced relaxation processes appearing in polymers drastically reduce the applicability of the Cauchy-Born rule.

Hybrid Molecular Dynamics-Continuum Mechanical Method
In the present section, we introduce our approach to accelerate MD simulations by predicting the particle positions via affine deformations. This technique has first been formulated in Ref. [37] and reduces the computational effort required for MD calculations by predeforming the particle-based system. This is done incrementally by applying load steps with subsequent equilibration phases in each increment.
Setup. Our approach combines molecular dynamics with aspects from continuum mechanics and is thus called the MD-CM method. As introduced in Section Affine deformation of the particle domain, the strain is applied in the deformation step to each single superatom of the polymer system and the box length in the loading direction is scaled according to the applied strain. The difference here is that the strain is applied stepwise in predefined load steps as depicted in Figure 3, whereby the number of load steps is defined by with a given load step size ~ε ls xx . Each load step is followed by an equilibration run of n equil time steps using an NPT ensemble at 100 K and at atmospheric pressure in lateral directions, while the box length in the loading direction is kept constant according to the Stepwise applied strain in the hybrid MD-CM approach compared to strain applied to the reference MD system for (a) time proportional uniaxial tensile loading with a maximum strain ε max xx ¼ 8% and a strain rate _ ε xx ¼ 1% ns -1 and (b) cyclic sinusoidal loading with a strain amplitude ε a xx ¼ 8% and a maximum strain rate j_ ε xx j � 6% ns -1 . applied deformation. As sketched in Figure 5, we define the ratio of equilibration time steps n equil to the total number of MD time steps possible per load step n ls;max as Λ ¼ n equil n ls;max ; (12) which is a coefficient characterizing the MD time step reduction. The strain rate is, as in the case of the pure MD simulations, controlled by adjusting Δt.
As sketched in Figure 4, particle I is located at its initial position R I , which can be formally written as the final position r H ls I end at load step ls ¼ 0 of the hybrid MD-CM approach, as indicated by the superscript "H." The subscript "end" refers to the final position after the MD equilibration phase of the present load step. We assume an instantaneous CM deformation at the beginning of each load step, in which the particles move affinely to the overall box deformation. Hence, the displacement increment ~U H ls I i of particle I at load step ls is proportional to the associated box length increment ~, ls i for each coordinate direction as with the particle coordinates r H lsÀ 1 I i end obtained after the MD equilibration in the previous load step ls À 1 and the associated simulation box lengths , lsÀ 1 i . The box length increments can be formally written as functions of the initial box lengths L i as ~, ls i ¼~ε ls ii L i ; i ¼ x; y; z; (14) such that the particle coordinates before the MD equilibration phase (indicated by the subscript "0") in load step ls read r H ls In this representation, the particle positions r H ls I 0 before the MD equilibration phase only depend on the strain increments ~ε ls ii , the initial box lengths L i , and the box lengths of the previous load step , lsÀ 1 I i . After the MD equilibration, which covers the time span from t ls 0 to t ls end , particle I attains its final position r H ls I end at the present load step. It should be remarked that this position is not the same as would be obtained by a purely affine deformation as discussed in Sections Kinematics of affine deformations and Affine deformation of the particle domain since in each load step, the particle positions relax accord- Figure 5. Stepwise filtering scheme in the hybrid MD-CM method. In each load step, the data points are filtered by averaging over a specific evaluation range defined by �. Furthermore, the equilibration range Λ is introduced as the ratio of equilibration time steps n equil to the total number of MD time steps n ls;max used in the reference pure MD simulation.
ing to the MD equilibration. In the subsequent load step, another displacement increment is applied and followed by an MD equilibration as before. The total displacement of particle I after load step ls with respect to its initial position R I is U H ls I ¼ r H ls I end À R I ; (16) such that the deviation with respect to the pure MD reference solution can be specified as For the sake of simplicity, the superscript indicating the current load step has been omitted in (17).
To conclude the present discussion, three special cases can be identified: 1. Pure MD simulation in equilibrium: If no continuum mechanical deformation is applied, the particles follow the MD equations of motion under the boundary conditions specified to the system. Here, the entire physical time t is considered, cf. Figure 4 for "pure MD".
2. MD simulation under deformation as described in Section Molecular dynamics set-up: In each load step, the particles first move according to the predefined deformation field, which is followed by a single MD time step. This mimics a continuous deformation and requires very small load steps since only a single time step is computed per load increment.
3. Purely affine (CM) deformation according to Section Affine deformation of the particle domain: In this case, no MD equilibration is carried out, such that the entire deformation follows from the continuum mechanical description with the material subjected to rather large load steps. Here, only specific instances in time t 1 0 , t 2 0 , . . . are considered, cf. Figure 4 for "pure CM". Evaluation method. In the MD part of the MD-CM method, the stress data points are sampled every 10 4 time steps. However, they are still subjected to significant fluctuations and, in addition, only the last data points of the equilibration phase of each load step are of specific interest. To be able to compare the results of the hybrid method with those of the reference system, an evaluation method is used as provided by Ref. [13] There, the stress data in each load step are filtered in the ordinate direction by averaging over n eval ¼ �n equil time steps with 0 < � � 1, cf. Figure 5. The result of this filtering procedure represents the stress value at the current load step.

Pure Molecular Dynamics Reference System
As shown in Ref. [13] , the material behavior of the present PS sample system can be considered isotropic, such that it is sufficient to consider uniaxial tension only in the xdirection. As expected for this load case, the components of the Cauchy stress tensor σ in lateral directions vanish except for negligible fluctuations as depicted in Figure 6.

Pure Continuum Mechanical Approach
In a first step, we compare the stress-strain curves of the pure MD reference system and of the pure continuum mechanical approach as introduced in Section Affine deformation of the particle domain. As is obvious from Figure 6, the stress response obtained from the CM approach is by far higher than the reference solution. Furthermore, the associated stress-strain curve is almost linear and tends to be slightly concave, whereas that of the pure MD approach shows a clear convex curvature. Apparently, the polymer chains are not able to rearrange and hence to relax the stress as they do in a pure MD simulation, but are only stretched according to the continuum deformation field. This fits well with the almost harmonic interaction potentials used here, cf. supplementary information, but does not capture the physical mechanisms taking place in the polymer.

Calibration of the Hybrid Method
Based on this, we thoroughly investigate our hybrid approach in the following. In this context, specific attention has to be paid to the parameters introduced in Section Hybrid molecular dynamics-continuum mechanical method, which are summarized in Table 1. In order to find an optimal set of parameters, we run uniaxial tension tests at a constant strain rate of _ ε xx ¼ 1%ns À 1 and obtain the parameter set, which renders the least deviation with respect to the pure MD reference solution. In a subsequent step, we validate our approach for various load cases and temperatures in Section Validation of the hybrid method.
Equilibration time. The dependence of the results on the equilibration time n equil after each load step is significant for the efficiency of the hybrid MD-CM method. A benefit compared to the pure MD calculation only exists if the number of equilibration steps is smaller than the maximum number of time steps in each load step used in pure MD simulations, i.e. if n equil < n ls;max . Figure 7 shows the dependence of the stress-strain behavior of the hybrid MD-CM system on the number of equilibration steps in terms of Λ ¼ n equil =n ls;max compared to the reference system. The other parameters are set constant. Good matching regardless of the choice of Λ for fairly small strains suggests that the MD-CM method is well applicable in the range of 0%-3% strain. For an equilibration time of Λ ¼ 0:5, the stress-strain behavior is reproduced almost perfectly over the entire deformation range. This optimum Λ corresponds to a CPU time of 50% compared to the pure MD reference calculation and thus reduces the computation time by half.
Load step size. The load step size ~ε ls xx affects the polymer system in two ways: On the one hand, a larger number of load steps imply a proportionally increasing number of intermediate equilibration stages, such that in the limit case, the hybrid method would eventually Figure 6. Results of uniaxial tension tests in the x-direction at a strain rate of _ ε xx ¼ 1%ns À 1 : Components of the Cauchy stress tensor from pure MD simulations and comparison of the tensile stresses σ MD xx (MD reference system) and σ CM xx (pure continuum mechanical approach).  Figure 7. Dependence of the stress-strain behavior of the hybrid MD-CM system on the equilibration range Λ at a load step size of ~ε ls xx ¼ 1%, a strain rate of _ ε xx ¼ 1%ns À 1 , and an evaluation range with � ¼ 0: 5. converge to the pure MD treatment. On the other hand, the load step size limits the maximum number of equilibration steps per load step. As the load step size becomes larger, the influence of Λ also increases, whereby large equilibration times reduce stresses, mimicking more a creep test than constant uniaxial deformation. These effects can be observed in Figure 8. For load steps up to 1%, the deviations up to approximately 4% strain are negligibly small. This coincides with the findings for the equilibration time, which also indicates only small deviations in the small strain regime. In the case of larger strains, it might appear strange that smaller load step sizes do not converge to the pure MD reference. This, however, is due to the fact that Λ and � are predefined as 0:5. Indeed, if they were chosen as 1, the pure MD results would be obtained for small ~ε ls xx . Based on the findings given in Figure 8, we identify ~ε ls xx ¼ 1% as the optimum load step size.
Evaluation range. Finally, the effect of the evaluation range � is examined. The parameter � has more the "taste" of a filtering and visualization tool, where the primary role is to filter the stress peak at the initiation of a load step and to prepare visually comparable graphs. Smearing over a wide range of data points might lead to falsified results regarding the actual system behavior. The parameter influence shall therefore play a minor role, which is shown in Figure 9. For small strains, once again, hardly any difference can be observed, and even for strains up to 4%, the curves tend to be within the range of their standard deviations. In summary, the evaluation range has the least impact on the results. Here, it is crucial to choose a suitable sampling time in MD and thus to generate a sufficient number of snapshots (data points) to minimize the influence of fluctuations. Focusing on the area of linear material behavior again, the findings from the preceding parameter studies are confirmed. For the subsequent validation, we choose � ¼ 0:5.

Validation of the Hybrid Method
In the present section, we demonstrate that the hybrid MD-CM method reproduces well deviating loading conditions with the simulation parameter set chosen based on the above findings for time-proportional tests with a constant load rate of 1%ns À 1 . In particular, we use an equilibration range of Λ ¼ 0:5, a load step size of ~ε ls xx ¼ 1%, and an evaluation range of � ¼ 0:5. To this end, we first run the same simulations as above, but with additional strain rates of 0:1%ns À 1 and 10%ns À 1 . Figure 10(a) depicts the associated results and clearly indicates that the hybrid MD-CM method renders appropriate results for the entire range of strain rates. In general, it can be stated that the material behaves stiffer for higher strain rates, suggesting a clear rate dependence and thus viscous material behavior. The specific material behavior has already been investigated in Ref. [13] , and hence, it will not be discussed in more detail here.
Next, we investigate the stress-strain-behavior at 80K and 120K, i.e., � 20K with respect to the calibration temperature of 100K. As depicted in Figure 10(b), both additional curves match very well the pure MD reference data and hence validate the presented MD-CM method also for deviating temperatures. Note that a relation between the dependencies on the strain rate, cf. Figure 10(a), and on the temperature, cf. Figure 10(b), Figure 8. Dependence of the stress-strain behavior of the hybrid MD-CM system on the load step size ~ε ls xx for _ ε xx ¼ 1%ns À 1 , Λ ¼ 0:5, and � ¼ 0:5. could be studied according to the time-temperature superposition principle, [38] which should not be exploited in more detail here.
Third, we consider cyclic loadings as an example for time-periodic loading conditions. We apply one cycle of a sinusoidal loading with different maximum strain amplitudes and a period of t p ¼ 8 ns each, cf. Figure 3 (b). The associated stress-strain behavior is depicted in Figure 11: Predictably, excellent results are obtained at the beginning of all loading cycles, whereas deviations occur in the subsequent stages of the loading and unloading process. Again, these deviations are almost negligible in the case of small maximum amplitudes and become more pronounced for larger amplitudes. Obviously, the parts of the load cycles where the load is reversed require particular attention, presumably a certain adaptivity of the load step size would have to be introduced in such parts with high load gradients. A detailed study of the origin of the deviations and possible remedies is a work in progress and thus beyond the scope of this contribution.
Finally, we focus on the relaxation behavior after the first load cycles are completed. To this end, we hold the system for 50 ns at a strain of ε xx ¼ 0% for preceding maximum strain amplitudes of ε a xx ¼ 2%, ε a xx ¼ 4%, and ε a xx ¼ 8%. By doing so, we eliminate elastic contributions to the stress, such that we can focus on the performance of the MD-CM method with respect to viscous and plastic material behavior. To smooth the data, a Savitzky-Golay filter [39] is used with cubic polynomials and a frame length of 51 data points. Figure 12 depicts the associated results: particularly, the relaxation tests Figure 9. Dependence of the stress-strain behavior of the hybrid MD-CM system on the evaluation range � for ~ε ls xx ¼ 1%, _ ε xx ¼ 1%ns À 1 , and Λ ¼ 0:5. Figure 10. Validation of (a) strain rate dependency and (b) temperature dependency of the reference MD polymer system and the MD-CM system with parameters optimized for _ ε xx ¼ 1%ns À 1 : for ε a xx ¼ 4% and ε a xx ¼ 8% show residual stresses even after long holding times, which indicates plastic material behavior. This is obvious for both the reference system and the hybrid system, but less pronounced for the systems using the hybrid MD-CM method. The overall curvature of the relaxation paths, which describes the viscous behavior, is almost the same for the pure MD and the MD-CM method, such that a good performance of the MD-CM method with respect to viscous effects can be assumed. A significant part of the shift between the relaxation curves of the pure MD and the MD-CM method might be traced back to the deviations occurring already during the preceding sinusoidal load cycles.
Concluding the present section, it can be stated that the MD-CM method introduced here only requires a rather simple calibration procedure under uniaxial tension and is able to produce good estimates for deviating load cases. In particular, excellent results could be obtained for uniaxial tension with different strain rates and temperatures, whereas in the case of time-periodic cyclic loading, the accuracy depends on the maximum load applied to the system. Given a reduction of the total CPU time by approximately 50%, the novel MD-CM method is a promising candidate to accelerate MD simulations especially in the elastic regime. The material considered here exhibits pronounced inelasticity even at low strains, which might constrain the applicability Figure 11. Validation of the hybrid MD-CM method: Cyclic sinusoidal loading with a maximum strain amplitude of (a) ε a xx ¼ 2% and a maximum strain rate of j_ ε xx j � 1:6%ns À 1 , (b) ε a xx ¼ 4% and a maximum strain rate of j_ ε xx j � 3:1%ns À 1 , and (c) ε a xx ¼ 8% and a maximum strain rate of j_ ε xx j � 6:3%ns À 1 . Figure 12. Validation of the hybrid MD-CM method: relaxation tests of the hybrid MD-CM system compared to the pure MD reference system after a full sinusoidal tension-compression cycle with maximum strain amplitudes of (a) ε a xx ¼ 2%, (b) ε a xx ¼ 4%, and (c) ε a xx ¼ 8%.
of this approach in its current form. For materials with less inelastic effects, we expect good results also for larger strain ranges. Beyond this, a detailed study of the interplay between inelasticity and the performance of the MD-CM method is still a work in progress.

Comparison of Particle Displacements
From the MD point of view, the particle movements and the associated trajectories are of specific interest. Hence, the present section compares the particle displacements obtained from our hybrid approach compared to the MD reference system. As a first qualitative estimate, Figure 13 depicts the particle displacements U MD x and U H x in the x-direction in the xy-plane (z ¼ 0) at a strain of ε xx ¼ 8%: For both approaches, increasing displacements with the increasing distance from the origin of the simulation box are obvious, but deviate at a detailed view, c.f. Figure 13(a).
To quantify the deviations, Figure 13(b) shows the absolute deviation ΔU H x in the tensile direction of the particle displacements obtained from the hybrid approach with respect to the pure MD reference simulation. It is obvious that individual particles show large deviations, whereas the majority displays rather small ones. Furthermore, a clear trend of large deviations in specific regions is not observable. This makes it difficult to determine a relative error measure since extremely high relative deviations may occur around x ¼ 0.
Hence, we use the mean squared displacement (MSD) calculated as an ensemble average of the entire system, as a quantitative measure, with N being the total number of superatoms. At a specific load step ls, the MSD gives the deviation of the particle positions with respect to their initial ones. For better comparison, we use the normalized ΔMSD with respect to the pure MD reference MSD MD .
Since the load step sizes differ between the pure MD and the hybrid approaches, we assume that the individual ΔMSD is evaluated at load steps at which the same overall deformation is applied to the systems. Figure 14 compares the ΔMSD for various Poisson's ratios ν in tensile and lateral directions. We display only ΔMSD y as we did not observe any significant deviations between yand z-directions. The ΔMSD curves indicate little deviation, converging toward −1 for the tensile direction and fluctuating around À 2 � 2 for the lateral direction, for all Poisson's ratios.
The deviations obvious from Figure 14 are in the same range of the statistical fluctuations of the pure MD systems. Thus, the hybrid approach almost perfectly captures the pure MD reference independent of the chosen Poisson's ratio ν. This demonstrates that our hybrid method does not require any constitutive formulation to reproduce the MD reference data. This makes it a very promising candidate for accelerating MD simulations for a variety of materials.

Computational Time
Our basic aim is to reduce the computational cost necessary for deformation simulations using molecular dynamics. To this end, we partly predefine the particle motion by continuum mechanics and relax the system in each load step by a subsequent MD equilibration. In Section Hybrid molecular dynamics-continuum mechanical method, we identify the classical load application in a pure MD setting as a special case of our novel hybrid method. As sketched in Figure 4, the hybrid approach applies an instantaneous CM deformation to the system and thus, roughly spoken, "jumps" over a certain number of time steps the pure MD computations would have to consider. Compared to the numerical effort of the MD steps, the computational time of the initial CM deformation is negligible such that the number of subsequent MD time steps governs the computational cost of such a hybrid load step. As demonstrated above, the limit case of pure CM deformation, which would be the least expensive choice without any MD equilibration phase, leads to severe deviations in stresses and particle positions and is thus not feasible. Hence, an appropriate compromise has to be found between the load step size and the equilibration range per load step specified by the parameter Λ. This factor is the ratio of MD time steps effectively carried out to those required in a pure MD simulation, such that Λ is directly proportional to the CPU time required for the hybrid method compared to that for the pure MD approach. As demonstrated in Section Stress-strain response, this factor has to be chosen appropriately in combination with the other parameters of the hybrid MD-CM method in order to balance between reduction of computational cost and maintaining accuracy of the results. As shown in this contribution, hybrid simulations with Λ ¼ 0:5 render reasonable results and bisect the computational time necessary for pure MD simulations.

Summary and Outlook
In this work, accelerated MD simulations of an amorphous thermoplastic under uniaxial loading conditions have been performed using polystyrene in coarsegrained resolution as a model system. In order to ensure statistical reliability, five slightly deviating initial configurations were generated by a random walk algorithm and have been subsequently equilibrated in a pure MD setting. These systems serve as a starting point for further investigations under uniaxial deformations at various loading conditions. We have demonstrated that the particle trajectories do not follow a homogeneous displacement field as one would expect in a pure continuum mechanical setting, and thus, the local and global system behavior differ considerably. Hence, we introduced a hybrid molecular dynamics-continuum mechanical (MD-CM) approach, which appropriately reproduces the global stress-strain behavior particularly in the range of small strains. The correlation between the system response of the pure MD, the pure continuum mechanical, and the novel hybrid setting has been examined in terms of the local (i.e. particle displacements) and the global (i.e. stressstrain response) system behavior. Based on a thorough parameter study, we found that the novel approach reduces the number of MD time steps by a factor of two, such that the overall computational time decreases by approximately 50%.
The validation of the parameter set identified has shown that comparably good results can also be obtained for other load cases in the small strain regime, including cyclic sinusoidal loading and different strain rates as well as different temperatures. However, the authors are aware that further investigations are necessary, particularly not only with respect to further quantities relevant for MD applicants but also in view of an extension to larger strain ranges. Nevertheless, the MD-CM method introduced here already provides a suitable means to significantly accelerate MD simulations for rather complex polymer systems while maintaining relevant local and global aspects of the system behavior. In this regard, our contribution can stimulate further adaptations of this novel approach based on the specific needs.
In our fields of interest, this procedure will become a vital part of the adaptive domain decomposition Capriccio method for multiscale simulations of polymers, which requires to generate a specific deformation state of an MD system without time-consuming recalculations. With such a technique at hand, multiscale fracture simulations tracking the crack tip at fine resolution as envisioned in Section Multiscale simulations of polymers will become possible and most likely will provide new insights into the processes taking place during fracture of polymers.