Acceleration of Mean Square Distance Calculations
with Floating Close Structure in Metadynamics
Simulations

arXiv:1801.02362v1 [cs.CE] 8 Jan 2018

Jana Pazúrikováa,∗, Jaroslav Oľhaa , Aleš Křeneka , Vojtěch Spiwokb
a

Institute of Computer Science, Masaryk University, Botanická 554/68a, 602 00 Brno,
Czech Republic
b
Department of Biochemistry and Microbiology, University of Chemistry and Technology,
Technická 5, 166 28 Prague 6, Czech Republic

Abstract
Molecular dynamics simulates the movements of atoms. Due to its high cost,
many methods have been developed to “push the simulation forward”. One
of them, metadynamics, can hasten the molecular dynamics with the help
of variables describing the simulated process. However, the evaluation of
these variables can include numerous mean square distance calculations that
introduce substantial computational demands, thus jeopardize the benefit of
the approach. Recently, we proposed an approximative method that significantly reduces the number of these distance calculations. Here we evaluate
the performance and the scalability on two molecular systems. We assess
the maximal theoretical speed-up based on the reduction of distance computations and Ahmdal’s law and compare it to the practical speed-up achieved
with our implementation.
Keywords: approximation, acceleration, molecular dynamics simulation,
metadynamics

1. Introduction
Molecular dynamics [1] is a well-known computational method which simulates the behavior of a molecule at the atomic scale by integrating Newton’s
∗

Corresponding author
Email address: pazurikova@ics.muni.cz (Jana Pazúriková)

1

equations of motion for all atoms over time. Due to a high frequency of movements at this scale, the integration step must not exceed femtoseconds, which
yields realistic computational speed to be at most nanoseconds of simulated
time per day of computation time. On the other hand, biologically relevant
phenomena take milliseconds and more at the macromolecular level, hence
their complete simulations are still not feasible nowadays. The necessity of
long simulation timescales can be traced down to high energy barriers between the distinct states of the macromolecule, which require a considerable
momentum to be gained at the nanoscale before they are overcome.
However, various classes of problems do not require a realistic trace of
the simulation, it is sufficient to explore the state space in which the simulation would travel eventually [2]. Therefore, the simulation run can be
“pushed forward” by additional steering, e.g. in the method called metadynamics by adding bias potential in already visited areas, or by “filling energy
minima with computational sand”, as illustrated in Figure 1.

Figure 1: Filling potential energy minima with a bias potential. The simulation starts with
x = 0, and it accumulates the bias potential in this area first. After 20 steps the energy
barrier at x = −1.75 is overcome and the bias potential starts filling the minimum at
x = −4. After approx. 300 steps the state space is explored. Taken from [3], image
copyright (2002) National Academy of Sciences.

2

Many methods have been developed to accelerate molecular dynamics
simulations by artificial intervensions [4]. Unfortunately, the most promising
ones often introduce a non-negligible computational overhead to calculate
the bias potential, hence jeopardize the benefit of the whole approach. Recently, we proposed a method [5] which considerably reduces the number
of required calls to the most expensive routine of such calculation—mean
square distance (MSD) between two molecular shapes, while not affecting
the overall accuracy. Here we introduce the ready-to-use implementation of
the method (replacing a rather inefficient prototype from [5]) and evaluate
the practical benefits—speed-up of the real computation, including the use
of multiple CPU cores.
2. Background
Simulations of molecular dynamics (MD) calculate the movements of
atoms caused by their interactions [1]. First, the potential energy of the system has to be computed. In this work, we modeled the interatomic potential
using molecular mechanics model, considering atoms as mass charged points
and approximating both bonded and nonbonded interactions with a series
of empirical functions and parameters. Second, the forces (determined by
the potential) are integrated with Newton’s law of motion to obtain the positions of atoms in the next step.
Computational demands of molecular dynamics grow with the number
of atoms and simulation steps. The calculation of the interatomic potential
scales as O(Natoms log(Natoms )) with common methods. The parallelization
through the spatial domain makes it possible to simulate even millions of
atoms. On the other hand, the timescales of simulations are more difficult
to handle. A tiny timestep of the integration scheme (in the order of 10−15 s)
makes it hardly feasible to simulate more than a few microseconds [6], a
millisecond at most [7], although many biologically or chemically relevant
processes take longer than that. We need longer simulation times because
we need to explore the whole energy landscape, ideally all possible states of
the molecular system. However, molecules tend to stay in the local energy
minima and they cross the high-energy barrier between them only once in
a while.
Many methods have been developed to accelerate the rare events of those
crossings, for an overview see [8, 4, 9]. They enhance the sampling of the energy landscape and explore it much faster using various artificial interven3

tions to the simulation. One of the most common, metadynamics [3] (MTD),
fills the minima in the energy landscape with an artifical bias potential. It
“spreads the sand into the valleys”, and thus makes it easier to cross the peaks.
To properly function, metadynamics needs a few collective variables (CV)—
preselected degrees of freedom that differentiate between the states of the simulated process and encompass all the relevant motions. Moreover, they have
to be limited in number, otherwise, both the performance and the accuracy
of metadynamics will decline [10, 11]. Collective variables are usually computed from the coordinates of the molecular system at the given time, further
referred to as the current structure or x.
Metadynamics accumulates a bias potential in the form of Gaussian hills
centered in the values of collective variables, see Equation (1). In Figure 1, the sum of Gaussian hills forms curvy lines depicting the borders
of the summed bias potential in various time points. Well-chosen CVs direct
to the right spot for “the sand to be spread on”. The bias potential expresses
as the bias force, so collective variables need to be differentiable with respect
to coordinates.

Vbias (S(x), t) =

X

Y

t0 =τG ,2τG ,...
t<t0

i



0
(Si (x(t)) − sti )2
w exp −
2σ 2
(1)
Fbias (t) =

∂Vbias (S(x), t)
∂x

where x are coordinates of the current structure, S(x) is the vector of functions computing collective variables; τG is the frequency of Gaussian hills
0
addition; w and σ is the height and the width of the Gaussian hill; sti is
the value of i-th CV at time t0 , i.e. Si (x(t0 )).
The choice or design of CVs requires physico-chemical expertise and naturally gets more complicated with complex processes. A distance between
specific atoms, angles or their combinations work well for simple processes.
More complex ones, like protein folding, require a more sophisticated approach such as path CVs [12] or their generalization, property map CVs [13].
They are based on comparison to reference landmark structures, a series of
snapshots capturing various states of the process. Reference structures can
be generated with any non-continuous simulations, non-physical actions, high
temperatures, annealing techniques or even manual construction in a visualization software.
4

The value of the property map collective variable corresponds to a weighted
average of properties of the closest reference structures, where the weight is
determined by the distance, see Equation (2). For example, if we take an index of the reference structure as a property, the value of CV would be close
to the index of the closest reference structure.
N
P

S(x) =

qi exp (−λD(x, ai ))

i=1
N
P

(2)
exp (−λD(x, ai ))

i=1

where N is the number of reference structures, ai are coordinates of i-th reference structure, qi is an arbitrary property of the given reference structure,
D() is the distance function, λ is a tuning parameter.
Clearly, the evaluation of property map CVs requires as many distance
computations in every step of the simulation as there are reference structures.
For smooth processes, we can reduce the number of distance computations
by employing a neighbourlist (NL). Since the contribution to a CV’s value decreases exponentially with the distance, without sudden changes in system’s
structure it suffices to compute distance only to the few closest reference
structures in the NL and update the list regularly. However, in processes with
abrupt changes, we need to compute the distances of the current structure to
all the reference structures in each step (referred to as no neighbourlist). In
both cases, numerous distance computations introduce high computational
cost.
3. Floating Close Structure
The standard distance measure between molecular systems in computational sciences is the (root) mean square distance between superimposed,
fitted structures [14, ch.16]. The fitting between structures minimizes the distance and consists of two tasks: translation and rotation. Shifting both
structures to the same location in the coordinate system means subtracting
coordinates of the center of mass from coordinates of all atoms. As this is
both simple and fast, we omit it from further equations for simplicity and
assume that it was done in advance for all the structures. On the contrary,
obtaining the rotation matrix presents the most expensive task in the distance computation. It involves demanding matrix operations, such as the
5

matrix diagonalization from quaternions as in [15, 16] or the singular value
decomposition as in [17].
To get a closer look on the bottlenecks, we analyzed the performance of
simulations with a profiling software. We ran simulations using Gromacs [18]
as the molecular dynamics software and Plumed [19] as the metadynamics
tool, details in section 5. Even with the neighbourlist employed, distance
computation takes up to 43% of miniprotein Trp-cage simulation walltime.
For small 24-atom cyclooctane simulation without the neighbourlist, it takes
93% of walltime, see Table 1. The majority of that time is spent in BLAS
function DSYEVR, diagonalizing matrix in Kearsley’s method [15] implemented in Plumed. As BLAS is already a finely-tuned routine, we assume no
further code optimizations would bring a substantial improvement. Offloading to GPU brings its own issues, as we outlined in [20]. Therefore, we aim
to modify the entire method in a way that reduces the number of demanding
distance computations and thus decreases the computational cost.
3.1. Exact MSD Calculation
The original method for exact MSD computation calculates accurately
but expensively the rotation matrix for all the reference structures in the neighbourlist in every step, see Algorithm 1 and Figure 2. Then it uses rotation
matrices to compute the distance, see Equation (3), and derivatives with
respect to coordinates (to evaluate the bias force), see Equation (4).

x

ai

Figure 2: Original method expensively computes distances from the current structure
(x) to reference structures in neighbourlist (ai ) in every step. Arrows denote passing
simulation time.

6

Algorithm 1 Computation of the collective variables in the original method
1: function compute_cvs(x, ref_structures)
2:
loop through ai in neighlist
3:
Rxai ← compute_R(x, ai )
4:
distance ← Equation (3)
5:
derivatives ← Equation (4)

D(x, ai ) =

NX
atoms

wj kdj k2

(3)

j

where dj = xj − Rai x aij
∂D(x, ai )
= 2wk dk + wk0
∂xk

!
X

−2wj dj
!

X

−2wj dj × aij

j

+

j

⊗

∂Rai x
∂xk

(4)

where x are coordinates of the current structure, ai are coordinates of i-th
reference structure, Rai x is the rotation matrix fitting the reference structure
ai onto the current structure x, w and w0 are weights, ⊗ is element-wise
multiplication.
Two sets of weights influence the computation: displacement w that describes the contribution of an atom’s displacement to the total distance, and
alignment w0 that describes the contribution to the center of mass which
influences the translation.
During a usual step, we need to compute the expensive rotation matrix
M times, where M is the neighbourlist size. During an update of NL, every L steps, we calculate the distance (thus also the rotation matrix) to
all the reference structures. More details about the computational cost in
section 6.2.
3.2. Approximative MSD Calculation
Our method is based on the assumption that the molecular system moves
reasonably smoothly in dynamics simulations due to a tiny timestep. As its
7

coordinates alter only slowly, the values in the rotation matrices should also
change only slightly every step. Therefore, we reuse them in a few consecutive
steps, until the current structure moves too far.
At the beginning of the simulation, we start as usual, see Figure 3 and
Algorithm 2. We assign current coordinates to a close structure y and save
computed rotation matrices Rai y . In most steps, we reuse these saved matrices to approximate Rai x ≈ Rxy Rai y while the current structure x and
the close structure y stay close to each other. If the distance D(x, y) exceeds
a given threshold ε, we reassign the close structure, y ← x and recompute
the matrices Rai y .
Algorithm 2 Approximative computation of the collective variables with
the close structure
1: function compute_cvs(x, ref_structures)
2:
Rxy ← compute_R(x,y)
3:
Dxy ← Equation (3)
4:
∂Rxy /∂xk ← Equation(4)
5:
if Dxy > ε then
6:
loop through all ai
7:
Rai x ← compute_R(x, ai )
8:
save Rai x
9:
distance ← Equation (3)
10:
derivatives ← Equation (4)
11:
else
12:
loop through ai in neighlist
13:
Rai x ← Rxy Rai y
14:
distance ← Equation (5)
15:
derivatives ← Equation (6)

We use exact MSD calculation, Equations (3) and (4), in the reassignment
steps and for the computation of D(x, y). In the majority of steps, we use
Equation (5) for distance and Equation (6) for derivatives. Rxy Rai y approximates Rai x , omitting the demanding accurate computation of the rotation
matrix.

8

D(x,y) >
y⬅x

y⬅x

x

ai

Figure 3: Close structure method expensively computes distances only between the current
(x) and the close structure (y) and then during reassignments when the distances to all
reference structures (ai ) are accurately calculated. Arrows denote passing simulation time.

D̃(x, ai ) =

NX
atoms

wj kd̃j k2

(5)

j

d̃j = xj − Rxy Rai y aij
∂ D̃(x, ai )
= 2wk d˜k + wk0
∂xk

!
X

−2wj d̃j

j

!
X

−2wj d̃j × Rai y aij

j

+
(6)

∂Rxy
⊗
∂xk

where y is the close structure, Rxy is the rotation matrix that fits the close
structure to the current structure, Rai y is the saved rotation matrix that fits
the close structure to the reference structure ai .
The neighbourlist technique, not detailed in Algorithm 2 for simplicity,
can be applied as usual in conjunction with the close structure method. However, we do not recompute the distance to all reference structures in the NL
update step, as it is in the original method. Instead, we recompute during
the reassignment of the close structure. Therefore, in the close structure
method, we select the structures during the NL update according to the approximated distance.
9

Close structure approximation reduces the number of distance computations in a usual step to one—evaluating D(x, y). When reassigning the close
structure, we expensively calculate the distance to all N reference structures.
More details about the computational cost in section 6.2.
4. Implementation
We implemented the close structure method into Plumed [19], a standard
tool for metadynamics. Plumed works as a plugin into many molecular dynamics softwares, we used it with Gromacs [18]. We changed two classes, one
that represents property map and path CVs, and another that encompasses
RMSD computation. Since the prototype implementation of our method, we
have dealt with issues regarding code optimization, data structures, memory access patterns, communication through MPI and scalability issues. The
code ready for production use is available in the version v2.4 of Plumed, see
https://github.com/plumed/plumed2.
For performance evaluation, we combined Gromacs v5.1.4 with Plumed
v2.3. We applied the same level of minor code optimizations, enforcing
the vectorization of loops and inlining of template functions, both on our
modified and the original code in parts regarding the evaluation of property
map collective variables. We did not include any aggressive optimizations
that would transform data structures since such heavy changes to software
design would make it difficult to incorporate the code into the official software [21]. We compiled the code with Intel compiler version 17 and SIMD
instruction set AVX_256.
5. Computational Details
In this paper, we focus on the evaluation of the acceleration and the
scalability. We reuse two molecular systems and their simulation setup from
[5] where we have thoroughly evaluated the accuracy. Thus, we describe
datasets only shortly, for further details see [5].
The first molecular system, a non-symmetric trans,trans-1,2,4-triflourocyclooctane (referred to as cyclooctane throughout the paper) can form different
conformations with rapid transitions. As it contains only 24 atoms and is simulated in the vacuum (ergo no further molecules of the solvent), the usually
high computational cost of molecular dynamics stays low for this molecule.
However, abrupt and rare changes in structure require the computation of
10

distances to all 521 reference structures in every step, because the neighbourlist would not be able to keep up. As distance computations take a vast
majority of the computation time, it serves as an example of a molecule that
should greatly benefit from the acceleration of CV computation.
Second, 304-atom miniprotein Trp-cage, PDB ID 1L2Y, simulated in
an implicit solvent smoothly diffuses between several conformations in its
slow folding. The molecular dynamics, especially long-range electrostatic interactions, takes the majority of the computation time. Even with the neighbourlist, the closest 50 structures chosen from all 2120 reference structures
every 50 steps, the contribution of metadynamics is not negligible, it takes
almost half of the computation time. The trp-cage simulation shows a more
common example where the acceleration of metadynamics speeds up the whole
simulation only modestly. However, reducing computational overhead brought
by metadynamics enables us to simulate at the speed of classical molecular
dynamics with the advantage of rare events’ acceleration with metadynamics.
We have set up the simulation as described in [5]. All simulations ran
for 100,000 steps. Metadynamics simulation has been done without a neighbourlist for cyclooctane, i.e.distances to all 521 reference structures calculated every step, with 50-item NL updated every 50 steps for Trp-cage. In
all simulations, the bias potential has been computed with the grid [22],
diminishing the cost of a Gaussian hill summation growing with longer simulations. In close structure simulations, we reassign the close structure if
D(x, y) > ε, where ε = 0.01 nm2 .
Experiments for this paper differ from those in [5] in versions of software
and our method’s implementation. Accuracy evaluation there has been done
on Gromacs 4.5.7 and Plumed 2.1 with the prototype implementation of our
method. Performance evaluation here has been done on Gromacs 5.1.4 and
Plumed 2.3 with ready-to-use, optimized implementation of our method. All
performance evaluation experiments were executed on one machine with 2x
8-core Intel Xeon E5-2650 v2 2.6 GHz, even multiple MPI processes were
assigned to the same machine.
6. Results and Discussion
The main aim of the close structure method is to accelerate metadynamics
simulations by reducing the number of distance computations. Therefore, in
its evaluation, we focused both on the theoretical and the practical speed-up.

11

We have evaluated the accuracy of the close structure method meticulously
in [5], here we only shortly outline the main findings.
6.1. Accuracy
Because of the chaotic character of molecular dynamics, even the smallest changes caused by the close structure approximation cause subsequent
significant changes in trajectory. Since the step-by-step comparison of trajectories would be meaningless, we inspected the similarity of explored energy
landscapes. First, we evaluated both Equation (3) and (5) on a long series
of x and corresponding Rxy , finding almost perfect correlation. Second, we
visualized property map collective variables computed during the simulation
with the original method and with the close structure method. Both showed
a great resemblance. Finally, we analyzed essential coordinates [23] of both
trajectories. We found a match not only between the essential movements
of the system but also in minor harmonic oscillations. We concluded that
the trajectory of the close structure simulation explored the same energy
landscape as the simulation with the original method.
6.2. Theoretical Speedup
For theoretical evaluation of speed-up, we inspected the number of expensive distance computations and conditions for speed-up.In the next section,
we calculate the theoretical speed-up according to Ahmdal’s law for both
evaluated datasets.
By an expensive distance computation we mean finding the rotation matrix with Kearsley’s method [15] and calculating the distance and derivatives
with Equations (3) and (4). In the original method, all distances are calculated expensively. In the close structure method, the majority of distances
are just cheaply approximated with saved rotation matrices using Equations
(5) and (6).
In the original method, in a usual step, we need M (size of the neighbourlist, or M = 0 in case of no neighborlist) distance computations. In
a neighbourlist update step, every L steps (L = 1 in case of no neighborlist),
we compute N (number of reference structures) distances. On average, we
expensive distance computations per step.
need M + N −M
L
In the close structure method, in a usual step, we need only one MSD
computation, the rest is approximated. Furthermore, in a close structure
reassignment step, we expensively recompute N rotation matrices, so that
we keep the accuracy of distance computations for neighborlist updates. On
12

N
expensive distance computations per step, where
average, we need 1 + K
the close structure is reassigned (on average) every K steps.
The close structure method will accelerate metadynamics if M + N −M
>
L
N
1 + K . From our experiments, if ε is not too small, then K  L, as
the close structure was reassigned on average every few thousand steps for
ε = 0.01 nm2 (10000 for cyclooctane, 3000 for Trp-cage) in the 100,000step simulations we did. The average number of reassignments depends on
how often and how fast the structure of molecule changes, so it varies with
the threshold ε, but also with the length of the simulation, the starting structure, temperature and other simulation parameters.
For the simulations without the neighbourlist, the comparison is reduced
N
which clearly applies as K  1. Therefore, the number of
to N > 1 + K
distance computations is reduced approximately by an order of N .
For the simulations with the neighbourlist, the specific values of M , L, N
and K determine the acceleration. It is safe to assume M  N and L  N .
If K > L, which is quite easily achievable for reasonable values of ε, then
N −M
N
> K
. Thus, the number of distance computations in simulations with
L
the NL can be reduced up to an order of M .
Naturally, the speed-up of distance computation accelerates the whole
simulation only to some extent. This can be assessed using Ahmdal’s law,
see Equation (7). The portion p represents the percentage of time spent in
functions we have accelerated, i.e. the distance computations, in the original
implementation. The speed-up s denotes the acceleration of this portion in
our modified implementation. The whole computation can achieve only the
speed-up S. For example, even an infinite speed-up of distance computation would accelerate the whole simulation only ten times, if the distance
computations originally took up 90% of computation time.

S=

1
(1 − p) +

p
s

(7)

where p is the original portion of time spent in the given function, s is
the speed-up of the accelerated function.
More general performance behavior can be deduced from the analysis
above. First, the acceleration brought by the close structure method decreases with the increased cost of molecular dynamics due to Ahmdal’s law.
This applies especially in the case of larger molecules. The cost of computing MSD with the original method also grows with the number of atoms,
but to a much lesser extent than the demands of MD. Therefore, a higher
13

acceleration is expected with smaller molecules. With larger molecules, the
best we can achieve is to get closer to the performance of molecular dynamics
without metadynamics. Second, the acceleration brought by the close structure method increases with the decreasing stride of neighborlist’s update (L),
but does not depend much on its size. Therefore, a higher acceleration is expected in the case of no neighborlist. And finally, the acceleration brought by
the close structure method decreases with the increased number of reference
structures, especially with the neighborlist updated often.
6.3. Theoretical Speedup for Evaluated Datasets
We evaluated the theoretical speed-up expected in simulations with our
datasets. Moreover, we assessed the speed-up of the whole simulation using
Ahmdal’s law and information from profiling.
For cyclooctane, we need 512 MSD computations each step with the original method and an average of 1 + 512/10000 = 1.05 with the close structure
method. The theoretical speed-up for MSD computation reaches 512/1.05 =
488.
For Trp-cage, we need on average 50 + (2120 − 50)/50 = 91.4 MSD computations with the original method and 1 + 2120/3000 = 1.7 with the close
structure method. The theoretical speed-up for MSD computation reaches
91.4/1.7 = 54.
The theoretical speed-up for the whole simulation can be assessed using
Ahmdal’s law and information from profiling. We combined the information from log files (time spent with whole metadynamics and evaluation of
collective variables) and from Intel VTune Profiler (time spent computing
the distance) to assess the percentage of computation time spent in various
functions in cyclooctane and Trp-cage simulations with one MPI process and
one OpenMP thread, see Table 1. Metadynamics (MTD) includes the evaluation of collective variables (CV), initialization, logging and addition of
Gaussian hills. CV evaluation requires distance computations (MSD) before
the calculation of its own value. MSD computation consists of matrix diagonalization (DSYEVR), the calculation of the distance and the derivatives.
The speed-up SM SD represents the maximal theoretical speed-up of the whole
simulation according to Ahmdal’s law (see Equation (7)) taking into account
the speed-up of MSD computation as assessed above.
As expected, cyclooctane simulation would greatly benefit from the acceleration due to a smaller size and no neighborlist. Trp-cage simulation can

14

Table 1: Time portions spent in various metadynamics functions and the maximal speedup of the whole simulation

MTD
CV MSD DSYEVR SM SD
cyclooctane 99.7% 98.5% 93%
70%
14
Trp-cage
53%
50% 43%
23%
1.7
achieve only modest improvement as the most of the computational cost lies
in MD code.
6.4. Practical Speedup and Scalability
For the practical evaluation of speed-up and scalability, we ran simulations without metadynamics (Gromacs), with metadynamics and the original
method for MSD computation (Gromacs+Plumed) and with metadynamics modified with the close structure method (Gromacs+modified Plumed).
Furthermore, we executed them for various numbers of MPI processes and
OpenMP threads. As a speed measure, we took ns/day stated by Gromacs
in its logs. As a reference point for the speed-up, we considered Gromacs
with original Plumed ran with one MPI process and one OpenMP thread.
Out of four runs for each variant and MPI/OpenMP configuration, we
considered the minimal running time, i.e. the highest speed, to eliminate
the random interference of the operating system’s background activity.
The upper parts of Figures 4 and 5 show the speed-ups for the basic case
of one MPI process and one OpenMP thread. For both molecular systems,
the close structure method almost reaches the maximal theoretical speed-up
as calculated in the previous subsection.
Cyclooctane simulations with the close structure method achieved a speedup of an order of magnitude, notice the logarithmic scale on the x-axis of
Figure 4. The close structure method almost reached to the maximal theoretical performance, albeit it still lags behind Gromacs without metadynamics.
That is caused by Plumed overhead like initialization, logging, Gaussian hills
addition and CV evaluation.
Trp-cage simulations with the close structure method achieved a more
modest speed-up due to the major cost of molecular dynamics. Furthermore, the practical and theoretical speed-up differ more than in cyclooctane
simulation. We suggest it is caused by a smaller portion of time spent in matrix diagonalization function, as shown in Table 1. The distance computation
takes 43% of computation time. Out of that, only half is spent in diagonal15

Figure 4: Speedup and strong scalability of cyclooctane simulations. Notice the logarithmic scale of x axis, the reference point is the blue bar with the original method. The grey
bar with molecular dynamics without metadynamics puts into perspective how metadynamics increases the cost.

ization, as opposed to three-quarters in the case of cyclooctane. The rest is
spent in the computation of distance and especially derivatives, which scales
with the size of the molecule. We reduced the number of calls to DSYEVR

16

Figure 5: Speedup and strong scalability of the Trp-cage simulations. The reference point
is the blue bar with the original method. The grey bar with molecular dynamics without
metadynamics puts into perspective how metadynamics increases the cost.

with the close structure method, but the other half of computation remained
unchanged: both Gromacs+Plumed and Gromacs+modified Plumed contain
the same code optimizations. Nevertheless, we got much closer to the performance of Gromacs alone, thus gaining the advantage of enhanced sampling
17

for only a small cost.
The lower parts of Figures 4 and 5 show the speed-ups of strong scaling
for various combinations of MPI processes and OpenMP threads.
First, compare the performances with the same total number of threads,
i.e. 1/4 with 2/2, 1/8 with 2/4 and 1/16 with 2/8. In all cases for the original method, two MPI processes perform better than one with twice as many
threads. The reason for that stems from the implementation of CV calculation in Plumed. The loop over reference structures in the neighbourlist
(line 2 in Algorithm 1) is parallelized only with MPI. As the close structure
method accelerates the computation within that loop, the second MPI process does not bring any advantage. In the case of Trp-cage, it even slows
down the computation due to the communication and synchronization cost.
Decreasing speed with more resources (8 total threads for cyclooctane, 16
for Trp-cage) suggests that the overhead caused by parallelization exceeds
its benefits.
Second, compare performances for one MPI process, i.e. 1/4, 1/8 and
1/16. Any speed-up here can come only from faster molecular dynamics.
For cyclooctane simulations, a small amount of computation saturates even
a single core, thus performance stays about the same for all three variants.
For Trp-cage simulation, some improvement appears as molecular dynamics
requires more computation due to the protein’s size. The simulations with
the close structure method also exhibit this trend.
Overall, the close structure method reduces the gain of another MPI
process but does not interfere with the general scalability behavior.
7. Conclusion
In this work, we address high computational demands of distance calculations in path/property map collective variables in metadynamics. The original method expensively evaluates the rotation matrices needed to superimpose the current structure to the reference structures in each step for many
reference structures. By introducing a close structure, we reuse those matrices in next steps and cheaply approximate the vast majority of distance
computations.
We thoroughly evaluated accuracy elsewhere [5] and concluded that simulations with the close structure method explore the same energy landscape
as simulations with the original method, despite different trajectories.

18

Here, we presented the tuned implementation of our method ready for
production use and evaluated its performance. We calculated the number
of distance computations for both the original and close structure method
and assessed the theoretical speed-up of the whole simulation. The number
of expensive MSD computations decreased approximately by the order of
the neighbourlist size or the number of all the reference structures in case of
no neighborlist. Clearly, the accelerated distance computation is expressed
in the speed-up of the whole simulation only to the extent proportional to
the time spent on distance computation in the original method. That depends on the neighborlist size and update stride, but mostly on the cost of
the molecular dynamics part, determined by the molecule’s size.
The experiments showed that the actual speed-up closely approaches
the theoretical one for both evaluated molecular systems. Apart from the close
structure approximation, minor code optimizations in distance computation
contributed to this, including enforcing of loop vectorization and inlining of
frequently called functions. These optimizations did not help the original
method, as it spent the majority of time in the finely-tuned BLAS routine. However, unvectorized loops became prominent after the employment
of the close structure, so even minor code optimizations pushed the practical
speed-up further.
The close structure method affects the scalability of Plumed code, diminishing the gain of an additional MPI process. This was expected as we
reduced the amount of computation so that even one core suffices.
Overall, the positive effects of the close structure approximation are
prominent especially in smaller molecules with frequently updated neighborlist or no neighborlist at all. With larger molecules, the method diminishes the overhead of metadynamics and approaches the performance of sole
molecular dynamics simulation with the advantage of enhanced sampling.
The error brought by the approximation of the close structure method does
not negatively influence the simulation, the explored landscapes remain the
same.
Here we presented a typical application of the close structure method,
molecular dynamics with metadynamics using property map collective variables to calculate the bias potential. The close structure method would
benefit any CV-based method of enhanced sampling, provided it uses mean
square distance as a measure. Our method significantly reduces the overhead of metadynamics in simulations, getting closer to the ideal situation:
molecular dynamics with the advantage of enhanced sampling—but without
19

significant additional cost.
Acknowledgments
This work was supported by Czech Science Foundation (15-17269S) and
LM2015047 Czech National Infrastructure for Biological Data. Computational resources were provided by the CESNET LM2015042 and the CERIT
Scientific Cloud LM2015085, provided under the programme “Projects of
Large Research, Development, and Innovations Infrastructures”.
References
[1] F. Jensen, Introduction to computational chemistry, 2nd Edition, John
Wiley & Sons Ltd, Great Britain, 2007.
[2] R. C. Bernardi, M. C. Melo, K. Schulten, Enhanced sampling techniques
in molecular dynamics simulations of biological systems, Biochimica et
Biophysica Acta (BBA) - General Subjects 1850 (5) (2015) 872–877.
doi:10.1016/j.bbagen.2014.10.019.
[3] A. Laio, M. Parrinello, Escaping free-energy minima., Proceedings of the
National Academy of Sciences of the United States of America 99 (20)
(2002) 12562–12566. arXiv:0208352, doi:10.1073/pnas.202427399.
URL http://www.pnas.org/content/99/20/12562
[4] P. Tiwary, A. van de Walle, A review of enhanced sampling approaches
for accelerated molecular dynamics, in: Springer Series in Materials
Science, Vol. 245, Springer International Publishing, 2016, pp. 195–221.
doi:10.1007/978-3-319-33480-6_6.
[5] J. Pazúriková, A. Křenek, V. Spiwok, M. Šimková, Reducing the number
of mean-square deviation calculations with floating close structure in
metadynamics, Journal of Chemical Physics 146 (11) (2017) 115101.
doi:10.1063/1.4978296.
URL http://aip.scitation.org/doi/10.1063/1.4978296
[6] J. D. Durrant, R. M. Bush, R. E. Amaro, Microsecond Molecular Dynamics Simulations of Influenza Neuraminidase Suggest a Mechanism
for the Increased Virulence of Stalk-Deletion Mutants, Journal of Physical Chemistry B 120 (33) (2016) 8590–8599. doi:10.1021/acs.jpcb.
6b02655.
20

[7] K. Lindorff-Larsen, P. Maragakis, S. Piana, D. E. Shaw, Picosecond
to Millisecond Structural Dynamics in Human Ubiquitin, Journal of
Physical Chemistry B 120 (33) (2016) 8313–8320. doi:10.1021/acs.
jpcb.6b02024.
[8] V. Spiwok, Z. Sucur, P. Hosek, Enhanced sampling techniques in
biomolecular simulations (nov 2015). doi:10.1016/j.biotechadv.
2014.11.011.
[9] A. Laio, F. Gervasio, Metadynamics: a method to simulate rare events
and reconstruct the free energy in biophysics, chemistry and material science, Reports on Progress in Physics 71 (12) (2008) 126601:1–126601:22.
doi:10.1088/0034-4885/71/12/126601.
[10] A. Laio, A. Rodriguez-Fortea, F. L. Gervasio, M. Ceccarelli, M. Parrinello, Assessing the accuracy of metadynamics, Journal of Physical
Chemistry B 109 (14) (2005) 6714–6721. doi:10.1021/jp045424k.
[11] A. Barducci, M. Bonomi, M. Parrinello, Metadynamics, Wiley Interdisciplinary Reviews: Computational Molecular Science 1 (5) (2011)
826–843. doi:10.1002/wcms.31.
[12] D. Branduardi, F. L. Gervasio, M. Parrinello, From A to B in free energy
space, Journal of Chemical Physics 126 (5) (2007) 054103. doi:10.
1063/1.2432340.
[13] V. Spiwok, B. Králová, Metadynamics in the conformational space nonlinearly dimensionally reduced by Isomap, Journal of Chemical Physics
135 (22) (2011) 1–7. doi:10.1063/1.3660208.
[14] J. Gu, P. E. Bourne, Structural bioinformatics, Wiley-Blackwell, 2009.
[15] S. Kearsley, On the orthogonal transformation used for structural comparisons, Acta Crystallographica Section A Foundations of Crystallography 45 (2) (1989) 208–210. doi:10.1107/S0108767388010128.
[16] E. A. Coutsias, C. Seok, K. A. Dill, Using quaternions to calculate
RMSD, Journal of Computational Chemistry 25 (15) (2004) 1849–1857.
arXiv:arXiv:1011.1669v3, doi:10.1002/jcc.20110.

21

[17] W. Kabsch, A solution for the best rotation to relate two sets of vectors, Acta Crystallographica Section A 32 (5) (1976) 922–923. arXiv:
05677394, doi:10.1107/S0567739476001873.
[18] S. Páll, M. J. Abraham, C. Kutzner, B. Hess, E. Lindahl, Tackling exascale software challenges in molecular dynamics simulations with GROMACS, Lecture Notes in Computer Science (including subseries Lecture
Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 8759
(2015) 3–27. arXiv:1506.00716, doi:10.1007/978-3-319-15976-8_
1.
[19] G. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, G. Bussi,
PLUMED 2: New feathers for an old bird, Computer Physics Communications 185 (2) (2014) 604–613. arXiv:1310.0980, doi:10.1016/
j.cpc.2013.09.018.
[20] J. Filipovič, J. Pazúriková, A. Křenek, V. Spiwok, Accelerated RMSD
Calculation for Molecular Metadynamics, in: J. Évora-Gómez, J. J.
Hernandéz-Cabrera (Eds.), 30th European Simulation and Modelling
Conference, Reproduct NV, Ghent, Belgium, 2016, pp. 278–280.
URL https://www.eurosis.org/cms/index.php?q=taxonomy/term/
58
[21] Plumed2.
URL https://github.com/plumed/plumed2
[22] User documention on feature METAD.
URL https://plumed.github.io/doc-v2.2/user-doc/html/_m_e_
t_a_d.html
[23] A. Amadei, A. B. Linssen, H. J. Berendsen, Essential dynamics of proteins., Proteins 17 (4) (1993) 412–25. doi:10.1002/prot.340170408.

22

