ARL-TR-7989 •APR 2017

US Army Research Laboratory

Theory of Wavelet-Based Coarse-Graining Hierarchies for Molecular Dynamics

by Berend Christopher Rinderspacher, Jaydeep P Bardhan,

and Ahmed E Ismail

Approved for public release; distribution is unlimited.
NOTICES

Disclaimers

The findings in this report are not to be construed as an official Department of the

Army position unless so designated by other authorized documents.

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the use thereof.

Destroy this report when it is no longer needed. Do not return it to the originator.
ARL-TR-7989 •APR 2017

US Army Research Laboratory

Theory of Wavelet-Based Coarse-Graining Hierarchies for Molecular Dynamics

by Berend Christopher Rinderspacher

Weapons and Materials Directorate, ARL

Jaydeep P Bardhan

Department of Mechanical and Industrial Engineering, Northeastern University, Boston,

MA

Ahmed E Ismail

Department of Chemical and Biomedical Engineering, West Virginia University, Morgantown, WV

Approved for public release; distribution is unlimited.
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704‐0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the

data needed, and completing and reviewing the collection information.  Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing

the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704‐0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202‐

4302.  Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a

currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD‐MM‐YYYY) 2. REPORT TYPE 3. DATES COVERED (From ‐ To)

4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S) 5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION REPORT

NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)

11. SPONSOR/MONITOR'S REPORT NUMBER(S)

12. DISTRIBUTION/AVAILABILITY STATEMENT

13. SUPPLEMENTARY NOTES

14. ABSTRACT

15. SUBJECT TERMS

16. SECURITY CLASSIFICATION OF:   17. LIMITATION

       OF  

       ABSTRACT

18. NUMBER

       OF  

        PAGES

19a. NAME OF RESPONSIBLE PERSON

a. REPORT b. ABSTRACT c. THIS PAGE 19b. TELEPHONE NUMBER (Include area code)

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18

Apr 2017 Technical Report

Theory of Wavelet-Based Coarse-Graining Hierarchies for Molecular Dynamics

Berend Christopher Rinderspacher, Jaydeep P Bardhan, and Ahmed E Ismail

ARL-TR-7989

Approved for public release; distribution is unlimited.

October 2013-September 2016

FY14-WMR-038

US Army Research Laboratory

ATTN: RDRL-WMM-G

Aberdeen Proving Ground, MD 21005-5066

primary author’s email: berend.c.rinderspacher.civ@mail.mil

We present a multiresolution approach to compressing the degrees of freedom (DoFs) and potentials associated with molecular

dynamics (MD). We suggest a systematic way to accelerate large-scale MD with more than 2 levels of coarse-graining,

particularly for simulation of polymeric materials. We derive explicit models for linear polymers and iterative methods to

compute large-scale wavelet decompositions from fragment solutions. This approach does not require explicit preparation of

atomistic-to-coarse-grained (CG) mappings, but instead uses diffusion wavelets for graph Laplacians to develop system-specific

mappings. Our methodology leads to a hierarchy of system-specific CG DoFs that provide a conceptually clear and rigorous

framework for modeling chemical systems at relevant model scales. The approach is capable of automatically generating as

many CG model scales as necessary, that is, to go beyond the 2 scales in conventional CG strategies. Furthermore, the

wavelet-based CG models explicitly link time and length scales. Finally, a straightforward method to introduce omitted DoFs is

presented, which plays a major role in maintaining model fidelity in long-time simulations and capturing emergent behaviors.

molecular dynamics, multiresolution, emergent physics, simulation, diffusion wavelets, polymers

56

B. Christopher Rinderspacher

Unclassified Unclassified Unclassified UU 410-306-2811

ii
Approved for public release; distribution is unlimited.

Contents

List of Figures v

Acknowledgments vii

1. Introduction 1

2. Methodology 4

2.1 Basic Multiresolution 6

2.1.1 Properties of ∆ ̃ as a Filter 8

2.1.2 How Molecular Information Influences ∆ ̃ 9

2.1.3 Multiresolution Construction from Fragments 13

2.1.4 Example: Linear Homopolymers 13

2.2 Computational Methods 14

2.3 Simulation Results and Analysis 14

2.4 Adaptive Multiresolution 20

2.4.1 Reconstruction Theory 20

2.4.2 Coarse-Graining Hierarchies 21

2.4.3 The Wavelet Hierarchy 22

2.4.4 A Priori Approximations for Reconstruction 22

3. Operator Wavelets Theory 25

4. Trotter Factorization for Use with Operator Wavelets 25

5. Solutions to "Harmonic" Oscillators 27

6. Hermite Functions for Approximating Potentials 28

7. Koopman for MD 30

8. Derivation of Wavelet Spaces 32

9. Error Bounds on Scales 32

iii
Approved for public release; distribution is unlimited.

10. Invariant Subspaces 33

11. Perturbation Theory for Molecular Fragments 33

12. Derivation of Homopolymer Wavelets 34

13. Derivation of Mixed Resolution from Separable Coarsenings 35

14. Conclusions and Outlook 36

15. Publications and Presentations 37

16. Bibliography 38

List of Symbols, Abbreviations, and Acronyms 43

Distribution List 45

iv
Approved for public release; distribution is unlimited.

List of Figures

Fig. 1 Example of a weighted graph Laplacian derived from a simple weighted

graph.....................................................................................5

Fig. 2 Dyadic trees generated by multiresolution formalism using a filter T.

Shaded spaces are final subspaces of the multiresolution. All other spaces

are intermediates of the construction. a) Generic scheme; and b) an

example using butane. In the final spaces in gray, red and blue denote

opposite signs of weights in the construction from the finer scale. .........7

Fig. 3 Chemical structure of polytetrazine ............................................. 10

Fig. 4 Heat map of selected wavelet DoFs in 1,4-D-glucose tetramer. Blue

denotes positive coefficients, red negative. From left to right: fine-grained

wavelet DoF covering a single repeat unit; CG wavelet DoF covering one

half of the oligomer; coarser wavelet DoF covering the oligomer isolating

repeat units with sign changes; coarsest wavelet. ............................ 15

Fig. 5 Selected reduced wavelet representations of a polyethylene (PE) crystal.

From a) to d) wavelet information is successively added. a) Coarse

representation highlights the anisotropy in 1-D chain averaged over all

chain segments; b) Coarse representation of the 2-D chain plane; c) 3-D

representation; d) full resolution. ............................................... 15

Fig. 6 Fourier transform of the y-component of 1,000 atoms in crystalline PE

(100,800 atoms) excluding zero frequency to allow detail at other

frequencies. MD at 500 K and 1 atm. Left: Individual power spectra per

atom. Right: Power spectrum of magnitude of optimal representation. . 16

Fig. 7 Top row: Fourier transform of the y-component of a 100,800 atom

crystalline PE sampled at 1 fs. 3 scales are shown: 1st (left, finest scale),

5th (middle), and 12th scale (out of 25, right). .............................. 16

Fig. 8 Fourier transform of the z-component of alanine dipeptide in vacuum

excluding zero frequency to allow detail at other frequencies. MD at

500 K and 1 atm. Left: Individual power spectra per atom. Right: Power

spectrum of magnitude of optimal representation. ........................... 17

Fig. 9 Fourier transforms of alanine dipeptide DoF time series in vacuum

excluding zero frequency to allow detail at other frequencies. MD at

500 K and 1 atm. a) z-component of optimal representation; b) second

finest optimal wavelet DoF z-component; c) third finest optimal wavelet

DoF z-component; d) fourth finest optimal wavelet DoF z-component. 18

Fig. 10 Top: Heat maps of the DoF correlations of a water dimer in vacuum.

Bottom: Heat maps of the DoF correlations of alanine dipeptide in

vacuum. Left: Correlations of the cartesian DoFs. Right: Correlations of

the wavelet DoFs. .................................................................. 19

v
Approved for public release; distribution is unlimited.

Fig. 11 Dyadic tree generated by separable coarsenings. Shaded spaces are final

subspaces of the multiresolution. All other spaces are intermediates of the

construction with γ1 = μA ⊕ μ

⊥

A and γ2 = μA ⊕ μB......................... 21

vi
Approved for public release; distribution is unlimited.

Acknowledgments

Dr BC Rinderspacher would like to thank Dr Bryan Glaz for his insights and advice

in determining the links between molecular dynamics and Koopman operators. He

also wishes to thank Drs Bardhan and Ismail whose advice and dedication was

indispensable.

vii
Approved for public release; distribution is unlimited.

INTENTIONALLY LEFT BLANK.

viii
Approved for public release; distribution is unlimited.

1. Introduction

Among the challenges faced in molecular simulations is balancing the competing

needs of computational complexity and fidelity to the underlying chemical and

physical phenomena being studied. The challenge is made more acute when the

properties under investigation evolve over time, as many orders of magnitude exist between the fundamental time scale of molecular motions and the time scale

of collective material processes such as heat transfer, diffusion, or elastic deformations. Coarse-graining methodologies have been frequently used to reduce the

complexity of molecular dynamics (MD) simulations, by reducing the number of

“particles” being studied1–5. Although the elimination of degrees of freedom (DoFs)

often leads to minor gains in the time steps that can be used in dynamic simulations6

, such improvements have not led to significant breakthroughs in simulation

capabilities. More advanced coarse-graining techniques are required to enable simulations of systems exhibiting structure on the micron or submillimeter scales, such

as semicrystalline materials.

In this report, we present a new approach toward that end, focusing on polymeric

materials composed of homopolymers or block copolymers. Our approach is based

on the concept of diffusion wavelets 7,8, which both automatically identifies chemical structures that can be reduced into coarse-grained (CG) units and also allows

for repeated application, thus providing much greater customization of the level of

simulation. Such a technique enables the potential reduction of the number of DoFs

by orders of magnitude, as well as for a vastly increased time step, which, when

taken together, permits vastly extended time scales to be simulated using currently

available computational resources.

Among the shortcomings of typical coarse-graining approaches has been that they

normally involve only 2 levels of description, the “original” atomistic level and the

coarse representation, which has typically been based on the developer’s judgment

or intuition. These techniques rely on partitionings of the atoms as the foundation

to deriving coarse DoFs, which typically replace the atomistic groupings with a single “bead”-like entities9–15. While approaches such as the force-matching method

of Voth et al.16 and the reversible coarse-graining method for phenolic polymers

by Kremer et al.9,10,17 have greatly impacted accessible simulation scales for these

materials, they do not fundamentally change the underlying computational process.

1
Approved for public release; distribution is unlimited.

Furthermore, these CG approaches introduce the difficult problem of generating

consistent atomistic reconstructions, because particle-like CG beads lose information about the particles they subsume, which has to be recaptured by other means

(e.g., the use of dummy variables as in dissipative-particle dynamics with energy

conservation18).

The present work focuses on polymeric materials, which are the target application

for many classes of coarse-graining methodologies. In particular, polymers’ long

chain lengths and low defect concentrations compound the problems of 2-scale representations and fine-scale reconstruction, by requiring a large number of particles

and their associated DoFs. Moreover, the equilibration time for polymeric liquids

is itself computationally challenging: a melt whose chains contain N beads each

will require O(N3) time to equilibrate, making simulation essentially impossible

without advanced simulation approaches 19. However, coarse-graining need not be

restricted only to polymers; even relatively “simple” fluids, such as water, have been

the subject of CG models.

Many of the previously mentioned problems can be alleviated through the selection of an appropriate basis set for describing the internal structure of individual

molecules. Consider the analogy of a time-varying quantity f(t). One naturally

represents f(t) in terms of an infinitely local basis (Dirac delta functions). A more

sophisticated approach, suitable for analyzing certain average properties of the signal, might employ the Fourier basis, sines of varying frequency. However, the latter

basis is global in nature, challenging processing and storage for signals that are

infrequent or have sharply varying features only over small durations in time. In

signal processing, wavelets are often used as a basis to differentiate between local

and increasingly global features of the signal, because wavelet bases can be flexibly

defined to efficiently capture features of varying localization20,21. Few other models provide the on-the-fly adaptivity required for important problems, which in the

structural modeling sense might be problems including crack initiation, crack propagation, and interfacial phenomena.22–24. By analogy to the time-varying signal,

purely atomistic models and standard 2-scale CG models represent infinitely local

Dirac distributions, which are expensive; coarse-graining methods based on globally periodic functions (e.g., plane waves) are inefficient for modeling localized

properties. In the context of time-varying signals, the mathematical inefficiency can

be quantified precisely: except for limited special cases, either basis requires in2
Approved for public release; distribution is unlimited.

finitely many coefficients to approximate a signal to within tolerance of a given

error metric (e.g., the l2 or l∞ metrics). In contrast, the number of coefficients required for a wavelet-based representation is usually O(log2 N), where N is the size

of the signal (e.g., number of time samples).

Wavelet ideas have already been used extensively to analyze time-series data 21.

Our application of wavelet ideas to structural representation extends the work of

Ismail25,26. In particular, the earlier work employed Monte Carlo and did not capture dynamical information. The work here develops an approach suitable for MD

simulation by interpreting the wavelet transform of Ismail in terms of the equations

of motion. The method provides a consistent and systematic framework to derive

multiple levels of model resolution while also reducing simulation complexity.

This approach has numerous advantages, whose theoretical basis we address in this

report; follow-on work will illustrate more concrete applications. The general theme

of the advantages is that for the dynamical and nonequilibrium metrics of interest,

this approach especially captures molecular information relevant to both thermodynamics and kinetics. In particular, application of diffusion wavelets to the chemical

topology of the molecule leads to the identification of what we call collective action

modes (CAMs) that represent coordinated motions within the molecule at various

length and time scales. Our approach is rigorously tied to the underlying physics

and has the potential to increase simulation size and duration by several orders of

magnitude. Moreover, the approach is agnostic to the kind of material being studied

and can be applied both to structures of arbitrary chemical complexity, including

both relatively simple molecules such as water or benzene, as well as more complicated chain molecules such as polymers and biopolymers. Finally, it is also capable

of capturing minor effects such as mass effects from chemical substitution (e.g.,

partial deuteration or fluorination). This may be of special importance for materials

design, where the task is to link macroscopic behavior (e.g., Young’s modulus) to

the atomic structure of the monomer unit.

We proceed by introducing the type of classical Hamiltonian model we seek to

accelerate, as well as the wavelet methodology and its basic properties. We then

discuss how wavelets can be applied to uncover the CAMs for a given molecule

and how to use CAMs to systematically reconstruct finer resolutions, as well as

deriving mixed resolutions. We present several examples showing the application

3
Approved for public release; distribution is unlimited.

of the method for small molecules such as water and hydrogen cyanide (HCN)

before summarizing our conclusions.

2. Methodology

We begin by analyzing the equations of motion in MD for a physics-based wavelet

construction. The foundations of MD lie in the application of Newtonian mechanics

to the energy functional

E =

1

2

tr(  ̇x

TMx ̇) + V (x), (1)

where x ∈ R

N×3

are particle positions, x ̇ ∈ R

N×3

are particle velocities, N is

the number of particles, M is the diagonal matrix of particle masses, and V is

the potential. For the macromolecular systems we are interested in, V is usually

partitioned as

V = Vharmonic + Vnon−harmonic, (2)

where Vharmonic(x) = P

i6=j Kij (kxi − xjk − r

(ij)

0

)

2/2, Kij is the force constant of

the harmonic oscillator, xk is the position of particle k, and r

(ij)

0

is the equilibrium

distance of particles i and j. The atoms and the bonds between them define a graph

in which the atoms are the vertices and an edge between atoms i and j has weight

Kij . The maximum number of edges for a vertex in organic materials is 4, and

even in organometallic complexes, the number of edges is unlikely to exceed 6 (the

typical maximum coordination number). Consequently, the matrix representation

of this graph, defined by K, should be highly sparse.

Our multiresolution approach begins with the graph Laplacian of the weighted

graph defined using particles as vertices and the bonds as edges weighted by the

harmonic force constant. This graph Laplacian, denoted ∆, is here a matrix that can

be defined as

∆ := diag(K1) − K, (3)

∆ij =

(P

j Kij , i = j

−Kij , i 6= j

(4)

where 1 denotes a vector of all ones and diag(v) denotes a diagonal matrix defined

4
Approved for public release; distribution is unlimited.

in terms of a vector v, so that if A = diag(v) then Aii = vi. An example of the

weighted graph Laplacian ∆ for a linear triatomic molecule can be found in Figure

1.

Using Eq. 1 and our definition of ∆, we derive the equations of motion

Mx ̈ = −∇V (x) = −∆x − V

0

(x), (5)

where V

0

(x) = ∇V (x) − ∆x is the force for x not due to ∆. ∆ is the graph

Laplacian of the weighted graph with particles as vertices and the bonds as edges

weighted by the harmonic force constant (Figure 1).

1 2

3

K12 = 1.3

K23 = 2.1 → ∆ =

1 2 3

1 1.3 -1.3 0

2 -1.3 3.4 -2.1

3 0 -2.1 2.1

Fig. 1 Example of a weighted graph Laplacian derived from a simple weighted graph

By letting r = M1/2x, p =  ̇r, and q = ∆ ̃ 1/2M1/2x where ∆ =  ̃ M−1/2∆M−1/2,

the harmonic energy term can be expressed as (kpk

2 + kqk2

) /2 and the equations

of motion become

p ̇ = −∆ ̃

1

2 q − V

0

(x) = −∆ ̃

1

2 q − V ̃ (r), (6)

where V ̃ (r) = V

0

(M−1/2r) is the effective potential in terms of r. If V ̃ = 0, the

system can be solved analytically,

y(t) = e

Υt

y0, (7)

where

Υ = 

0 ∆ ̃

1

2

−∆ ̃

1

2 0

!

, y(t) = 

q(t)

p(t)

!

, and y0 =

 

q0

p0

!

. (8)

It can be shown that ∆ ̃ is positive semi-definite, and that the eigenvectors of the

exponential operator in Eq. 7 are (Uj, ±iUj ), where Ujis an eigenvector of ∆ ̃ .

5
Approved for public release; distribution is unlimited.

The solutions to Eq. 7 oscillate with frequencies ±ωj, respectively, where ω

2

j

is

the eigenvalue of ∆ ̃ corresponding to the eigenvector Uj. For this analytical case,

any eigensolution of Υ can be propagated through time independently of any other

solution. Unfortunately, this simplicity is in general broken by V ̃ , which nonlinearly couples all of the eigensolutions to one another. As a consequence, all solutions have to be simulated concurrently according to the highest frequency associated with an eigenvector. Hence, the transformation from particle space coordinates

(x, x ̇)

†

to the harmonic solution coordinates offers no advantage. To circumvent

this issue, a basis that isolates the coupling effects from high- and low-frequency

components is needed.

The key motivation for our work is the recognition that the mass-weighted graph

Laplacian ∆ ̃ relates spatial coordinates to temporal frequencies, which suggests that

its eigenvector matrix is a promising basis for compression. The weighted graph

Laplacian ∆ ̃ and its matrix of eigenvectors are then analogous to the Laplace operator and the Fourier transform, respectively, in conventional wavelet theory.

2.1 Basic Multiresolution

Here we introduce the wavelet transform used in this work and its derivation. We

use the multiresolution analysis for diffusion wavelets as introduced by Coifman

and Maggioni 7. In essence, the multiresolution decomposition partitions the eigenvalues and eigenvectors of ∆ ̃ , effectively strongly coupling high frequencies in the

time domain to high-frequency eigenvectors of ∆ ̃ in the “particle” domain. This

is an important point for the applicability of our approach: not only can DoFs be

reduced, but the time step may also be increased considerably, approximately by a

factor of 2 at each subsequent resolution.

The multiresolution scheme (Figure 2) relies on a positive-definite low-pass filter

T with kTk∞ = 1 (i.e., an operator that suppresses high-frequency vectors) and an

accuracy operator Pε that projects eigenvectors of a matrix X ∈ span{T

2

n

|n ∈ N}

with associated eigenvalue less than a given accuracy ε > 0 to zero. The effect and

purpose of the filter is to project out DoFs associated with high frequencies, thereby

producing a hierarchy of CAMs in respective vector spaces. After each application,

the associated frequencies roughly halve, and concomittantly, the minimum time

step size roughly doubles.

6
Approved for public release; distribution is unlimited.

a) V0

V1

V2

V3

Pε(T

2

)V2

W3

P

⊥

ε

(T

2

)V2

Pε(T)V1

W2

P

⊥

ε

(T)V1

Pε(T)V0

W1

b) CH3CH2CH2CH3

•

•

•

•

•

•

•

•

• •

• • •

•

•

•

•

•

•

•

•

H

H

H

,

H

H

H

,

H H

,

H H

,

H

H

H ,

H

H

H ,

C

H

H

H

,

C

H

H

H ,

C

H H

, C

H H

,

Fig. 2 Dyadic trees generated by multiresolution formalism using a filter T. Shaded spaces are

final subspaces of the multiresolution. All other spaces are intermediates of the construction.

a) Generic scheme; and b) an example using butane. In the final spaces in gray, red and blue

denote opposite signs of weights in the construction from the finer scale.

7
Approved for public release; distribution is unlimited.

The recursively defined vector spaces Vn = Vn+1 ⊕ Wn+1, where

Vn+1 = Pε(T

2

n+1 )Vn (9)

and

Wn+1 = ker Pε(T

2

n+1 |Vn

) (10)

are iteratively associated with orthonormal bases via QR-decompositions,

T

2

n

= QnRn, (11)

where Qn is unitary, Rn is upper triangular, and both are dependent on the basis

used to express T

2

n

.

One advantage of this approach is the inherent invariance of the wavelet spaces with

respect to permutations of the graph vertices and its ability to deal with arbitrary

matrices T with kTk∞ ≤ 1. The wavelet space at scale n (Wn) is approximately

spanned by eigenvectors of T whose eigenvalues obey

2

−n+1 ln ε < ln λ ≤ 2−n

ln ε. (12)

Only a limited degree of mixing with eigenvectors associated with eigenvalues outside of these bounds is possible (see Section 9 for details). Since T is positive semidefinite, T can be rewritten as an exponential e

ln T

. The application of the previous

procedures are equivalent to doubling the spectrum of ln T followed by projecting

out the high-frequency components of the spectrum of ln T.

2.1.1 Properties of ∆ ̃ as a Filter

With an efficient means of computing the wavelet transform in hand, we consider

exclusively in the following the low-pass filter T = I −∆ ̃ /C, where C is a constant

sufficiently large to render T positive semi-definite. To minimize the number of

matrix-matrix multiplications, in particular, with ker Pε(T

2

n

)|Vn = ∅, the normalization constant C = k∆k∞ would generally be optimal as Ψ1 6= ∅. Eigenvalues

for graph Laplacians, such as ∆, are known to lie in [0, 2 maxi ∆ii].

27 Hence, C is

chosen to be between maxi ∆ ̃

ii and 2 maxi ∆ ̃

ii for all numerical examples.

As implied by Eq. 12, the frequency range ∆ ln λ = 2−nln ε for each wavelet space

Wn is drawn ever tighter with each iteration while the eigenvalues λ approach 1.

8
Approved for public release; distribution is unlimited.

Thus generally, more CAMs correspond to scales of high frequency than scales of

low frequency. As a result, each successive scale corresponds to an increase in the

minimum time step size of that scale in MD simulations of approximately ε

−2

n

as

well as a significant reduction of CAMs. Furthermore, unless the effective potential

V ̃ couples involved scales strongly, sufficiently coarse scales (n  m) are quasistatic compared to a given scale m, while sufficiently fine scales (n  m) only

influence the target scale via their static mean. Therefore, only the relevant scales

need to be propagated through time, reducing the DoFs and allowing us to raise the

time step size to match the scales.

Assuming that eigenvalues of T are distributed approximately exponentially, ε =

e

−1/2 would lead to log2 N scales, where N is the dimensionality of ∆ ̃ . This leads

to many DoFs per scale due to issues discussed in the following sections. Instead,

a higher resolution of ε

1/2

machine is used. Although this wastes some computation on

the first few iterations because ker PεT

2

n

|Vn = ∅, it is equivalent to choosing a

tolerance on the scale of δmax, εeffective = ε

2−m

machine ≈ 1 − δmax.

2.1.2 How Molecular Information Influences ∆ ̃

Since the filter T and the weighted graph Laplacian ∆ ̃ derived from the MD potentials share the same eigenvectors albeit with reversed order of eigenvalues, we

discuss its properties in greater detail. These properties have a major impact on the

performance of the wavelet transform. In the following, we discuss shortly the conditions under which separate groups of atoms are strictly independent of each other

leading to highly localized CAMs.

Due to the small number of bonds an atom generally has, vertices are also generally of low degree. As a result, there are highly localized modes due to invariant

subspaces of chemical motifs. For example, any methylene (CH2) group has an associated medium-frequency, highly localized eigenvector of ∆ ̃ . This follows from

the fact that the hydrogens are leaves on the graph (i.e., the weighted subgraph

Laplacian

∆ ̃

CH2 =





2KCH

12 + o −

K

√CH

12 −

K

√CH

12

−

K

√CH

12 KCH 0

−

K

√CH

12 0 KCH



 , (13)

where o collects the contributions from further bonds with the carbon atom, shows

9
Approved for public release; distribution is unlimited.

only contributions from CH2 for the hydrogens). If the hydrogens are identified with

indices i and j on the full graph Laplacian, then (ei −ej )/

√

2 is an eigenvector with

an associated frequency pKCH/mH of a CH vibration.

A slightly more involved example that also shows that this phenomenon is not restricted to leaves on the graph is the repeat unit of the energetic polymer poly1,2,4,5-tetrazine (Figure 3), which has an invariant subspace spanned by 2 independent eigenvectors. The repeat unit block of the weighted graph Laplacian ∆ ̃ is

A =





2X+Z1

12 − √

X

12·14 0 − √

X

12·14

− √

X

12·14

X+Y

14 − √

Y

12·14

0 − √

Y

12·14

X+Y

14 − √

X

12·14

− √

X

12·14

X+Y

14 − √

Y

12·14 0

− √

Y

12·14

X+Y

14 − √

X

12·14

− √

X

12·14 0 − √

X

12·14

2X+Z6

12





, (14)

where X is the bond constant for the C–N bond, Y for the N–N bond, and Z1/6 for

the contributions of vertices outside the repeat unit. The 2 vectors(0, 1, 1, −1, −1, 0)

and (0, 1, −1, 1, −1, 0) span a subspace invariant under application of A but which

is mapped to zero for operators B with entries Bij 6= 0 ⇐⇒ i/j 6∈ {2, 3, 4, 5}.

Hence, these vectors are highly localized (see Section 10 for details of the general

case).

N N

 N N !

n

Fig. 3 Chemical structure of polytetrazine

In large linear homopolymers, discussed in later sections, the invariant subspaces

represent highly degenerate eigenvalues due to the polymer’s repetitive structure.

Degeneracies in ∆ ̃ reduce the effectiveness of the wavelets because no choice of

accuracy can be used to separate them into subscales. In such cases, it may be

possible to incorporate more information from the potential V ̃ , but this is outside of

the current scope.

As an example, butane (H(CH2)4H) has 14 atoms and, therefore, the graph Lapla10
Approved for public release; distribution is unlimited.

cian associated with H(CH2)4H is the 14-dimensional square matrix

∆ ∝





1 −1

−1 3 + κ −1 −1 −κ

−1 1

−1 1

−κ 2 + 2κ −1 −1 −κ

−1 1

−1 1

−κ 2 + 2κ −1 −1 −κ

−1 1

−1 1

−κ 3 + κ −1 −1 −1

−1 1

−1 1

−1 1





,

(15)

where κ = 44/31 is the ratio of force constants for a CC-bond to a CH-bond in

the polymer-consistent force field (PCFF)

28 and the units of kcal/mol/Å2 have been

subsumed in the proportionality constant. The individual CH2 repeat units have

been boxed for emphasis. Since 660 kcal/mol/Å2is an upper bound for ∆, using

the filter T = I − ∆/(660 kcal/mol/Å2) and εmachine as the cutoff criterion yields 6

scales. The first 4 applications (T, T2, T4, T8, T16) did not lead to any unit vectors

below the threshold. At n = 5, the 4 highest frequency modes of ∆ (λ = 517.1,

524.0, 556.3, and 556.3 kcal/mol/Å2, respectively),





H C H H C H H C H H C H H H

−0.50 0.46 −0.50 −0.50

0.33 −0.53 −0.53 0.21 −0.36 −0.36

−0.24 0.35 0.35 0.35 −0.52 −0.52

−0.46 0.50 0.50 0.50





,

(16)

approximately map to zero for T

32, that is, the logarithm of the expectation values

with T do not exceed 2

−5

ln ε ≈ −1.13 (log(1 − λ/660) < −1.53) . These CAMs

correspond to the symmetric stretches of the methyl groups.

The 6 second-highest frequencies (λ =440 kcal/mol/Å2for all, log(1 − λ/660) ≈

11
Approved for public release; distribution is unlimited.

−1.10),





H C H H C H H C H H C H H H

0.81 0 −0.39 −0.43

0 0.72 −0.69

0.10 0 −0.70 0.70 0

0 −0.71 0.71

0 0.82 −0.41 −0.41

0.71 −0.71





,

(17)

are computed for n = 6 and their expectations with T are between 2

−5

ln ε ≈=

−1.13 and 2

−6

ln ε ≈= −0.56. These CAMs correspond to individual HH stretches.

The CAMs on the first and last 2 rows cover the degeneracy between the 3 hydrogens within the respective methyl groups, whereas the third and fourth rows show

isolated HH vibration modes.

The next 2 powers of T (T

128, T256) do not filter out any new spaces. The third

wavelet subspace is spanned by a single vector (λ =72.7 kcal/mol/Å2, log(1 −

λ/660) ≈ −0.12),

 

H C H H C H H C H H C H H H

0 −0.22 0 0 0.59 0.20 0.20 −0.59 −0.20 −0.20 0.22 0 0 0 !

,

(18)

with exponent n = 9 for T

2

n

. The logarithm of its expectation value with T is between 2

−6

ln ε ≈= −0.56 and 2

−9

ln ε ≈= −0.070. This CAM isolates the asymmetric stretch between the terminal carbons and the centers of mass of the bridging

methylene groups.

The fourth wavelet subspace is spanned by another vector (λ =42.1 kcal/mol/Å2,

log(1 − λ/660) ≈ −0.066),

 

H C H H C H H C H H C H H H

0.13 0.43 0.13 0.13 −0.46 −0.15 −0.15 −0.46 −0.15 −0.15 0.43 0.14 0.14 0.14!

,

(19)

with exponent n = 10 for T

2

n

= T

1024. The logarithm of its expectation value with

T is between 2

−9

ln ε ≈= −0.070 and 2

−10 ln ε ≈= −0.035. This CAM captures

the symmetric stretch between the center of mass of the bridging methylenes and

the centers of mass of the terminal methyl groups.

12
Approved for public release; distribution is unlimited.

The last nontrivial wavelet subspace also is spanned by a single vector (λ = 12.2

kcal/mol/Å2, log(1 − λ/660) ≈ −0.019),

 

H C H H C H H C H H C H H H

−0.17 −0.59 −0.17 −0.17 −0.24 0 0 0.24 0 0 0.59 0.17 0.17 0.17!

,

(20)

with exponent n = 11 for T

2

n

= T

2048. The logarithm of its expectation value

with T is between 2

−10 ln ε ≈= −0.035 and 2−11 ln ε ≈= −0.018. This CAM

captures the symmetric stretch of the bridging carbons and the centers of mass of

the terminal methyl groups.

The coarsest level (0 kcal/mol/Å2) is described by

 

H C H H C H H C H H C H H H

0.13 0.45 0.13 0.13 0.45 0.13 0.13 0.45 0.13 0.13 0.45 0.13 0.13 0.13!

.

(21)

This last CAM is just the center of mass.

2.1.3 Multiresolution Construction from Fragments

Since the previously mentioned invariant subspaces are inherent to molecular fragments and some molecular fragments are particularly common, the question of how

much can be gained by precomputing the internal scales of these fragments arises.

For example, proteins are long heterogeneous polymers, but they are mostly composed of only 20 amino acids. Hence, it is instructive to characterize the effects of

modifying bonds or substituting different atoms. Indeed, connecting fragments (see

Section 11 for details) affect a small portion of precomputed CAMs, such that only

a few CAMs need to be adjusted. Precomputed fragments can therefore speed up

the computation of scales considerably when the perturbations due to connecting

them are relatively small.

2.1.4 Example: Linear Homopolymers

Linear homopolymers are an important class of materials, whose graph Laplacians

exhibit convenient structures that we exploit in the following to derive their CAMs.

Linear homopolymers are a successive addition of edges between identical building blocks. We can derive computationally inexpensive algorithms to compute the

eigenvalues of their weighted graph Laplacians and thereby the successive construction of the respective wavelet spaces. The eigensystems are computed by exploiting

the recursive structure of ∆ ̃ to solve n much smaller eigensystems, where n is the

13
Approved for public release; distribution is unlimited.

number of repeat units in a single chain of the polymer. The graph Laplacian of

linear homopolymers can be ordered to have a block Toeplitz structure, where each

nonterminal block is a constant m × m matrix for the off-diagonal B and diagonal

A, respectively,

∆ =  ̃





A B∗ 0 B

B

.

.

.

.

.

. 0

0

.

.

. B∗

B∗ 0 B A





, (22)

where m is the number of particle DoFs in the repeat unit of the polymer. Furthermore, the off-diagonal block B generally consists of a single nonzero entry representing a single bond connecting repeat units, such as B1,m (off-diagonal blocks in

Eq. 15 for an example). This highly repetitive structure can be exploited to compute

eigenvectors and eigenvalues very efficiently and hence the CAMs (Section 12).

2.2 Computational Methods

The following simulation exemplifies the effectiveness of the approach. All MD

simulations were performed in the the Large-scale Atomistic/Molecular Massively

Parallel Simulator (LAMMPS)

29. A crystalline model of polyethylene (PE) consisting of 100,800 atoms in a single infinite chain was simulated in the NVT ensemble

under periodic boundary conditions with 1-fs time steps, PCFF28 at 500 K. The

monoclinic unit cell has edge vectors a = (103.432 , 0, 0), b = (0, 147.87 , 0), and

c = (73.88 , 0, 50.78 ). The initial structure is depicted in Figure 5.

2.3 Simulation Results and Analysis

We witnessed this behavior in computations of PE as well as polyglucose. Furthermore, when applied to the alanine dipeptide model, commonly used for testing

coarse-graining methods for proteins, the wavelet decomposition correctly identified the various chemical groups at the finer scales, while the coarse scales captured

the partitioning corresponding to single bond rotations. Figure 4 illustrates the decomposition for tetra-1,4-D-glucose and Figure 5 demonstrates the informational

content in the wavelet DoFs for PE.

14
Approved for public release; distribution is unlimited.

Fig. 4 Heat map of selected wavelet DoFs in 1,4-D-glucose tetramer. Blue denotes positive

coefficients, red negative. From left to right: fine-grained wavelet DoF covering a single repeat

unit; CG wavelet DoF covering one half of the oligomer; coarser wavelet DoF covering the

oligomer isolating repeat units with sign changes; coarsest wavelet.

a) b)

c) d)

Fig. 5 Selected reduced wavelet representations of a PE crystal. From a) to d) wavelet information is successively added. a) Coarse representation highlights the anisotropy in 1-D chain

averaged over all chain segments; b) Coarse representation of the 2-D chain plane; c) 3-D

representation; d) full resolution.

Figure 6 shows the superimposed Fourier-transformed time series of 1,000 atoms

from the 500 K trajectory, which due to the high temperature shows the highest

mobility of atoms. Although, the zero frequency is by far the most intense signal

(and is omitted from the figure for clarity), a wide range of other frequencies are

active, most importantly around 0.45 PHz, which limits the time step of atomistic

simulations of PE to 1.2 fs or less.

15
Approved for public release; distribution is unlimited.

Fig. 6 Fourier transform of the y-component of 1,000 atoms in crystalline PE (100,800 atoms)

excluding zero frequency to allow detail at other frequencies. MD at 500 K and 1 atm. Left:

Individual power spectra per atom. Right: Power spectrum of magnitude of optimal representation.

On the other hand, Figure 7 shows the effects of scales on time-series analysis. The

finest scale out of 25 still retains the high-frequency components (top, left) as may

be expected, but they are much less intense than the remaining modes. Traversing

the scales, it is noteworthy that a decreasing number of DoFs at the coarser scales

show significant peaks at all. At the medium scale (top, right), no high-frequency

components are found anymore. Therefore, medium scale and coarser DoFs are

quasi-static compared to the finer scales. Furthermore, the signals are clearly separated and very sharp despite the high temperature, which speaks for strongly decoupled modes and justifies dropping the finer scales, which in turn facilitates speed-up

by reducing not only DoFs but also increasing the propagation time step.

Fig. 7 Top row: Fourier transform of the y-component of a 100,800 atom crystalline PE sampled at 1 fs. 3 scales are shown: 1st (left, finest scale), 5th (middle), and 12th scale (out of 25,

right).

The same is witnessed for alanine dipeptide and polyglucose (not shown). Figure

16
Approved for public release; distribution is unlimited.

8 shows the power spectra for all 22 DoFs. Clearly, a large number of frequencies

are active for all atoms. The spectrum of optimal representation shows some clear

peaks around 0.13 and 0.10 PHz. But in Figure 9, the 0.13-PHz peak is absent after

the second scale and the 0.10-PHz peak including noise down to about 0.06 PHz

vanishes from the fourth scale on through the remaining 4 scales.

Fig. 8 Fourier transform of the z-component of alanine dipeptide in vacuum excluding zero

frequency to allow detail at other frequencies. MD at 500 K and 1 atm. Left: Individual power

spectra per atom. Right: Power spectrum of magnitude of optimal representation.

17
Approved for public release; distribution is unlimited.

Fig. 9 Fourier transforms of alanine dipeptide DoF time series in vacuum excluding zero

frequency to allow detail at other frequencies. MD at 500 K and 1 atm. a) z-component of optimal representation; b) second finest optimal wavelet DoF z-component; c) third finest optimal

wavelet DoF z-component; d) fourth finest optimal wavelet DoF z-component.

We investigated the correlation matrices of the various coordinate systems to quantify the extent of interdependence. Figure 10 shows the correlation matrices,

Cij =

R

x

∗

i

(t)xj (s − t) ds dt

R

x

∗

i

(t)xi(s − t) ds dt Rx

∗

j

(t)xj (s − t) ds dt 1

2

, (23)

where xi(t) is a time-dependent coordinate, as heat maps. As may be expected, all

particle coordinates are strongly correlated, but in both the water case as well as

the alanine dipeptide case, much less correlation is witnessed for the wavelet DoFs

(Figure 10, right), despite the few DoFs involved in these small systems. The extent

of overall correlation can be assessed by the `1-norm of the correlation matrices,

about 28 versus about 21 for water, and about 315 versus about 168 for alanine

dipeptide.

18
Approved for public release; distribution is unlimited.

Fig. 10 Top: Heat maps of the DoF correlations of a water dimer in vacuum. Bottom: Heat

maps of the DoF correlations of alanine dipeptide in vacuum. Left: Correlations of the cartesian DoFs. Right: Correlations of the wavelet DoFs.

A further indication of the appropriateness of the diffusion wavelet DoFs, is the

`1-norm of structures and forces in a given representation basis. The `1-norm is

bounded by the `2-norm for any orthonormal basis and there exists at least one

coordinate system in which kxk1 = kxk2. These coordinate systems optimally describe a state x.

We compared the `1-norm in the Cartesian coordinate system, the nonlossy wavelet

coordinate system, and the optimal coordinate system on each scale. For alanine

dipeptide, a reduction by a factor of about 2 was achieved by switching to wavelets

(3 times the optimum), while the optimal wavelet system managed to reduce to

twice the optimum. For polyglucose, the full wavelet DoFs reduce the `1-norm to

half of the Cartesian coordinate system, and the optimal wavelet system reduces it to

within 20% of the optimum. For crystalline PE, the full wavelet DoFs improve the

structural `1-norm by a factor of about 7, while the optimal wavelet representation

comes in below 1% of the Cartesian coordinate system to a factor of roughly twice

the optimal representation. Similar results were obtained for amorphous PE and

19
Approved for public release; distribution is unlimited.

nanocellulose.

2.4 Adaptive Multiresolution

In the previous sections, our goal was to identify scales of the physical system to

remove unimportant CAMs, thereby increasing time step size and reducing DoFs.

However, for some phenomena, such as phase transitions like melting of a polymer,

fine details that are unimportant at one point in time can play a major role at another.

So we now turn to the problem of reconstruction, that is, reintroducing dropped

CAMs. We find that reconstruction is systematically possible for numerical (e.g.,

derived from an iterative Boltzmann inversion) as well as analytical (where such

exists, e.g., quadratic potentials) coarse-graining hierarchies.

2.4.1 Reconstruction Theory

We begin by putting coarse-graining into a wider context. In general, a coarsening

γ : α → β is a continuous surjection between 2 topological state spaces α, the

fine space, and β, the coarse space, which can be parameterized by n > m state

variables, respectively. As an example, one can consider mapping the positions of

groups of particles to their centers of mass. In statistical mechanics, the fine-grained

state space α is associated with a probability measure Pα : {X ⊂ α} → [0, 1] and

a corresponding probability space. For systems in equilibrium, α would further follow a Boltzmann distribution. The coarsening γ thereby induces a probability space

on the coarse states in β as well, with the probability measure Pβ(k) = Pα(γ

−1

(k)),

where k ⊂ β and γ

−1

(k) ⊂ α is the preimage of k. Pβ constitutes the mean probability distribution of the coarsened DoFs. It is thus possible to select (reconstruct) a

precursor for a state b ∈ β by sampling γ

−1

({b}) with Pα via the conditional probability P(a|b) = Pα({a} ∩ γ

−1

(b))/Pα(γ

−1

(b)) (e.g., using Monte-Carlo sampling).

Since n > m, there is also a complementary coarsening γ

⊥ : α → ker γ with its

complementary probability measure.

In MD, the state spaces consist of the positions and their associated momenta and

thermodynamic state variables, such as temperature or pressure. The probability

distributions are Boltzmann distributions that depend on the studied thermodynamical ensemble. In a sequence of coarsenings (γn : βn → βn+1)n spanning several

scales, such as derived from the preceding multiresolution scheme, it is generally

not cost effective to sample fully in the largest space α =: β0 and analytical derivations for Pβn

are rarely available. Hence, approximations need to be made. Com20
Approved for public release; distribution is unlimited.

mon solutions in the MD community are probability measures from iterative Boltzmann inversions or (successive) force matching to generate effective potentials that

are themselves Boltzmann-distributed. Hierarchical iteration thereby produces not

only probability distributions on the coarser space, but also conditional probability

distributions for a fine state mapping to a coarse state. Furthermore, the probability distributions can be used to indicate when a previously undersampled coarse

state subspace is encountered at some state x (e.g., using an expected improvement

measure30 or a sensitivity analysis of the potential − ln Pβ(x) with respect to the

sampled points for which the trust boundaries can be precomputed).

2.4.2 Coarse-Graining Hierarchies

In the context of MD, mixed resolutions are necessary for the concurrent treatment

of, for instance, gross rigid protein orientations and flexible active site residue interactions with a binding molecule. State spaces with mixed levels of fine and coarse

CAMs arise naturally from the multiresolution scheme laid out in Eqs. 9 and 10

(see Section 13 and Figure 11 for details).

βn

βAA⊥

β

⊥

A

μ

⊥

A

βA

μA

μA ⊕ μ

⊥

A

AB

μA ⊕ μB

βBB⊥

βB

μB

β

⊥

B

μ

⊥

B

μB ⊕ μ

⊥

B

AB⊥

Fig. 11 Dyadic tree generated by separable coarsenings. Shaded spaces are final subspaces of

the multiresolution. All other spaces are intermediates of the construction with γ1 = μA ⊕ μ

⊥

A

and γ2 = μA ⊕ μB.

The sequence of coarsenings (γn : Vn → Vn+1)n generates a full dyadic tree, due

to the complementary coarsenings γ

⊥

n

(i.e., each state space can be split into ker γ

⊥

n

and ker γn). A sequence of separable coarsenings (γn, βn, μAn, μBn) thereby induces a hierarchy of coarsenings and associated probability measures, induced as

described previously.

21
Approved for public release; distribution is unlimited.

2.4.3 The Wavelet Hierarchy

A minimal set of separable coarsenings that generates a full given hierarchy is

fundamental and characterizes said hierarchy. One such sequence for coarsenings

based on CAMs for Vm = (Lb2−mNc

n=1 W2mn−1)

log2 N

m=0 is

μAm :

Lb2−m−1Nc

n=1 W2mn−1 →

Lb2−m−2Nc

n=1 W2m+1n−1, (24)

μBm =

Lb2−mNc

n=b2−m−1Nc+1 W2mn−1 →

Lb2−m−1Nc

n=b2−m−2Nc+1 W2m+1n−1, (25)

where N is the number of particles in the finest resolution.

We note that these intermediate probability distributions are available both analytically and numerically, since an accurate fundamental coarsening has to include

proper statistics for the intermediate state space to be consistently sampled. Recursive application of conditional probabilities enables concurrent mixed resolutions.

Since the construction of CAMs from the multiresolution analysis produces a hierarchy of frequencies, it induces a hierarchy of coarsenings by dropping successively

higher-frequency CAMs (i.e., by applying the low-pass filter Pε(T

2

n

)).

2.4.4 A Priori Approximations for Reconstruction

Implementation of reconstruction algorithms as discussed above requires a starting

point. In the following, methods are proposed for finding good starting points based

on ∆ ̃ and other molecular information that is available prior to simulation (i.e.,

without the need for MD or Monte Carlo sampling). To second order, a quadratic

potential around the equilibrium positions of the transformed coordinates approximates the full potential. We assume dominance of harmonic terms, both in the

orginal and transformed basis of the potential,

Vbond(M− 1

2U

T

r ̃) ≈

1

2

( ̃r − r ̃0)

TU∆ ̃UT

( ̃r − r ̃0), (26)

where r ̃ = UM1/2x is the position vector in the wavelet basis, and UM1/2is the

wavelet transformation matrix. From statistical mechanics, the root mean square

deviation from equilibrium of a harmonic oscillator is p

kT/λ, where λ is the force

constant. In other words, the higher frequency components are found increasingly

close to their energy minima. This implies that finer scales only have small deviations from their equilibrium positions, while coarser scales may access a much

larger space.

22
Approved for public release; distribution is unlimited.

We start by approximating kxi − xjk by a Taylor expansion around r

(0)

ij . This transforms the bond potential Vbond into

Vbond ≈

1

8

X

ij

Kij

r

(0)2

ij



kxi − xjk

2 − r

(0)2

ij 2

. (27)

Hence to find equilibrium distances for r ̃ = UM1/2x, we solve the minimization

problem

min

r ̃0

XCij 

k(M− 1

2U

T

r ̃0)i − (M− 1

2U

T

r ̃0)jk

2 − r

(0)2

ij 2

, (28)

where Cij = Kij/r(0)2

ij .

For example, the 2 nonzero eigenvalues of H2O correspond to a unique solution for

reconstructing H2O. The harmonic Laplacian for H2O,

∆ ̃ H2O =





2KOH

mO

−

KOH

m

1

2

O

−

KOH

m

1

2

O

−

KOH

m

1

2

O

KOH 0

−

KOH

m

1

2

O

0 KOH





, (29)

shares the same structure as CH2. The eigenvalues λ0,1,2 of this simple matrix

are 0, KOH, and (1 + 2/mO)KOH, respectively, with corresponding eigenvectors

(m

1/2

O

, 1, 1), (0, 1, −1), and (2m

−1/2

O

, −1, −1). The transformation matrix from r ̃ to

x is

W =





m

− 1

2

O

1

1









m

1

2

O

0 2m

− 1

2

O

1 1 −1

1 −1 −1









(mO + 2)− 1

2

2

− 1

2

(2 + 4/mO)

− 1

2



 .

(30)

The corresponding harmonic potential for water is

VH2O(x0, x1, x2) = 1

2

KOH (kx0 − x1k − rOH)

2 + KOH (kx0 − x2k − rOH)2

.

(31)

23
Approved for public release; distribution is unlimited.

This attains its minimum when

kx0 − x1k

2 =kr ̃2k2



1

2

+ m−1

O



−

q

1 + 2m−1

O

r ̃

T

1

r ̃2 +

1

2

kr ̃1k

2 = r2

OH, (32)

kx0 − x2k

2 =kr ̃2k2



1

2

+ m−1

O



+

q

1 + 2m−1

O

r ̃

T

1

r ̃2 +

1

2

kr ̃1k

2 = r2

OH, (33)

where we have used that x0−x1 =  ̃r2

p

1/2 + 1/mO−r ̃1/

√

2, x0−x2 =  ̃r2

p

1/2 + 1/mO+

r ̃1/

√

2. Hence, r ̃1 and r ̃2 must be perpendicular to each other (subtracting Eqs. 32

and 33) and kr ̃2k

2

(1/2 + 1/mO) + 1

2

kr ̃1k

2 = r2

OH. As these solutions are exact, the

hydrogens are always at a distance of rOH from the oxygen. Since there is no angular potential, the solution is indeterminate with 2 extreme solutions, the first being

kr ̃2k = 0. In this case, the molecule is linear, while for kr ̃1k = 0 the 2 hydrogens

overlap.

Another instructive example is HCN. Its harmonic Laplacian,

∆ ̃ HCN =





KCH+KCN

mC

− √

KCH

mCmH

− √

KCN

mCmN

− √

KCH

mCmH

KCH

mH

0

− √

KCN

mCmN

0

KCN

mN



 , (34)

is no longer as simple as for H2O, nor are the eigenvalues except 0 simple functions

of the variables in Eq. 34; for the generalized amber force field (GAFF), they are

14.7 and 42.8. Eq. 28 was numerically solved. The numerical GAFF CAM distances

in one dimension are ( ̃r

(0)

1

, r ̃

(0)

2

) = (3.26, 0.66) and r ̃

(0)

1 = 2.61, r ̃

(0)

2 = −1.32.

Knowing that the molecule is linear at equilibrium selects the first solution to reconstruct the equilibrium structure.

In both examples it was necessary to include angle information to make the best

choice. The numerical solution to Eq. 28 can be computed efficiently using a variety of nonlinear least-squares algorithms, but more direct methods are still under

investigation. Similar derivations are possible for angle potentials and under current

investigation.

24
Approved for public release; distribution is unlimited.

3. Operator Wavelets Theory

To perform efficient wavelet representations of DoFs, there also needs to be an

appropriate effective operator that can exploit the wavelet structure of the CAMs.

Ideally, this operator can itself be constructed in a hierarchical fashion. The following explores the properties of operator hierarchies that arise from doubling and shift

operators. Repeated application of an invertible doubling operator D : L1 → L1

and an invertible shift operator σ : L1 → L1 on an operator space L1 generates a

scaling operator Φ analogous to a scaling function in wavelet theory on functions.

Then,

D(I + σ)Φ = Φ. (35)

Since the identity is subject to D and Dσ, there are associated operators d ∈ L1

and ds ∈ L1, respectively. If furthermore D and σ are endomorpisms on L1, then

(d + ds)Φα = (DI + DσI)(DΦ + DσΦ)α = Φα ⇒ (36)

(d + ds)α = α (37)

for any eigenfunction α of Φ. Hence the eigenfunctions of Φ are scaling functions

under the doubling operator d and shift operator s. The corresponding operator

wavelets are

Ψnk = D

nσkD(I − σ)Φ, (38)

ψα,n,k = d

n

s

k

(d − ds)α. (39)

An approximate procedure to develop a hierarchy that is not based on endomorphism on L1 is to use the doubling operator Φ 7→ 2Φ and use polynomial filtering

on the spectrum (i.e., Φ 7→ F(Φ) = P

i FiΦ

i where the polynomial F approximates a low-pass filter). Under this construction any scaling operator is unitarily

equivalent to any other (i.e., they share the same spectrum).

4. Trotter Factorization for Use with Operator Wavelets

A hierarchical decomposition of the infinitesimal generator of time propagation (Liouville operator) into operator wavelets enables the separation of time scales rigorously using Trotter factorizations for numerical propagation. We split the Liouville

25
Approved for public release; distribution is unlimited.

operator into

H =

X

N

i=0

Hi +

X

i,j∈{0,...,N}

Hij , (40)

where Hi are individually solved subsystems and Hij are the interactions between

these, such as the operator wavelets introduced in Section 3. Preferably, these partitions coincide with inherent time scales ti and tij . Let the time scales be in ascending order and readjusted to the closest integer multiple (i.e., t ̃I+1 =



tI+1/t ̃I



t ̃I )

and τ = max t ̃I . If t ̃I+1 = t ̃I , it is useless to separate HI+1 and HI . So then,

e

iH τ ≈

Y

N

I=0



e

iHI t ̃I /2

τ /t ̃I Y

0

I=N



e

iHI t ̃I /2

τ /t ̃I

. (41)

When the time propagator e

iHI t ̃I has analytical solutions, such as with harmonic oscillators, then the analytical solution e

iHI τ /2

can be used instead of the τ/t ̃I applications of the approximation. If τ/t ̃I is large, considerable savings may be expected.

In particular, consider a quadratic potential with matrix ∆ ̃ (i.e., iHq = −(x −

x)

T ∆ ̃ d

dp + p

TM−1 d

dx ), then

e

iH τ ≈ ei(H −Hq)τ /2

e

iHqτ

e

i(H −Hq)τ /2

. (42)

The middle term,

 

x(τ )

p(τ )

!

= e

iHqτ

 

x

p

!

= e





I

−∆ ̃



τ 

x(0) − x

p(0) !

+

 

x

0

!

, (43)

can be solved analytically, so any time interval τ can be computed exactly, whereas

H − Hq =



−∇T V + (x − x)

T ∆ ̃

 d

idp (44)

requires only a single force evaluation similar to other reversible reference system

propagator algorithms.

To leverage this factorization, an appropriate origin x needs to be assessed either

by running short bootstrapping trajectories or analytical expressions specific to the

system. For high frequencies, the propagation by Hq may wildly oscillate through

26
Approved for public release; distribution is unlimited.

the origin. Which origin is optimal may actually depend on the amount of energy

in the quadratic term.

5. Solutions to "Harmonic" Oscillators

The Trotter factorization in Section 4 is most useful if analytical solutions exist for

the involved operators. The Liouville operator L for bond oscillators in MD is

iL =

X

ij

−kij 

kxi − xjk − r

(0)

ij  xi − xj

kxi − xjk

∂

∂pi

+

X

i

pi

mi

∂

∂xi

. (45)

We set η = Ex, where E is the incidence matrix (i.e., η are the edges of the graph

spanned by the bonds). Furthermore, let E

T

θ = p. Then

∇x = E

T ∇η, (46)

E∇p = ∇θ. (47)

Inserting these into Eq. 45, gives

iL =

X

i

− ̃kI (kηIk − r

(0)

I

)

ηI

kηIk

∇θI + θ

TEM−1ET ∇η. (48)

Assuming a simple ansatz that

−  ̃kI (kηIk − r

(0)

I

)

ηI

kηIk

∇θIfI + θ

T

I EIjM−1

jj EIj∇ηI

fI = λIfI , (49)

then

iL g(η, θ) = θ

T

(EM−1E

T − D)∇ηg + Λg, (50)

where g =

PcIfI (ηI , θI ) and Λg =PcIλIfI (ηI , θI ).

The Liouville operator L is rotationally invariant (i.e., L commutes with any rotation R(φ)). Hence, there are simultaneous eigenfunctions of L and η

TRˆ∇η +

θ

TRˆ∇θ, where Rˆ =

 

0 1

−1 0!

. Eigenfunctions are compositions of (1 ± iRˆ)η and

(1 ± Rˆ)θ with eigenvalues ±i.

The following equations are some helpful identities that were used in the previous

27
Approved for public release; distribution is unlimited.

derivations:

η = ηr

 

cos φ sin φ

− sin φ cos φ

!

v (51)

θ = θr

 

cos φ − sin φ

sin φ cos φ

!

v (52)

∇v = ηrR(φ)

T ∇η + θrR(φ)∇θ (53)

d

dηr

= v

TR(φ)T ∇η (54)

d

dθr

= v

TR(φ)∇θ (55)

d

dφ = v

T

dR

dφ (φ)

T ∇η + vT

dR

dφ (φ)∇θ (56)

−k(ηr − r0)v

TRT

(φ)∇θ + θrv

TR(φ)M∇η (57)

6. Hermite Functions for Approximating Potentials

In the decomposition of the Liouville operator in Sections 3 and 4, a suitable basis

of potential functions is needed. This basis is preferably sufficiently localized to

avoid the need for global integrals. The Hermite functions,

ψn =(−1)n



2

nn!

√

π

− 1

2 e

1

2

x

2 d

n

dx

n

e

−x

2

=(−1)n



2

nn!

√

π

− 1

2 Hn(x)e

− 1

2

x

2

, (58)

ψ2n =(−1)2n



2

2n

2n!

√

π

− 1

2 e

1

2

x

2 d

2n

dx

2n

e

−x

2

=(−1)2n



2

2n

2n!

√

π

− 1

2 H2n(x)e

− 1

2

x

2

, (59)

are an orthonormal complete set of basis functions that obey a 3-term recursion

formula,

ψn+1(x) = r

2

n + 1

xψn(x) −

r n

n + 1

ψn−1(x), (60)

ψ2n+2 =

(2x

2 − 4n − 1)

p

(2n + 2)(2n + 1)

ψ2n −

s

2n(2n − 1)

(2n + 1)(2n + 2)ψ2n−2. (61)

The recursion formula enables a fast transform of functions. Since "harmonic" po28
Approved for public release; distribution is unlimited.

tentials are merely functions that capture 3 aspects of bonding, it is possible to use

Hermite functions in their stead. The boundary conditions are

f(0) = r

2

0 − D, (62)

df

dx

(r0) = 0, (63)

f(r0) = −D, (64)

d

2

f

dx

2

= k (65)

f(x) = f(−x), (66)

where D is the dissociation energy, k is the force constant, and r0 is the equilibrium

distance. These conditions result in a linear system for coefficients in an expansion

of Hermite functions,

f =



α + β

x

σ

2

+ γ

x

σ

4



e

− 1

2 (

x

σ )

2

. (67)

A simple decomposition converts this equation into the coefficients of the Hermite

functions of scale σ! The Liouville operator for such a potential would be

iL = − ∇T

i

f(kxi − xjk)

d

dpi

− ∇T

j

f(kxi − xjk)

d

dpj

+ m−1

i

p

T

i

d

dxi

+ m−1

j

p

T

j

d

dxj

=

xi − xj

σ

 

f +

 

2β + 4γ













xi − xj

σ













2

!

e

− 1

2










xi−xj

σ










2

! 

d

dpi

−

d

dpj



+

m−1

i

p

T

i

d

dxi

+ m−1

j

p

T

j

d

dxj

. (68)

A further point of note is the fact that this expansion only includes quadratic terms

with respect to the DoFs, which allows much simpler decomposition into component terms after linear transformation of the DoFs into CAMs. Finally, geometries

may be defined by distances obviating the need for angular and dihedral potentials.

29
Approved for public release; distribution is unlimited.

7. Koopman for MD

In the following, the relationship between operator wavelets and the Koopman

modes of the Liouville operator is demonstrated. The Liouville operator L is a

Koopman operator for MD systems. Let xˆ(x, p) = P

{λ}

xλKλ(x, p), where Kλ is

an eigenfunction of L and xˆ : (x, p) 7→ x. Then,

x(t) = e

−iL txˆ(x0, p0) = Xeiλj txλjKλj

(x0, p0), (69)

where λjis the corresponding eigenvalue to Kλj

. A decomposition of L using

an operator hierarchy as in Section 3 naturally separates the scales by the wavelet

operators. Therefore, Kλj

is in the kernel of all nontrivial operator wavelets save

one. The eigenvalues λj are generally degenerate. Grouping eigenfunctions Kλ,k

with the same eigenvalue

x(t) = Xe

iλj Xxj,kKj,k(x0, p0) = XAjeiλj

, (70)

reduces the number of actually free variables, where Aj are the effective free variables. Assume 2 observables evolve in time

x(t) = e

−iL txˆ(x0, p0) = Xeiλj txλjKλj

(x0, p0), (71)

y(t) = e

−iL t

yˆ(y0, q0) = Xe

iλj t

yλjKλj(y0, q0). (72)

Then, cross-correlation indicates

hx|yi =

Xx

∗

λj

yλjKλj(x0, p0)Kλj(y0, q0). (73)

It follows that the eigenvectors of the correlation matrix of a set of observables is

the optimal linear representation of same observables. This includes the original

DoFs as well as CAMs.

We can further characterize the eigenfunctions Kλ of L for physical systems by its

distributive property,

L (f · g) = g · L f + f · L g. (74)

30
Approved for public release; distribution is unlimited.

Hence, any power of an eigenfunction of L is itself an eigenfunction of L via

L Kλ = λKλ ⇒ L Kn

λ =

Xn−1

m=0

Km

λ

(L Kλ)Kn−m−1

λ = nλKn

λ

. (75)

This motivates a simple shift operator K−1

λ

and a doubling operator D : f(Kλ) 7→

f(K2

λ

), where f is a formal series in polynomial powers of Kλ.

The effect of the Koopman operator can be captured by bootstrapping. With an

initial trajectory, the effective frequencies can be found by solving

min

{x ̃j},{wj}

X

N

n=0
















x(n∆t) −

X

M

j

e

iwjn∆tx ̃j
















2

+ C

Xkx ̃jk, (76)

where N is the number of time steps in the initial run, ∆t is the time step, M is the

number of frequencies to capture, and C is a cost factor to control overfitting with

regularizer Pkx ̃jk. By comparison with Eq. 71, x ̃j ≈ xwj. Then ∆x := x(t) − PM

j

e

iwj tx ̃j

is the residual due to ∆L := L − L ̃, where L ̃ is the linear operator

with eigenfunctions and eigenvalues (Kwj, wj ) such that f 7→ 0∀f 6∈ span {Kwj}.

It follows that L ̃ and ∆L are orthogonal operators and commute. Thus

e

iL t = ei∆L t+iL ̃t = ei∆L t

e

iL ̃t = eiL ̃t

e

i∆L t

. (77)

At time t, the action of ∆L on x is ∆p := p(t)−

Piwjx ̃je

iwj t

, whereas the action

on ∆p is −∇V (x(t)) + Pw

2

jx ̃je

iwj t

.

Using the usual Velocity-Verlet factorization going from time t = 0 to t = ∆t,

∆x1 = ∆p0∆t/2 + ∆x0, (78)

∆p2 = ∆p0 −



∇V (∆x1 +  ̃x(∆t/2)) − p ̃ ̇(∆t/2)t, (79)

∆x3 = ∆x1 + ∆p2t/2 = ∆x0 + (∆p0 + ∆p2)t/2, (80)

x(∆t) = ∆x3 +  ̃x(∆t), (81)

p(∆t) = ∆p2 +  ̃p(∆t). (82)

This approach is equivalent to introducing new variables αj and βj:=  ̇αj and

31
Approved for public release; distribution is unlimited.

changing the Liouville operator to

iL = tr h(∇T

x V )



x +

Xαj



∇p

i

+ tr pTM−1∇x +

Xβ

T

j ∇αj − ω

2

jα

T

i ∇βj

.

(83)

8. Derivation of Wavelet Spaces

The wavelets are derived iteratively from a filter T using the QR decomposition. At

each iteration, T is cast in the “Q” basis of the previous iteration,

 Yn

i=0

Qi

!T

T

2

n+1 Yn

i=0

Qi

!

= Qn+1Rn+1. (84)

Therefore, repeated application (with infinite precision) is equivalent to the QR algorithm for finding eigenvectors. Since T is positive definite and kTk∞ = 1, the

squaring introduces a de facto projection operator Pε via the machine precision. It

is this projection that distinguishes a conventional QR algorithm for finding eigenvalues and eigenvectors from the wavelet decomposition into CAMs.

We separate Qn into a low-frequency submatrix Φn and a high-frequency submatrix

Ψn, where the latter are the columns Qn,i of Qn for which QT

n,iT

2

n

Qn,i < ε. The

matrix Φn collects the remaining columns of Qn. Thus,

 nY−1

i=0

Qi

!T

T

2

n+1 nY−1

i=0

Qi

!

≈ Φ

T

nQnRnR

T

nQ

T

nΦ

T

n ≈ (

Yn

i=0

Φ

T

i

)T

2

n+1 (

Yn

i=0

Φi).

9. Error Bounds on Scales

The contamination of the wavelet spaces Wn by eigenvectors U> of larger eigenvalues λ

2

n

> ε of the filter operator T is limited by

S

(n)

T

:=

ε − ωn,max

ωn,min − ε

, (85)

where ωn,min = min{λ

2

n

> ε|λ ∈ σ(T)}, ωn,max = max{λ

2

n

≤ ε|λ ∈ σ(T)},

and σ(T) denotes the spectrum of T. In classical wavelet theory, wavelets are localized both in real and reciprocal space (e.g., the Fourier transform of the Mexi32
Approved for public release; distribution is unlimited.

can hat wavelet transform W f(x) = −σ

−2∆

R

e

−(y−x)

2/(2σ2

)

f(y)dy is W ̃  ̃f(ω) =

ω

2σ−2

e

−ω

2/(4σ2

)  ̃f(ω)). The sensitivity S

(n)

T

of T at scale n measures the localization

of Wn in the frequency domain. In practical terms, small S

(n)

T

implies that the bases

Φn and Ψn are generally sparse matrices if T

2

n

is sparse, and the wavelet transform

can thus be computed efficiently.

10. Invariant Subspaces

∆ ̃ is positive semi-definite and ∆ ̃

JJM1/21J = 0 for each set of indices J, where

∆ ̃

JJ is the square submatrix of ∆ ̃ with indices in J, and 1J is the vector of ones

on indices in J and 0 otherwise. It is possible to block tridiagonalize ∆ ̃ using only

transpositions with diagonal blocks ∆ ̃

J,J and off-diagonal blocks ∆ ̃J,K, where J

and K are disjoint index sets and J, K, J ∪ K are contiguous index sets. Without

loss of generality, we let j < k∀j ∈ J, k ∈ K. The rank of ∆ ̃

JK is generally

low due to the low maximum degree of a vertex in ∆ ̃ . If Λ = span ∆ ̃ KJ and ∆ ̃

JJ

has a nontrivial, invariant vector space Γ perpendicular to Λ, then Γ is localized

to indices preceding K. Examples include linear homopolymers discussed in detail in Section 2.1.4 and Section 12, but also disconnected graphs from individual

molecules. Invariant subspaces {Γ} are computationally convenient because they

allow a separation of the problem into independent smaller subproblems.

11. Perturbation Theory for Molecular Fragments

Changing the mass of atom i by δ leads to M0−1/2 = M−1/2 + δeie

T

i

, and similarly

∆ ̃ becomes

∆ ̃ 0 =(M− 1

2 + δeie

T

i

)∆(M− 1

2 + δeie

T

i

) (86)

=∆ + ̃ δm

1

2

i

(eie

T

i ∆ + ̃ ∆ ̃ eie

T

i

) + δ

2mi∆ ̃

iieie

T

i

. (87)

For an eigenvector ν of ∆ ̃ contained in some wavelet space Wn with nondegenerate eigenvalue ην, the first-order correction to the eigenvalue is m

1/2

i

δ|νi|

2

(2ην +

m

1/2

i ∆ ̃

iiδ), and the first-order correction to the eigenvector is P

μ6=ν μ·δμiνim

1/2

i

(ην+

ημ + δmi∆ ̃

ii)/(2ημ − 2ην), where {μ} is an orthonormalized set of eigenvectors of

∆ ̃ . Hence the scale of ν is unaffected if |νi|

2

is 0 or sufficiently negligible. Corrections to the eigenvector are only of significance if some eigenvector μ outside of

Wn additionally has a significant component μi.

33
Approved for public release; distribution is unlimited.

By a similar argument, first-order corrections can be computed for a perturbation of

the edge weights (e.g., changing bond spring constants),

∆ ̃ 0 =∆ + ̃ δM− 1

2 (eie

T

i + eje

T

j − eie

T

j − eje

T

i

)M− 1

2 , (88)

η

(1)

ν =δ|m

− 1

2

i

νi − m

− 1

2

j

νj|

2

, (89)

ν

(1) =

X

μ6=ν

δ

(m

− 1

2

i μi − m

− 1

2

j μj )(m

− 1

2

i

νi − m

− 1

2

j

νj )

ημ − ην

μ. (90)

Hence, a wavelet subspace Wn only changes if ∃ν ∈ Wn : kT νk2 + δkνk∞ ≥

ε

1/2

n

or kT νk2 − δkνk∞ ≤ ε

1/2

n−1

.

12. Derivation of Homopolymer Wavelets

The computation of wavelet spaces for linear homopolymers with n repeat units can

be subdivided into n subproblems. The recursive structure of linear homopolymers

implies that ∆ ̃ for a linear homopolymer of n repeat units can be reordered by a

permutation κ of indices such that κ

T ∆ ̃ κ = A ⊗ In + B ⊗ (ΣT

n + ene

T

1

) + BT ⊗

(Σn + e1e

T

n

), where In is the n × n identity matrix, and Σn is the n × n-matrix with

all ones on the first subdiagonal only and zeros elsewhere. Then, κ

T ∆ ̃ κ takes the

simple form





A1,1In A1,2In · · · A1,mIn + B1,m(Σn + e1e

T

n

)

A2,1In A2,2In · · · A1,mIn

.

.

.

.

.

.

.

.

.

Am,1In + B1,m(ΣT

n + ene

T

1

) · · · · · · Am,mIn





.

(91)

Let U = (e

i2πkl/m/

√

m)kl be a unitary matrix that diagonalizes A1,mIn+B1,m(ΣT

n +

ene

T

1

), then U := Im ⊗ U transforms κ

T ∆ ̃ κ into





A1,1In A1,2In · · · D∗

A2,1In A2,2In · · · A1,mIn

.

.

.

.

.

.

.

.

.

.

.

.

D Am,2In · · · Am,mIn





= U

∗κT ∆ ̃ κU, (92)

where D is a complex diagonal matrix with eigenvalues Dk,k = A1,m +B1,me

ikπ/m.

34
Approved for public release; distribution is unlimited.

Applying the similarity transform by κ

T produces a block diagonal matrix (κU∗κT

∆ ̃ κUκ

T

) with square blocks A ̃

i = A − (A1,m + Di,i)e1e

T

m − (A1,m + D∗i,i)eme

T

1

.

Hence, the computation of the wavelet spaces has been reduced to n problems of

size m instead of one problem of size nm. Furthermore, any eigenvector ν of A

with e

T

mν = e

T

1

ν = 0 has an n-fold degenerate eigenvalue.

13. Derivation of Mixed Resolution from Separable Coarsenings

Let β

τ

n

: R

mn → βn denote a parameterization of the mn-dimensional state space

βn (e.g., the positions and momenta of particles). A coarsening γn can be separated

into components if there exist parameterizations β

τ

n

: An ⊕ Bn → βn and β

τ

n+1 :

An+1⊕Bn+1 → βn+1, such that An⊕Bn = R

mn , An+1⊕Bn+1 = Rmn+1 , and there

exist continuous surjective mappings μAn: An → An+1 and μBn: Bn → Bn+1

with

β

τ

n

(μ

−1

An

(tA), μ−1

Bn

(tB)) = γ

−1

n

(β

τ

n+1(tA, tB)) (93)

for tA ∈ An+1, tB ∈ Bn+1. μA and μB map a fine parametrization to a coarse

parameterization. An example parametrization is the representation of Vn as Vn+1⊕

Wn+1 with fine-to-coarse mappings as per Eqs. 9 and 10. See Figure 11.

A separable coarsening for which neither μAn nor μBnare bijections induces intermediate coarsenings. The state space An × Bn+1 is an intermediate state space

with the coarsenings γ : β

τ

n

(tAn, tBn) 7→ (tAn, μBn(tBn)) ∈ An × Bn+1 and

γ

0

: (tAn, tBn+1 ) 7→ β

τ

n+1(μAn

(tAn), tBn+1 ) ∈ βn+1. Finally, mixed resolution

spaces An × Bn+2 are induced via function composition, γ : β

τ

(tAn, tBn) 7→

(tAn, μBn+1 ◦ μBn(tBn).

The reconstruction from a fundamental coarsening can be achieved by reconstructing its separable components separately and independently via the conditional probabilities P(tAn

|tAn+1 , tBn+1 ) = Pβn(({tAn }∩μ

−1

An

(tAn+1 ))×μ

−1

Bn

(tBn+1 ))/Pβn+1 (tAn+1 , tBn+1 ).

The same can be achieved for the complements analogously.

35
Approved for public release; distribution is unlimited.

14. Conclusions and Outlook

We have characterized a coarse-graining procedure for accelerating molecular simulations through a systematic, hierarchical algorithm based on multiresolution diffusion wavelets. The proposed wavelet-CG approach goes beyond conventional approaches based on expert knowledge: because our proposed method can accelerate calculations of novel classes of molecules without requiring extensive expert

insight and model parameterization. This advantage is especially important for inverse problems in materials design, wherein the materials engineer aims to optimize

material performance in an essentially infinite design space (the chemical space of

polymer repeat units). Our demonstration of the perturbation theory for chemical

variations in the repeat unit illustrates this key advantage for the wavelet CG approach.

Importantly, these advantages are obtained in a framework that automatically recapitulates the physical insights underlying existing coarse-graining methods such

as united-atom models. On the other hand, diffusion-wavelet CG models are simultaneously more general (they do not require a priori expert modeling and parameterization) and more specific. In fact, the diffusion-wavelet CG approach leads

to system-specific CG models derived automatically from the system’s underlying

bonding topology and atomistic force field, without further input other than an error tolerance. The systematic and purely algorithmic basis offers the opportunity for

adaptive error control, whose obvious importance has motivated significant analysis

already31–33.

Future work will establish the relationship between time steps and simulation accuracy34. Currently, the bootstrapping procedure for optimal Koopman mode approximation requires embedding in the propagation scheme and linking to stochastic dynamics. The operator wavelet approach has implications for developing new operator decompositions, which require derivation and solution of analytical subsystems.

Meanwhile, a proper effective operator decomposition using Hermite functions or

conventional harmonic, angle, and dihedral functions as a basis for fitting can be

achieved.

36
Approved for public release; distribution is unlimited.

15. Publications and Presentations

The following are previously published presentations and publications that document work in this effort:

1. Rinderspacher BC, Bardhan J, Ismail AE. Diffusion wavelet-based decompositions for coarse-graining of polymer chains. Paper presented at: American Institute of Chemical Engineers. 2015 Nov10; Salt Lake City, UT.

2. Rinderspacher BC, Bardhan JP, Ismail AE. Diffusion wavelet decomposition

for coarse-graining of polymer chains. MRS Proceedings. 2015;1753.

3. Rinderspacher BC, Bardhan J, Ismail AE. Multiresolution of molecular dynamics in the particle domain. Paper presented at: SIAM Mathematical Aspects of

Materials Science Conference. 2016 May 8; Philadelphia, PA.

4. Rinderspacher BC, Bardhan JP, Ismail AE. Diffusion wavelet decomposition

for coarse-graining of polymer chains. Paper presented at: Materials Research

Society Fall Meeting. 2014 Dec 1; Boston, MA.

5. Rinderspacher BC, Bardhan J, Ismail AE. Diffusion wavelet-based decompositions for coarse-graining of polymer chains. Paper presented at: University of

Delaware, Applied Mathematics. 2015 Oct 5; Newark, DE.

6. Rinderspacher BC, Bardhan J, Ismail AE. Wavelet analysis for molecular dynamics. Paper presented at: MEDE. 2014 Oct 16; Towson, MD.

7. Rinderspacher BC, Bardhan J, Ismail AE. Diffusion wavelet decomposition for

coarse-graining of polymer chains. Paper presented at: American Chemical Society National Meeting. 2015 Aug 16; Boston, MA.

8. Rinderspacher BC, Bardhan J, Ismail AE. Diffusion wavelet decomposition for

coarse-graining of polymer chains. Paper presented at: Hopkins Extreme Materials Institute. 2015 Jul 30; Baltimore, MD.

9. Rinderspacher BC, Bardhan JP, Ismail AE. Theory of wavelet-based coarsegraining hierarchies for molecular dynamics. Physical Review E. ; submitted.

37
Approved for public release; distribution is unlimited.

16. Bibliography

1. Akkermans RLC, Briels WJ. Coarse-grained interactions in polymer melts: a

variational approach. J Chem Phys. 2001;115:6210-6219.

2. Baschnagel J, Binder K, Doruker P, Gusev AA, Hahn O, Kremer K, Mattice WL,

Müller-Plathe F, Murat M, Paul W, Santos S, Suter UW, Tries V. Bridging

the gap between atomistic and coarse-grained models of polymers: Status and

Perspectives. 2000;152.

3. Lyubartsev AP, Karttunen M, Vattulainen I, Laaksonen A. On coarse-graining

by the inverse Monte Carlo method: dissipative particle dynamics simulations

made to a precise tool in soft matter modeling. Soft Materials. 2003;1:121-137.

4. Tan RKZ, Petrov AS, Harvey SC. YUP: a molecular simulation program for

coarse-grained and multiscaled models. J Chem Theory Comput. 2006;2:529-

540.

5. Clementi C. Coarse-grained models of protein folding: toy models or predictive

tools? Current Opinion in Structural Biology. 2008;18:10-15.

6. Depa PK, Maranas JK. Speed up of dynamic observables in coarse-grained

molecular-dynamics simulations of unentangled polymers. J Phys Chem.

2005;123(9):094901.

7. Coifman R, Maggioni M. Diffusion wavelets. Appl Comput Harm Anal.

2006;21(1):53–94.

8. Bremer J, Coifman R, Maggioni M, Szlam A. Diffusion wavelet packets. Appl

Comput Harm Anal. 2006;21(1):95–112.

9. Abrams C, Delle Site L, Kremer K. Dual-resolution coarse-grained simulation of the bisphenol-A-polycarbonate/nickel interface. Phys Rev E.

2003;67:021807.

10. Abrams C, Kremer K. Combined coarse-grained and atomistic simulation of

liquid bisphenol A-polycarbonate: Liquid packing and intramolecular structure.

Macromolecules. 2003;36:260.

38
Approved for public release; distribution is unlimited.

11. Marrink SJ, de Vries AH, Mark AE. Coarse grained model for semiquantitative lipid simulations. J Phys Chem B. 2004;108(2):750-760.

12. Marrink SJ, Risselada HJ, Yefimov S, Tieleman DP, de Vries AH. The martini force field:âAL’ coarse grained model for biomolecular simulations. J Phys  ̆

Chem B. 2007;111(27):7812-7824.

13. Monticelli L, Kandasamy SK, Periole X, Larson RG, Tieleman DP, Marrink SJ. The MARTINI coarse-grained force field: extension to proteins. J Chem

Theory Comput. 2008;4:819-834.

14. Gautieri A, Russo A, Vesentini S, Redaelli A, Buehler MJ. Coarse-grained

model of collagen molecules using an extended MARTINI force field. J Chem

Theory Comput. 2010;6:1210-1218.

15. Fan B, Maranas J. Coarse-grained simulation of cellulose iÎš with application

to long fibrils. Cellulose. 2015;22(1):31-44.

16. Lu L, Izvekov S, Das A, Andersen HC, Voth GA. Efficient, regularized, and

scalable algorithms for multiscale coarse-graining. J Chem Theory Comput.

2010;6(3):954–965.

17. Harmandaris VA, Reith D, Van der Vegt NFA, Kremer K. Comparison between coarse-graining models for polymer systems: Two mapping schemes for

polystyrene. Macromolecular Chem Phys. 2007;208(2109).

18. Espanol P. Dissipative particle dynamics with energy conservation. Europhys.

Lett. 1997;40(6):631-636.

19. Binder K, editor. Monte Carlo and molecular dynamics simulations in polymer

science. Oxford (UK): Oxford University Press; 1995.

20. Daubechies I. Ten lectures on wavelets. Philadelphia: SIAM; 1992. Proceedings of the CBMS-NSF Regional Conference Series in Applied Mathematics;

vol. 61.

21. Strang G, Nguyen T. Wavelets and filter banks. Cambridge (MA): WellesleyCambridge Press; 1996.

39
Approved for public release; distribution is unlimited.

22. Heyden A, Truhlar DG. Conservative algorithm for an adaptive change of resolution in mixed atomistic/coarse-grained multiscale simulations. J Chem Theory Comput. 2008;4:217-221.

23. Praprotnik M, Delle Site L, Kremer K. Adaptive resolution moleculardynamics simulation: changing the degrees of freedom on the fly. J Chem Phys.

2005;123(224106).

24. Peter C, Delle Site L, Kremer K. Classical simulations from the atomistic

to the mesoscale and back: coarse graining an azobenzene liquid crystal. Soft

Matter. 2008;4:859-869.

25. Ismail AE, Rutledge GC, Stephanopoulos G. Topological coarse graining of

polymer chains using wavelet-accelerated Monte Carlo. I. freely jointed chains.

J Chem Phys. 2005;122:234901.

26. Ismail AE, Stephanopoulos G, Rutledge GC. Topological coarse graining of

polymer chains using wavelet-accelerated Monte Carlo. II. Self-avoiding chains.

J Chem Phys. 2005;122:234902.

27. Chung FRK. Spectral graph theory. American Mathematical Society; 1997.

28. Sun H, Mumby SJ, Maple JR, Hagler AT. An ab initio cff93 all-atom force

field for polycarbonates. J Am Chem Soc. 1994;116(7):2978–2987.

29. Plimpton S. Fast parallel algorithms for short-range molecular dynamics. J

Am Chem Soc. 1995;117(1):1–19; .

30. Vazquez E, Bect J. Convergence properties of the expected improvement algorithm with fixed mean and covariance functions. Journal of Statistical Planning

and Inference. 2010;140(11):3088 - 3095; .

31. Katsoulakis MA, Plechác P, Rey-Bellet L, Tsagkarogiannis DK. Mathematical ˇ

strategies in the coarse-graining of extensive systems: error quantification and

adaptivity. J Non-Newtonian Fluid Mech. 2008;152:101-112.

32. Shell MS. The relative entropy is fundamental to multiscale and inverse thermodynamic problems. J Chem Phys. 2008;129(14).

40
Approved for public release; distribution is unlimited.

33. Katsoulakis MA, Plechác P. Information-theoretic tools for parametrized ˇ

coarse-graining of non-equilibrium extended systems. J Chem Phys.

2013;139(7).

34. Mansour AA, Ortoleva PJ. Implicit time integration for multiscale molecular dynamics using transcendental Padé approximants. J Chem Theory Comput.

2016;(in press).

41
Approved for public release; distribution is unlimited.

INTENTIONALLY LEFT BLANK.

42
Approved for public release; distribution is unlimited.

List of Symbols, Abbreviations, and Acronyms

CAM collective action mode

CG coarse-grained

DoF degree of freedom

GAFF generalized amber force field

LAMMPS Large-scale Atomistic/Molecular Massively Parallel Simulator

MD molecular dynamics

PCFF polymer-consistent force field

PE polyethylene

43
Approved for public release; distribution is unlimited.

INTENTIONALLY LEFT BLANK.

44
Approved for public release; distribution is unlimited.

1

(PDF)

DEFENSE TECHNICAL

INFORMATION CTR

DTIC OCA

2

(PDF)

DIRECTOR

US ARMY RESEARCH LAB

RDRL CIO L

IMAL HRA MAIL & RECORDS MGMT

1

(PDF)

GOVT PRINTG OFC

A MALHOTRA

4

(PDF)

DIR USARL

RDRL WML B

B RICE

RDRL WMM G

J ANDZELM

T SIRK

J LENHART

45
Approved for public release; distribution is unlimited.

INTENTIONALLY LEFT BLANK.

46