Challenges in Atomistic-to-Continuum Coupling
This paper is concerned with the design, analysis, and implementation of concurrent coupling approaches where different (atomic and continuous) models are used simultaneously within a single simulation process. Thereby, several problems or pitfalls can happen, for example, the reflection of molecular movements at the “boundary” between the atomic and continuum regions which leads to an unphysical increase in energy in the atomic model. We investigate the problems with the aim of giving an introduction into this field and preventing errors for scientists starting their research towards multiscale methods.
Classical molecular dynamics provides often an adequate description of a physical phenomenon but is computationally expensive. On the other hand, a depiction by a continuum mechanical model combined with the finite elements method is a smaller computational burden but does not go into atomistic details and can be considered less accurate. Multiscale methods which aim at combining fine scale method such as MD with coarse models such as continuum mechanics have been developed. In particular, for elastic systems these methods have reached a considerable interest, since they promise to be a powerful tool in many engineering problems, where only in a highly localized region a precise model (fine scale) is needed and on a larger (surrounding) domain a less accurate model (coarse scale) suffices.
The development of different multiscale methods in different fields started decades ago and has been accelerated in the last few years. Along with this expansion several survey articles have been published in order to classify the multiscale methods by different aspects [1–5].
These methods vary not only in scope and the underlying assumptions but also in their approach to broader questions such as a hierarchical and concurrent multiscale approach. In the first class, the computations are performed on each scale separately. Often, the scale coupling is done by transferring problem parameters; that is, the results obtained on one scale determine the parameters for the computational model on another scale [6, 7]. Thus, for instance, a continuum model can be derived from the atomic information . Another approach is pursued in the concurrent coupling techniques, where the behavior at each length scale depends strongly on the others and an appropriate model is solved on each scale simultaneously, while smooth coupling between the scales is introduced. In the following, we confine our considerations to concurrent methods.
2. Demands on Multiscale Methods and Domain Decompositions
In the following, we show that Domain Decomposition (DD) methods serve as a good motivation for a classification of multiscale methods. To do so, we briefly explain some basic concepts from DD methods.
In the following, we use a different approach by defining DD as certain numerical methods that split the computational domain into two or more subdomains. Although DD methods have been developed for the purpose of achieving concurrency, they can be used in sequential as well as concurrent computations.
In their origin DD techniques have been developed as a powerful iterative method for solving systems of algebraic equations stemming from the discretization of partial differential equations (PDE), that is, from a continuum description; see [11, 12] for an overview. Therein, DD is considered as a decomposition of the finite element space into a sum of subspaces. Then these subproblems are solved by a direct or iterative method. In a next step, projection operators are developed for the information transfer between the subspaces. As a matter of fact, the quality of the approximation on each subdomain depends on the corresponding properties of the approximation subspace. On the contrary, the DD method allows taking benefit of the presence of the subdomains in order to choose the discretization method, which is best adapted to the local behavior of the solution of the PDE which has to be approximated. Thus, the shape of the subdomains and their magnitude of overlap (interface) can be chosen problem dependent. Summing up, the choice of overlapping or nonoverlapping domains and the choice of the transfer operators deeply influence the performance of the DD method. In particular, the choice of an overlapping or nonoverlapping method directly influences the choice of the transfer operator.
Let us now consider the continuum/molecular coupling. Then, the following three points have to be clarified:(i)Definition of the domain.(ii)Design of the coupling region.(iii)Design of suitable transfer operators.
Let us elucidate this DD motivated classification.
2.1. Definition of the Domain
Depending on the type of problem and on the considered domain, the regions with highly local interest have to be defined. For example, cracks and similar defects can involve a global simulation by finite elements and highly localized regions with strong deformations (e.g., crack tips) which are resolved by a molecular dynamics simulation. More precisely, the domain is decomposed bywhere in a fine resolution down to the atomistic scale is employed and in a coarser representation by finite elements is applied. As mentioned previously, since the simulation on is a higher computational burden, the size of plays an important role in the overall performance of any multiscale scheme. As a particular advantage of DD methods the size and shape of the subdomains and thus the interface or handshake region can be chosen arbitrarily. For the different types of interfaces we can distinguish between three cases (cf. Figure 1); these are the following:(i)Interface Methods: The two scales are separated by an interface, which is often a manifold of dimension ().(ii)Handshake Region Methods: The atomistic and the continuum part are matched in an overlap region , , and .(iii)Completely Overlapping Methods: The continuum mechanic description is on the whole domain and the molecular part is a portion of it; that is, the atomistic simulation is everywhere accompanied by the continuum simulation ().
The design of the coupling regions has consequences for the coupling between the molecular dynamics and the continuum mechanics. In nonoverlapping methods, the separation of the fine and the coarse scale is defined by the interface. In other words, the interface is the border between coarse and fine scale such that coexistence in some coupling region can be excluded.
2.2. Design of the Coupling Region
In the coupling region the two descriptions of the matter have to be matched. However, in general, the molecular dynamics is based on a description in some Euclidean space, whereas the continuum mechanics is often given in a weak sense, that is, in a function space. Thus, it is a priori not clear whether the coupling space is a subset of the Euclidean space or of a functions space.
2.3. Space and Transfer Operator
Here, two entirely different models (finite elements and molecular dynamics) are coupled. However, the relationship between their parameters is usually not direct and, thus, care must be taken in the construction of suitable transfer operators in order to relate them. Consequently, any operator has to deal with the incapability of the finite element discretization to resolve displacements down to the atomistic scale. Moreover, the choice of the underlying space for seamless coupling has to be made carefully.
Summing up, the domain of the problem under consideration is decomposed in a coarse and a fine scale part; that is,In the engineering literature, these methods are also sometimes termed as domain-decomposition techniques. In general we expect, for an efficient multiscale or atomistic-to-continuum method (AtC method), that . One prominent example is the crack propagation, where a zone around the crack tip is simulated by molecular dynamics, whereas in the surrounding domain a finite element simulation is applied.
The AtC methods have been applied with success in various engineering applications [13–15]. The different overlapping AtC methods vary in their design of the overlap and the treatment of the interaction between the continuum and the atomistic scale. One of the probably first overlapping AtC methods has been introduced by Mullins and Dokanish, where the stresses in the finite element simulation are computed by an interatomic potential . Kohlhoff et al.  developed the FEAt method, where an atomistic model is surrounded by a finite element mesh with a small overlap region enforcing boundary condition on the atomistic and on the continuum domain. In the Quasicontinuum method [18–20], the total domain is considered and the atomistic region is systematically coarsened out by using kinematic constraints for predefined subsets of the domain. In the Bridging Domain method [21–24], the overlapping is defined via a convex combination of the two descriptions, such that in the handshake region the energy terms are not counted twice. Closely related to this approach is the Arlequin method , which has been initially introduced in the context of continuum-to-continuum blending (cf. Mortar methods ). This approach has been further investigated in the context of atomistic-to-continuum blending . The general blending aspect has been investigated by . In  the coupling has been realized in a function space concept allowing for a transfer between the coarse and the fine scale by using less degrees of freedom in the handshake region . The Bridging Scale method [30, 31], inspired by the variational multiscale method , considers the fine scale part as a subset of the coarse domain (). In order to identify the coarse and the fine scale part in the overlapping region, the total displacement field is decomposed into a fine scale and a coarse scale part. The coarse scale part is the projection of the total displacement onto the finite elements, and the fine scale part is the residue of this projection; that is, the fine scale is that part of the total displacement field which has not been captured by the projection onto the finite element basis. The AtC coupling has also been examined from the mathematical point of view, for example, [6, 33–35]. Recently, coupling methods have been proposed, where the microscale is a canonical ensemble, that is, AtC with finite temperature [24, 36–38].
As a matter of fact, coupling of different physical models is not free from artifacts. One well known and well described phenomenon is the spurious reflections on the interface between the pure atomistic and the mixed region [39, 40]. Here, in this paper, we elucidate the problems or challenges, when coupling continuum mechanics with atomistic models. This paper is structured as follows. In Section 3, we introduce the mathematical setting for the fine scale and for the coarse scale; thereby, we work out that the natural space for continuum mechanics is a function space, whereas for the molecular dynamics the underlying space is the Euclidean space. In Section 4, we derive the different dispersion relations for the atomistic and the continuum scale; thereby, we see that they differ for large wave numbers. In Section 5, we show that a larger mesh size of the finite elements, which is not the same as the distance between the particles, also can cause reflections; thereby, we show a relationship with Shannon’s sampling theorem.
3. Challenge 1: Atomistic versus Continuum Description
When coupling continuum mechanics and molecular dynamics one has to account for the different fundamentals of each theory. Whereas in the atomistic theory every individual particle is considered, in the continuum mechanics an average is considered. In this section we briefly introduce each theory.
3.1. Molecular Dynamics
The material behavior on the microscale is modeled by means of an isolated system of atoms or molecules of a crystalline solid. We identify each of the atoms in their reference position with a point , , where is an index set. Under the influence of external and internal forces, the atoms are displaced in space. The position of the th atom in a deformed configuration is then given aswhere is the displacement of atom .
The atomic displacements are assumed to obey Newton’s law of motionwhere and are the internal and external forces. With each atom , we associate the mass , such thatis the mass matrix on the microscale. If the system is conservative a potential is given; the internal forces acting on a conservative system can be obtained as .
3.2. The Continuum Mechanics Setting
The most common way to classical elasticity is the axiomatic introduction of the relevant quantities and the equations of motion (see, e.g., ). Therein, the atomistic structure of the solid is neglected and the motion of the body is described by functions. Let us consider . In classical elasticity it is assumed that the motion, due to external and internal forces, is given by the mappingsuch that the deformed configuration at time of the body is given bywith the displacement , , .
In the forthcoming we denote points of the deformed configuration byand write instead of and instead of whenever possible. Furthermore, it is assumed that is locally injective in , sufficiently smooth, and orientation preserving, which can be expressed by the pointwise inequalityLet us denote the deformation gradient by ; that is,
The mass is given bywhere is the mass density. Here, we assume that the mass is constant in time, which means that there are no material flows through the boundary of a material subdomain and we do not consider mass to energy conversions. For the sake of simplicity the density is given as a function on by extending by zero outside of .
The change of square length of the line segment due to the deformation is given by the Green strain tensor: . Moreover, the Cauchy stress is defined by , where is the normal and is the traction on any surface segment . For computational and analytical reasons, we use the Lagrangian description and thus define the first Piola-Kirchhoff tensor , which is a second-order tensor, given bywhere is the surface segment in the undeformed (reference) configuration. The governing equations are then derived by conservation of mass, momentum, and angular momentum. They are given bywhere . Here, is the body force density and and . Let us remark that, for the sake of simplicity, we neglected the boundary conditions in our introduction. However, depending on the material model and the pointwise smoothness of the boundary of , the solutions for may not be smooth enough to fulfill (13).
Thus, we introduce the weak formalism. To do so, let be measurable. We denote by the space of -Lebesgue-integrable functions on . Moreover, we denote the scalar product by and the normLet be a multi-index with . Then, for a bounded domain and for , , we denote by the Sobolev space, given byand equipped with the normWe assume that is a polygonal domain and is the disjoint union of the portions and . Moreover, we assume that , , , and the initial displacements . Then for the weak formulation of (13) we seek such that , , andfor all vanishing on , where denotes . Comparing (4) and (17) it becomes clear that there are different methods to combine these two expressions. The interpretation of the molecular description as points in the function space can be realized in different manners. Many of the existing AtC methods constrain the atomistic displacement in the handshake region by an interpolation of the continuum displacement (e.g., [21, 22, 30]). More precisely, the continuum displacement is approximated by finite elements and the atomistic displacement is projected onto the finite element basis functions. Thus, a matching between the two scales is realized in a pointwise sense. In [29, 42], this matching in the handshake region is realized in a function space oriented sense, by embedding the molecular displacement into the .
4. Challenge II: The Dispersion Relation
In the foregoing section, we introduced the molecular description (fine scale) and the continuum description (coarse scale). Here, in this section we consider the velocity at which waves propagate from to or from to . In detail, we want to guarantee that this velocity is conserved when the model and consequently the discretization change (i.e., cross ), since otherwise artefacts at the common boundary can spoil the simulation. In Figure 2, an example of a wave entering a region of slower (faster) wave propagation speed due to a change of the spatial discretization size is shown: The sketch reveals that the difference in the dispersion relation leads to a partial reflection of the wave. To do so, we derive the dispersion relation, giving the dependence of the frequency on the wave number. In this context, we then see that continuum mechanics, as introduced here, with a specific constitutive relation, only represents the linear part of this relation. This difference has fundamental impact for atomistic-to-continuum coupling schemes. More precisely, in an atomistic setting, the dispersion relation, giving the dependence of the frequency on the wave number, is nonlinear whereas in the continuum setting this relation is linear. As a consequence, in the atomistic region the waves with higher frequency move slower than waves with a lower frequency. A typical behavior of such a nonlinear dispersion relation is shown in Figure 3. On the contrary, in the continuum setting, both types of waves have the same velocity. When coupling continuum mechanics and molecular dynamics, this fact has to be taken into account.
In order to derive the dispersion relation, we first briefly introduce the traveling waves in crystals and show then the dispersion relation for the molecular dynamics setting as well as for the continuum setting. For the sake of simplicity, we confine our discussion to the case.
4.1. Traveling Waves in Crystals
Recalling that denotes the deviation of the atoms from their starting configuration a harmonic wave in the atomistic model is of the formwhere is the amplitude, is the phase, is the wave vector, and is the angular frequency. We also define the wave number and the wave length . In the forthcoming, we consider two different kinds of velocities. The phase velocity is defined asMoreover, for the other kind of velocity, we define the wave package as a wave whose amplitude is only in a bounded domain nonzero. The velocity of a wave package is the group velocity, defined as
We can decide between longitudinal and transverse vibrations. An example of a chain in is given by Figure 4.
Remark 1. In the wave vector has three components and points into the direction of propagation. The atomistic displacements associated with a wave where is the equilibrium position of an atom are simultaneously three-dimensional vectors. These displacement vectors may be parallel to (longitudinal), perpendicular to (transversal), or along a direction, that is, not directly related to the direction of .
For the sake of simplicity, we assume that the time dependent motion of the atoms is a linear superposition of harmonic waves. From linearized equation of motion, necessary conditions for nonlinear systems such as typical MD simulations can then be given. In each can be rewritten as , where and is the atomistic spacing. Then, for each , we have the superpositionand we can understand the behavior of the solution by examining the harmonic solutions independently; that is, for all but one .
Let us consider a mass spring system with lattice spacing , mass , and spring constant . Assuming that for this harmonic system a Hamiltonian is given, the respective equations of motion of can be stated aswhich then gives us the dispersion relation for the molecular case (cf. [43, 44])
One can see easily that for any . We thus only consider the case of , which is known as the first Brillouin zone. Note that the open interval is chosen, as for the group velocity equals zero, meaning that the solution is a standing wave. This physical phenomenon is also known as Bragg-Reflection [43, 44]. For symmetry reasons it is sufficient to be restricted to .
4.2. The Dispersion Relation in the Continuum Setting
The continuum mechanics and the molecular dynamics stem from the same basic physical laws, namely, the Hamiltonian or Lagrangian principles. In this section, we consider the dispersion relation for the continuum mechanics. Let us therefore reconsider the molecular dispersion relation given by (23) for small wave numbers . More precisely let be such thatthen for all we haveMoreover, the phase velocity is then given byA connection to the macroscopic elastic properties can be given by compressing the one-dimensional chain from Section 4.1 giving rise to the strain , such that the average distance between the atoms becomes , where .
With respect to the harmonic energy, given by , the energy of the strained chain becomesThus, the strain energy, which can be considered as the extra energy per atom, is given bywhere is the elastic constant. From the velocity given by (26) we haveSumming up, for the long wavelength limit (, ), we obtain
The meaning of the different dispersion relations (23) and (30) becomes clear when the phase velocity and the group velocity are considered. It can easily be seen that the molecular phase and the group velocity are given bywhereas their continuum counterparts are given byThus, for small we haveIn contrast, for large wave numbers, that is, , we have , which implies that in the atomistic setting waves with high frequencies propagate slower than waves with low frequencies, which is not the case in the continuum since its dispersion relation is linear (cf. (30)). This difference in the velocities of the respective models leads to spurious reflections on the interface between the atomistic and continuum description. We have thus motivated that the spurious reflection are not only an artefact of the discretization but a consequence of the difference in the models.
Therefore, each simulation or discretization has to account for this fact. To overcome this difference, absorbing boundary conditions or nonreflecting boundary conditions, like the time history kernel method combined with a scale decomposition , the variational boundary condition , blending methods , and other methods, for example, [24, 29], have been introduced.
5. Challenge III: Discretization of the Coarse Scale
So far, we have shown that the continuum and the molecular scale have different dispersion relations, which carries over to a difference in the velocities of waves. However, for long wavelengths these differences are insignificant. For the numerical dispersion, that is, the case when the continuum is discretized, the mesh size and the shape of the basis functions have to be accounted for.
In order to approximate the continuous displacement field , we now employ a finite element discretization of lower order. Let denote a mesh with mesh size parameter , such that the family is shape regular.
By using Lagrangian conforming finite elements of first order () for the displacements and denoting the set of all nodes of by , the finite element space is spanned by the nodal basisThe Lagrangian basis functions are uniquely characterized by the Kronecker-delta propertywhere is the Kronecker-delta. Any function can uniquely be written aswhere , , is the coefficient vector. We can identify each element of with its coefficient vector . In the forthcoming, we omit the superscript whenever possible. In space, (17) gives a nonlinear ordinary differential equation in which can be solved by, for example, the Verlet algorithm (e.g., ). In other words, we do a spatial discretization followed by a time discretization (method of lines).
Analogously to the atomistic case we consider a finite element approximation for the continuum Hamiltonian of a harmonic system in . By assuming that the material of the atomistic as well as of the continuum model has the same parameters, the respective equations of motion with a spatial discretization can be stated as Hence, if the mesh size equals the atomistic spacing we have that (22) and (37) coincide.
The harmonic wave given by (18) can be interpreted as a signal and the discretization of it as a sampling, which leads us directly to Shannon’s sampling theorem. It states that the sampling frequency of a signal should be at least twice the largest frequency contained in the signal.
In our context, the sampling rate is and the maximum wave length is given by such that Shannon’s theorem becomesThis shows that the mesh size has to be chosen small enough to capture the total atomistic displacement. However, doing so would dramatically increase the computational cost.
Let us assume that the solution of (37) is given componentwise bywhere is the amplitude and is the frequency.
Closely related to the Shannon sampling condition, the mesh size also influences the velocities of the wave. To see this let us consider the numerical dispersion relation for the finite element discretization
In Figure 5 the dependence of the velocity on the space discretization parameter and the wave number is shown. Here the velocity is the normalized phase velocityOne can clearly see that for the wave becomes slower.
Besides the size of the elements, the shape of its basis functions also plays a crucial rule. To see this, let us consider the piecewise linear basis functions.
In molecular dynamics the spectrum of the displacement is large in the sense that there are low and high frequent parts. This fact can cause trouble when applying finite element basis functions based on polynomials. To see this, we observe that the Fourier transformation of our basis function is -band limited; that is, its Fourier transform is less than outside some compact interval . As a consequence high frequent parts of a molecular spectrum being outside of the interval are only represented poorly. In the literature, methods can be found, which attempt to capture these high frequent parts by enrichment basis functions.
In Figure 6, we can see the Fourier transform of the linear finite element basis function (red line).
Obviously, waves with low wave number are better represented (close to ) than high frequencies.
For a given initial amplitude with high and low frequencies the Fourier spectrum is given in Figure 6 (blue, dashed). Moreover, the Fourier spectrum of the corresponding coarse displacement (finite element) is given (red). This problem has been treated in the existing literature differently. Either damping techniques have been applied such that no or damped high frequent waves enter the coarse scale [29, 30, 45] or constraints (Lagrange multiplier) in the handshake region hinder the high frequent waves to enter the coarse scale [21, 29] or the discretization on the coarse scale is enriched .
In this paper, we show that coupling molecular dynamics simulations with continuum mechanics is a challenging task. We describe the respective equations for the two scales and showed that it is not clear at all how to combine them in a sound mathematical way. The spurious reflections which have been observed by several AtC articles are described and their origin is explained: Firstly, the reflections occur due to the different physical models; that is, the dispersion relation, being linear in continuum mechancis and not linear in the atomistic model. Secondly, the different length scale of the two discretizations, finite element mesh size on the one hand and the distance of the atoms on the other hand, causes a reflection which can be explained by Shannon’s sampling theorem. This paper is directed to scientists working for the first time on atomistic-to-continuum methods and it intends to explain the main challenges.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
P. Bochev, M. Christon, S. Collins et al., “A mathematical framework for multiscale science and engineering: the variational multiscale method and interscale transfer operators,” Tech. Rep. SAND 2004-2671, Sandia National Laboratories, Livermore, Calif, USA, 2004.View at: Google Scholar
E. Weinan, X. Li, and E. Vanden-Eijnden, “Some recent progress in multiscale modeling,” in Multiscale Modelling and Simulation, P. Koumoutsakos and S. Attiner, Eds., vol. 39 of Lecture Notes in Computational Science and Engineering, pp. 3–21, Springer, Berlin, Germany, 2004.View at: Publisher Site | Google Scholar | MathSciNet
M. Arndt and M. Griebel, “Derivation of higher order gradient continuum models from atomistic models for crystaline solids,” SIAM: Multiscale Modeling and Simulation, vol. 4, no. 2, pp. 531–562, 2005.View at: Google Scholar
H. A. Schwarz, “Gesammelte mathematische Abhandlungen,” Vierteljahrsschrift der Naturforschenden Gesellschaft in Zürich, vol. 15, pp. 272–286, 1870.View at: Google Scholar
J.-L. Lions and E. Magenes, Nonhomogeneous Boudary Value Problems and Applications, Volume 1, Springer, 1972.
A. Quarteroni and A. Valli, Domain Decomposition Methods for Partial Differential Equations, Oxford Science Publications, 1999.View at: MathSciNet
A. Toselli and O. Widlund, Domain Decomposition Methods, Springer, 2005.
S. Kohlhoff, P. Gumbsch, and H. F. Fischmeister, “Crack propagation in b.c.c. crystals studied with a combined finite-element and atomistic model,” Philosophical Magazine A: Physics of Condensed Matter, Structure, Defects and Mechanical Properties, vol. 64, no. 4, pp. 851–878, 1991.View at: Google Scholar
E. B. Tadmor, M. Ortiz, and R. Phillips, “Quasicontinuum analysis of defects in solids,” Philosophical Magazine A, vol. 73, no. 6, pp. 1529–1563, 1996.View at: Google Scholar
H. Ben Dhia and G. Rateau, “The arlequin method as a flexible engineering design tool,” International Journal for Numerical Methods in Engineering, vol. 62, no. 11, pp. 1442–1462, 2005.View at: Google Scholar
C. Bernardi, Y. Maday, and A. T. Patera, “A new nonconforming approach to domain decomposition: the mortar element method,” in Nonlinear Partial Differential Equations and Their Applications: Collège de France Seminar Volume 11, H. Brezis, Ed., vol. 299 of Lectures Presented at the Weekly Seminar on Applied Mathematics, scientific edition, pp. 13–51, Longman, Harlow, UK, 1994.View at: Google Scholar
S. Prudhomme, R. Bouclier, L. Chamoin, H. Ben Dhia, and J. T. Oden, “Analysis of an averaging operator for atomic-to-continuum coupling methods by the arlequin approach,” in Numerical Analysis of Multiscale Computations, B. Engquist, O. Runborg, and Y. R. Tsai, Eds., vol. 82 of Lecture Notes in Computational Science and Engineering, pp. 369–400, Springer, Berlin, Germany, 2012.View at: Publisher Site | Google Scholar
A. Sadeghirad and F. Liu, “Coupled atomistic-continuum simulations at finite temperature using a three-layer-mesh bridging domain method (TBDM),” in Proceedings of the 55th AIAA/ASMe/ASCE/AHS/SC Structures, Structural Dynamics, and Materials Conference, January 2014.View at: Publisher Site | Google Scholar
H. S. Park, E. G. Karpov, and W. K. Liu, “Non-reflecting boundary conditions for atomistic, continuum and coupled atomistic/continuum simulations,” International Journal for Numerical Methods in Engineering, vol. 64, no. 2, pp. 237–259, 2005.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
P. G. Ciarlet, Mathematical Elasticity; Volume 1: Three Dimensional Elasticity, North-Holland, Amsterdam, The Netherlands, 1988.
C. Kittel, Einführung in die Festkörperphysik, Oldenburg, 2006.
L. Brillouin, Wave Propagation in Periodic Structures, Dover, 2nd edition, 1953.View at: MathSciNet
M. Griebel, S. Knapek, G. Zumbusch, and A. Caglar, Numerical Simulation in the Dynamics of Molecules. Numerics, Algorithms, Parallelization, Applications, Springer, Berlin , Germany, 2004.