Towards a priori uncertainty quantification in coarse-grained molecular dynamics:
Generalized multipole potentials
Abstract
In computational materials science, coarse-graining approaches often lack a priori uncertainty quantification (UQ) tools that estimate the accuracy of a reduced-order model before it is calibrated or deployed. This is especially the case in coarse-grained (CG) molecular dynamics (MD), where “bottom-up” methods need to run expensive atomistic simulations as part of the calibration process. As a result, scientists have been slow to adopt CG techniques in many settings because they do not know in advance whether the cost of developing the CG model is justified. To address this problem, we present an analytical method of coarse-graining rigid-body systems that yields corresponding intermolecular potentials with controllable levels of accuracy relative to their atomistic counterparts. Critically, this analysis: (i) provides a mathematical foundation for assessing the quality of a CG force field without running simulations; and (ii) provides a tool for understanding how atomistic systems can be viewed as appropriate limits of reduced-order models. Simulated results confirm the validity of this approach at the trajectory level and point to issues that must be addressed in coarse-graining fully non-rigid systems.
Towards a priori uncertainty quantification in coarse-grained molecular dynamics:
Generalized multipole potentials
Paul N. Patrone ^{†}^{†}thanks: Applied and Computational Mathematics Division, National Institute of Standards and Technology, 100 Bureau Drive, Gaithersburg MD, 20899, USA., Andrew M. Dienstfrey ^{*}^{*}footnotemark: *, Geoffrey B. McFadden ^{*}^{*}footnotemark: * |
Scientists have been slow to adopt so-called “bottom-up” coarse-grained (CG) molecular dynamics (MD) in materials development and integrated computational materials engineering (ICME) settings. From a practical perspective, this can be traced to the fact that most, if not all, CG-MD strategies lack a priori uncertainty quantification (UQ) tools that estimate the accuracy of a CG model relative to its atomistic counterpart [1, 2]. As a result, modelers do not know beforehand whether the cost to develop a CG model justifies the end result and have therefore been reluctant to deploy such tools without outside validation. Just as troubling, a variety of studies meant to provide exactly this validation have instead shown that many CG methods suffer from an inability to correctly predict more than one or two material properties, and even then, only after significant recalibration [3, 2, 1, 4, 5, 6, 7]. Seen collectively, these issues paint a confusing picture for the viability of CG-MD and suggest the need for a deeper mathematical understanding of how these methods work.
Conceptually, the lack of a priori UQ can be understood through comparison with coarse-graining strategies used in other fields, e.g. electrostatics. Multipole expansions provide an especially compelling example.[8] In this case, complicated charge distributions are treated as a weighted sum over much simpler geometric configurations (i.e. dipoles, quadrupoles, octopoles, etc.) under far-field conditions. The three main benefits of this approach are that: (i) multipole interactions are pre-computed, so that there is no need to evaluate complicated Coulomb interaction integrals; (ii) truncations (i.e. “coarse-grained” representations) of the multipole expansion have errors that are relatively straightforward to bound; and (iii) the associated expansion converges to the full interaction in as we keep successively more terms. These last two properties are especially important because they imply an analytical connection between the coarse-grained and full model that can be exploited to both quantify and control limitations of the former.
By contrast, many CG-MD strategies tend to rely on ad hoc coordinate reduction strategies that lack a notion of limiting procedure or related connection to the atomistic model.[7, 9, 10, 11] Rather, the key idea is to: (i) define a projection operator that maps the detailed atomistic system onto the desired CG representation; and then (ii) numerically compute a corresponding potential by optimizing an objective that compares the atomistic and CG predictions. As a result, the atomistic system cannot be shown to arise as a limiting case of the CG model, since the latter is chosen by the user and not defined through perturbative or limiting-type arguments. Moreover, notions of trajectory-level convergence between models is lost, since the projection is only enforced in an average sense. As a result, such methods had defied attempts to provide rigorous and a priori UQ to date, especially in the realm of non-equilibrium statistical mechanics.[2]
As a first step to addressing these problems, we present a coarse-graining strategy that reformulates atomistic potentials of rigid bodies in terms of generalized multipole expansions and treats CG models as truncations thereof. The main intuition underlying this analysis is the observation that a rigid body is characterized by only six degrees of freedom, namely its center-of-mass (COM) position and orientation. As a result, the corresponding intermolecular potential collapses to a sum of geometric moments of the rigid body that interact in an orientation-dependent manner. This has the benefit of allowing us to quantify errors in a CG potential without ever having to run a simulation by analyzing truncations of the atomistic potential. Moreover, it provides an analytical route to understanding the limitations of other CG-MD approaches that invoke alternate coordinate reduction strategies. Simulated results at the trajectory level confirm our main results and point to open problems.
We emphasize that our approach, while analytical in nature, has certain limitations. First, accuracy of the method relies on a far-field approximation, which is obviously not valid for nearby particles. Second, our method, as developed, only applies to rigid molecules, which is not a good approximation for all systems. In the discussion section, we pursue these issues in more detail and point towards the appropriate mathematical resolutions.
The rest of this paper is organized as follows. In the next section, we present our generalized multipole potentials for simple power-law systems. The following section presents simulated results, with a discussion and conclusion following thereafter.
I Generalized Multipole Potentials
For simplicity, we consider a system composed of molecules having identical sub-components.^{a}^{a}aHere we mean identical in the sense that all particles have the same interaction law, but not necessarily the same mass. Generalizations to more realistic systems are straightforward but require burdensome notation that is not helpful in an expository setting. To fix notation, we posit that the atomic positions in a molecule are given by
(1) |
where is the atomic position in a laboratory frame of reference, is the corresponding center of mass coordinate, is the coordinate of the atom in a fixed body reference frame, and is a rotation matrix connecting the body and lab frames. It is important to note that Eq. (1) is a coarse-grained representation of a molecule in the sense that is fixed, so that the atomistic coordinates are specified entirely in terms of and . Moreover, only has three independent degrees-of-freedom (DOF), so that each CG molecule has a total of six independent DOF.
To make use of this decomposition, we consider as an example the case in which the interatomic interaction is given by , where is the separation distance between two atoms on different molecules. The intermolecular interaction between molecules and is then given by
(2) |
where () indexes molecules in molecule (), , and , and . In order to make the following expressions more compact, we define and . By analogy with multipole expansions in electrostatics, Eq. (2) can then be expressed as a series expansion of the form
(3) |
where we assume that , i.e. the intermolecular separation is larger than the intramolecular distances. Note that is the gamma function. Interchanging sums over and , we find that the potential can be expressed in the form
(4) |
where the no longer require evaluations over the atomic coordinates.[12]
To see this last point in more detail, consider the term in Eq. (4). One finds that
(5) |
where is the centroid defined by
(6) |
Note that can be precomputed for all molecules, so that Eq. (5) is an inner product weighted by a difference of rotation matrices. The second order term is given by
(7) |
where is a sum of tensor products given by
(8) |
if is expressed as a column vector. Higher order corrections are more tedious to write explicitly and involve higher-order tensor products of the body-frame coordinates.[12]
It is important to note that in using Eqs. (4)–(7) in MD simulations, one requires explicit expressions for the force and torque, which likewise depend on and . Generalized formulas for these expressions have been computed in many works, to which we refer to reader for more details.[12, 13, 14] We also note that this analysis works for any smooth potential, not just power laws. In particular, if is an arbitrary but differentiable interatomic potential, the corresponding generalized multipole potential is given by
(9) |
where denotes the th derivative of with respect to . Simplifying the sums over yields the corresponding multipole expansions.
Ii Simulated Results
To illustrate some of the properties of the multipole potentials, we consider an idealized system composed of triangular rigid bodies with interaction sites located at the verticies. We consider isosceles molecules that have different masses at each vertex, so that the center-of-mass does not coincide with the centroid. In the body-frame, these verticies are located at
(10a) | ||||
(10b) | ||||
(10c) |
where determines the relative length of the two equal sides to the remaining side of the triangle, and is an overall scale factor determining the size of the molecule. Note that as , the molecule reverts to a point-particle, so that together with the “range” of the potential, determines the overall distance scale at which the multipole approximation is valid. For convenience, we set , (where is a mass scale) so that the center of mass is at in the body-frame of reference. Details of the simulation algorithm are provided in other manuscripts.[12, 13, 14]
Figure 1 shows the results of four different simulations of molecules in head-on collisions. We prescribe initial conditions as follows:
(11a) | ||||
(11b) | ||||
(11c) | ||||
(11d) | ||||
(11e) | ||||
(11f) | ||||
(11g) | ||||
(11h) |
where is an adjustable parameter, is a rotation about the lab-frame -axis by an angle , and and are as defined in Eqs. (10a)–(10c). We choose and set the interatomic potential to be a simple power law repulsion. Visually it is apparent that the zeroth-order potential (which is radially symmetric) does not induce rotation, entirely inconsistent with the physics of the exact potential. Increasing the approximation order of the potential yields better agreement with the exact trajectory as shown in the bottom two subplots. Figure 2 shows a related collision in which the two particles are in-plane during the collision. Here and . Notably, the zeroth-order potential fails to predict deflection off the -axis, which is inconsistent with the trajectory of the exact potential. Further analysis of these simulations, including an assessment of the phase-space convergence, will be provided in a manuscript in preparation.[12]
Iii Discussion and Conclusions
As discussed in the introduction, generalized multipole potentials are useful in that they provide an analytical method for coarse-graining rigid-body systems. Moreover, a salient feature of the method as it pertains to UQ is that error in the local force-field evaluations can always be estimated by the first omitted term in the expansion. For example, an approximation to second order in the ratio will yield errors that are , which may often be neglected for large enough separations. We anticipate that such asymptotic approximations may be helpful in constructing more formal estimates of uncertainties associated with bulk properties of many-body simulations.
This being said, usefulness of the multipole approximation also relies on the validity of certain assumptions, which we now discuss. In particular, its accuracy is only guaranteed under far-field conditions. In condensed-matter systems we always expect some subset of molecules to violate this condition, e.g. nearest-neighbors in close-packed systems. However, we note that the coarse-grained representation in terms of Eq. (1) contains all the atomistic information about the system, albeit in a compressed form. As a result, it is always possible to use the full potential given by Eq. (2) whenever the far-field conditions do not hold. In this way, the generalized multipole potential can be tailored to have prescribed levels of accuracy at different separations, so that the overall potential is accurate to within a desired energy threshold.
A second key assumption of our analysis is its restriction to rigid-body systems. In general, many of the most useful MD models incorporate bond-lengths that vibrate about some mean distance in order to represent internal kinetic energy of the system. Such effects are not captured by the potential described in Eq. (2), since the body-frame coordinates are assumed to be constant. Thus, an extension of our analysis is to add a perturbative correction associated with deviations from these average lengths. It is likely that such terms will contribute to higher-order corrections in the potential and couple to translational and rotational DOF in non-trivial ways. Analysis of this scenario is left for future work.
Despite these shortcomings, the generalized multipole potential provides a useful starting point for recasting coarse-grained models as approximations and/or truncations of a fully atomistic model. The key insight we wish to emphasize is the recognition that a priori uncertainty quantification becomes straightforward in such circumstances, as the models become fully connected through appropriate limiting processes. Moreover, such analyses provide benchmarks against which other coarse-graining methods can be understood, given their assumptions in comparison to the multipole potentials.
References
- Dunn, N. J. H., Foley, T. T., and Noid, W. G., “Van der Waals Perspective on Coarse-Graining: Progress toward Solving Representability and Transferability Problems,” Accounts of Chemical Research, Vol. 49, No. 12, 2016, pp. 2832–2840.
- Wagner, J. W., Dama, J. F., Durumeric, A. E. P., and Voth, G. A., “On the representability problem and the physical meaning of coarse-grained models,” The Journal of Chemical Physics, Vol. 145, No. 4, 2016, pp. 044108.
- Patrone, P. N., Rosch, T. W., and Phelan Jr., F. R., “Bayesian calibration of coarse-grained forces: Efficiently addressing transferability,” Journal of Chemical Physics, Vol. 144, No. 15, 2016, pp. 154101.
- Carbone, P., Varzaneh, H. A. K., Chen, X., and Müller-Plathe, F., “Transferability of coarse-grained force fields: The polymer case,” J. Chem. Phys., Vol. 128, No. 6, 2008, pp. 064904.
- Farah, K., Fogarty, A. C., Boehm, M. C., and Müller-Plathe, F., “Temperature dependence of coarse-grained potentials for liquid hexane,” Phys. Chem. Chem. Phys, Vol. 13, No. 7, 2011, pp. 2894–2902.
- Ghosh, J. and Faller, R., “State point dependence of systematically coarse-grained potentials,” Mol. Simul., Vol. 33, No. 9-10, 2007, pp. 759–767.
- Shell, M. S., “The relative entropy is fundamental to multiscale and inverse thermodynamic problems,” The Journal of Chemical Physics, Vol. 129, No. 14, 2008, pp. 144108.
- Jackson, J., Classical Electrodynamics, Wiley, 2012.
- Reith, D., Putz, M., and Muller-Plathe, F., “Deriving effective mesoscale potentials from atomistic simulations,” Journal of Computational Chemistry, Vol. 24, No. 13, 2003, pp. 1624–1636.
- Zhou, J., Thorpe, I. F., Izvekov, S., and Voth, G. A., “Coarse-Grained Peptide Modeling Using a Systematic Multiscale Approach,” Biophysical Journal, Vol. 92, No. 12, 2007, pp. 4289 – 4303.
- Peter, C. and Kremer, K., “Multiscale simulation of soft matter systems - from the atomistic to the coarse-grained level and back,” Soft Matter, Vol. 5, 2009, pp. 4357–4366.
- Patrone, P. N., Dienstfrey, A. M., and McFadden, G. B., “Multipole coarse-graining of rigid-body molecular dynamics I: Trajectory-level perspective on transferability and representability,” In Preparation, 2018.
- Leimkuhler, B. and Matthews, C., Molecular Dynamics with Deterministic and Stochastic Numerical Methods, Springer International Publishing, Switzerland, 2015.
- Dullweber, A., Leimkuhler, B., and McLachlan, R., “Symplectic splitting methods for rigid body molecular dynamics,” The Journal of Chemical Physics, Vol. 107, No. 15, 1997, pp. 5840–5851.