arXiv:1502.03234v1 [cs.DC] 11 Feb 2015

Master Thesis - Master of Computer Science

A Scalable, Linear-Time Dynamic Cutoff
Algorithm for Molecular Simulations of
Interfacial Systems

Paul Springer
RWTH Aachen University
Paul.Springer@RWTH-Aachen.de

Supervised by
Prof. Paolo Bientinesi1 , Ph.D. & Prof. Ahmed E. Ismail1 , Ph.D.
1

Institute for Advanced Study in Computational Engineering Science

Declaration of Independence
I hereby declare that I am the sole author of this thesis and that I only used the resources listed in the
reference section.
Aachen, Germany
Mai 31, 2014

Paul Springer

2

Abstract
This master thesis introduces the idea of dynamic cutoffs in molecular
dynamics simulations, based on the distance between particles and the
interface, and presents a solution for detecting interfaces in real-time.
Our dynamic cutoff method (DCM) exhibits a linear-time complexity
as well as nearly ideal weak and strong scaling. The DCM is tailored
for massively parallel architectures and for large interfacial systems
with millions of particles. We implemented the DCM as part of the
LAMMPS open-source molecular dynamics package and demonstrate
the nearly ideal weak- and strong-scaling behavior of this method on
an IBM BlueGene/Q supercomputer. Our results for a liquid/vapor
system consisting of Lennard-Jones particles show that the accuracy of
DCM is comparable to that of the traditional particle-particle particlemesh (PPPM) algorithm. The performance comparison indicates that
DCM is preferable for large systems due to the limited scaling of FFTs
within the PPPM algorithm. Moreover, the DCM requires the interface
to be identified every other MD timestep. As a consequence, this thesis
also presents an interface detection method which is (1) applicable in
real time; (2) parallelizable; and (3) scales linearly with respect to the
number of particles.

3

Acknowledgments
I like to express my deep gratitude to Prof. Paolo Bientinesi and Prof. Ahmed E. Ismail for their superb support throughout this thesis. Only their continued confidence in my work and their guidance
made this thesis possible. A special thanks goes to Prof. Benjamin Berkels who added the missing
image-segmentation knowledge and greatly improved the quality of this work. Moreover, I thank
Dr. Edoardo Di Napoli for granting me access to the JUQUEEN supercomputer at Forschungszentrum Jülich. Last but not least I am grateful for the valuable and enlightening discussions with
Dipl.-Ing. Rolf Isele-Holder and Dipl.-Ing. Daniel Tameling.

4

Contents

Contents
1. Motivation and Contribution

6

2. Background and Related Work

9

2.1. Interface Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.2. Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.2.1. Short-Range MD calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.2.2. Adaptive Cutoff . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.2.3. Long-Range Tail-Corrections . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.2.4. Particle-Mesh Ewald Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

3. Dynamic Cutoff Method

15

3.1. Design of Cutoff Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.2. Neighbor-list Build and Force Calculation . . . . . . . . . . . . . . . . . . . . . . . . . .

18

4. Interface Detection

20

4.1. Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

4.2. Removal of Small Connected Components . . . . . . . . . . . . . . . . . . . . . . . . .

23

4.3. Compute Distance to Interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

4.3.1. Fast Sweeping Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

4.3.2. Cutoff-based Fast Sweeping Method . . . . . . . . . . . . . . . . . . . . . . . . .

27

5. Parallelization

29

5.1. Neighbor-list Build and Force calculation . . . . . . . . . . . . . . . . . . . . . . . . . .

29

5.2. Removal of Small Connected-Components . . . . . . . . . . . . . . . . . . . . . . . . .

31

5.3. Cutoff-based Fast Sweeping Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

6. Results

36

6.1. Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

6.1.1. Forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

6.1.2. Density, Surface Tension and Energy . . . . . . . . . . . . . . . . . . . . . . . . .

39

6.2. Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

6.2.1. MPI + OpenMP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

6.2.2. Scalability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

7. Future Work

44

8. Conclusion

45

A. Parallel Fast Sweeping Method

50

5

1

MOTIVATION AND CONTRIBUTION

1. Motivation and Contribution
Molecular dynamics (MD) simulations are ubiquitous in computational chemistry, materials science
and biophysics. Their ability to study the dynamics of systems with an every growing number of
particles is fueled by the exponentially increasing compute power in modern computer architectures
[1]. The last decade has seen a paradigm shift from single-core CPUs to multi-core CPUs. This
paradigm shift makes good scalability an absolute necessity if we want to continue to push the limits
of molecular dynamics simulations to larger systems and longer time scales. Hence, we developed
the dynamic cutoff method (DCM) which is especially tailored for large MD simulations that require
massively parallel supercomputers to solve them.

Figure 1: Interfacial system with 16000 particles. Particles in the high-density region (center) belong to the
liquid phase while particles in the low-density region (left and right) belong to the vapor phase.
The transition from the liquid phase to the vapor phase describes the interfaces.
Molecular dynamics simulations can be roughly broken down into short-range and long-range
calculations. While long-range algorithms deal with force contributions from particles which are
far apart from each other, short-range algorithms restrict themselves to particle-particle interactions
within a certain distance, known as the cutoff. Despite the weak r −n scaling of short-range potentials
(see Section 2 for further details) with the particle-particle distance r and n being larger than the

dimension of interest (i.e., larger than 3 in most cases), it is not sufficient to use a classical shortrange algorithm for molecular systems with interfaces [2–5]. An example of an interfacial system is
illustrated in Figure 1.

Figure 2: Interfacial system using a static cutoff. Gray area denotes the cutoff. The arrows represent the
force acting on a particle.
Figure 2 illustrates an interfacial system using a static cutoff. The problem with this approach
becomes apparent when looking at the red and blue particles. On the one hand it is reasonable to use
a short-range method which neglects all the particle-particle interaction beyond the cutoff for the red
6

1

MOTIVATION AND CONTRIBUTION

particle for two reasons: First, as the distance ri,j between particles i and j increases the short-range
forces converge to zero rather quickly (see Section 2 for further details). Second, particles within the

bulk phase are arranged more or less homogeneously. As an example, accounting for the forces of

the green particles on the red particle (green arrows) would not increase the accuracy because these
forces have the same absolute value and point in opposite directions and thus cancel each other’s
contribution. For particles close to the interface (e.g., the blue particle), on the other hand, it is no
longer reasonable to assume a homogeneous distribution of particles. As a result, using a cutoff
method introduces a non negligible error because interfacial particles are missing the attractive forces
from particles outside of the cutoff which is not canceled by particles on the opposite side of the
cutoff (simply because there are none).
Figure 3 shows the particle density along the z-axis1 for the system shown in Figure 1 with respect

to three different cutoffs. Due to the neglected attractive forces of cutoff methods, small cutoffs tend
to yield lower (i.e., less accurate) densities (see Figure 3) than larger cutoffs. While this effect is
less observable for low temperatures (see Figure 3a), the missing force contributions from particles
outside of the cutoff can lead to completely wrong results for high temperatures (see Figure 3b). To
be more precise, a cutoff of 3.0σ is not sufficient to preserve the liquid-vapor interface and yields
completely wrong results.
0.9

0.9
cutoff 3.0σ
cutoff 5.0σ
cutoff 7.0σ

cutoff 3.0σ
cutoff 5.0σ
cutoff 7.0σ

0.8

0.7

0.7

0.6

0.6
Density

Density

0.8

0.5
0.4

0.5
0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0
0

10

20

30

40

50

60

70

80

90

0

10

20

30

z

40

50

60

70

80

90

z

(a) Low temperature.

(b) High temperature.

Figure 3: Density profile for different cutoffs and temperatures. Note that low temperatures result in higher
densities than high temperatures.
One solution to this problem is to use long-range solvers such as Ewald Sums [6–9] which scale
as O(N 3/2 ), where N is the number of particles, hierarchical methods [10] and mesh-based Ewald

methods [11–15], which scale as O(N log N), or fast-multipole methods [16] as well as the multi-level
summation methods [17, 18], which scale as O(N).

Despite the good asymptotic complexity of some of these solutions, they all share a common dis-

advantage: non-local communication. This disadvantage becomes progressively more severe as the
performance gap between floating-point and memory operations keeps widening. Moreover, as the
trend towards massively parallel architectures continues, the next generation of supercomputers will
be much more affected by communication overhead than today’s systems.
With these insights in mind, we developed the dynamic cutoff method (DCM) and — as a pre1

All particles in the x-y plane are averaged.

7

1

MOTIVATION AND CONTRIBUTION

requisite — a surface detection method within molecular simulations. DCM is designed for large
interfacial systems and exhibits a strictly-local communication pattern as well as a linear-time complexity O(N).

The main idea of DCM is illustrated in Figure 4: We keep the computational demand to a minimum

and increase the accuracy where it is required. Instead of applying a static cutoff to all particles, DCM
assigns large cutoffs to particles close to an interface and small cutoffs to particles in the interior. The
DCM is close in spirit to the classical (static) cutoff method [19, 20], while it exploits the fact that large
cutoffs are only required for particles close to an interface.

Figure 4: Interfacial system using a dynamic cutoff. Gray area denotes the cutoff. The arrows represent
the force acting on a particle. The yellow arrow depicts the additional force contribution due to a
larger cutoff.

Contributions. The key contributions of this master thesis can be summarized as follows:
• A scalable, linear-time MD algorithm for interfacial systems is introduced.
• A fair performance and accuracy comparison between DCM and PPPM, a state-of-the-art MD
algorithm, is performed.

• A scalable, fast, linear-time algorithm for detecting interfaces in molecular simulations is developed. In contrast to other interface-detection methods [21–23] our algorithm is so fast that it
does not affect the performance of the MD simulation.
• A cutoff-based Fast Sweeping Method (CFSM) is developed and parallelized. This algorithm
greatly outperforms the standard Fast Sweeping Method [24].

The remainder of this thesis is organized as follows: Section 2 gives background information about
MD simulations. Section 3 describes the details about DCM. Section 4 covers the topic of surface
detection within MD. Section 5 describes the implementation details of the parallelized DCM. Section
6 compares DCM to PPPM with respect to accuracy and performance. Section 7 outlines possible
directions for future work. We conclude the thesis with Section 8.

8

2

BACKGROUND AND RELATED WORK

2. Background and Related Work
This chapter relates the work of this thesis to existing approaches and gives an overview of molecular
dynamics. Moreover, we will describe short-range MD calculations in some detail and give a highlevel overview of mesh-based Ewald solvers and identify their scalability bottleneck.

2.1. Interface Detection
The identification of truly interfacial molecules (ITIM) algorithm introduced by Partay et al. [21] and the
intrinsic sampling method by Bremse et al. [22] allow to identify an interface at an atomic level. While
these algorithms are restricted to planar interfaces, the generalized ITIM (GITIM) [23] is able to detect
arbitrary interfaces; it is a mixture of the α-shapes algorithm [25] and ITIM. Willard et al. introduced

another surface detection algorithm for arbitrary interfaces in [26]. A comparison between different
methods is presented in [27].

Our interface-detection approach is motivated by the work of Berkels et al. [28] in the field of
image segmentation. As we show in this thesis, image segmentation is closely related to interface detection in molecular simulations (see Section 4). The implementation of the segmentation algorithm
by Berkels et al. [28] is part of the open-source, finite-element library QuocMesh2 [29]. This powerful
segmentation algorithm enables us to create a (1) highly scalable, (2) linear-time interface-detection
method (3) which yields a good approximation to the interface within a few milliseconds. These are
important properties since this algorithm will be executed every other MD timestep on thousands of
CPUs. In contrast to the aforementioned algorithms, our approach is not required to represent the
exact interface on an atomic level. However, we believe that the accuracy of this algorithm is quite
high; a detailed accuracy analysis for the interface detection remains future work.

2.2. Molecular Dynamics
Molecular dynamics (MD) is a classical simulation method for many-body systems. It ranges from
the simulation of the universe, stars, down to atomistic simulations of molecules and atoms. Despite
the different length scales differ, the underlying idea of advancing the particles’ position according
to Newton’s laws of motion remains the same:
ẋi (t) = vi (t),

(1)

v̇i (t) = ai (t),

for i = 1 . . . N with xi , vi , ai ∈ R3 being the position, velocity and acceleration of particle i, respec-

tively. Furthermore, ai can be expressed as the gradient of the interaction potential Vij :
ai (t) =

2

Fi (t)
,
mi

Fi (t) = −∇i V (t),
X
V(t, x0 , x1 , . . . , xN−1 ) =
Vij (rij )
i,j<i

QuocMesh also offers a variety of image processing algorithms based on finite differences.

9

(2)

2

BACKGROUND AND RELATED WORK
where rij denotes the distance between particle i and particle j at time t.

An example of an interaction potential Vij is the gravitational potential which is of importance for

the simulation of stars:

VijG (rij ) = −G

mi mj
rij

(3)

In molecular simulations, on the other hand, the Lennard-Jones (LJ) potential (see Figure 5) is
ubiquitous. It accounts for both the short-range repulsion of the Pauli exclusion principle as well as
the long-range attraction of the van der Waals forces. Due to its fast decay with r −6 this potential is

considered to be short-ranged3 . We use the LJ potential as the potential of choice for the remainder
of this thesis. However, we want to stress that the dynamic cutoff method is also applicable to other
short-range potentials as well as the short-range part of the mesh-based Ewald solvers; an extension
of DCM to these potentials is left as future work.
3
2.5

1

V LJ (rij ) = 4ε

V

LJ

/ε

1.5

Repulsive

2

0.5
0
-0.5



σ
rij

12

−



σ
rij

6 !

(4)

Attractive

-1
-1.5
1

1.5

2
r/σ

2.5

3

3.5

Figure 5: Lennard-Jones potential. (Left) Lennard-Jones potential (normalized by ε) over the particleparticle distance r (normalized by σ ). σ and ε being parameters of the LJ potential (see Equation
4).
To advance the position of each particle from one timestep to the next, equation 1, 2 and 4 need
to be solved. The naive solution to this problem requires the interaction of each particle with every
other particle, hence, resulting in a quadratic complexity O(N 2 ). Such a solution is not feasible for
large systems consisting of thousands of particles. This problem becomes even more dominant as a
typical MD simulation requires millions of timesteps to observe the events of interest (e.g., protein
folding [30]).
Figure 6 shows a rough overview of a typical MD simulation 4 . Each MD step (i.e., one iteration
of the loop) must compute the forces acting on each particle (a). Based on these forces the particle
positions can be updated (b) and the process starts over. We will refine Figure 6 in the following
subsections.

3
4

All potentials which decay with r −x with x being larger than the dimension are considered short-range potentials.
For readability reasons we do not show the processing of the computed data.

10

2

BACKGROUND AND RELATED WORK
Particle positions

(a) Compute forces

(b) Update particles

New Particle positions

Figure 6: MD simulation overview.
2.2.1. Short-Range MD calculations
Due to the rapid falloff of the LJ potential with r −6 (see Equation 4), it is reasonable to assume a
maximal interaction distance rc beyond which the interaction potential is set to zero (i.e., V (rij ) =

0, if rij > rc ). The truncation of the potential reduces the complexity from O(N 2 ) to O(N). Instead of

calculating the forces for each particle i with all the other particles, it is now sufficient to only interact
with those particles j within the cutoff rc (i.e., rij ≤ rc ).

Hockney et al. [19] introduced the link-cell method, which is based on the idea of binning the

particles into cells of size rc × rc × rc . The binning restricts the search for all neighbors of particle i to
the 26 neighboring boxes and the box of particle i itself. Hence, results in the desired complexity of

O(N). The interested reader is referred to [31] for the description of a more sophisticated version of

this algorithm.

Verlet et al. [20] were the first to propose the use of neighbor-lists (also referred to as verlet-lists) for

each particle. The key idea behind the verlet-list is to store all particle indices j which interact with

particle i (i.e., rij ≤ rc ) in the neighbor-list of particle i; and reuse this list over multiple timesteps.

To reuse this list a skin distance rs ∈ R>0 is introduced and instead of only storing particle j with

rij ≤ rc all particles with rij ≤ rc + rs are stored. The verlet-list has to be rebuild once a particle
has moved at least

rs
2

away from its old position at the last neighbor-list build [32]. This technique

drastically decreases the number of spurious distance calculations (i.e., distance calculations with
rij > rc ) at line 6 of Listing 2.1. For instance, on average 84.5%5 of the distance calculations of the

cell-list method are spurious. Using the verlet-list with a skin distance rs = 0.1rc reduces the value to
≈ 25% 6 . The choice of the skin distance rs has a strong impact on the performance: Low values of rs

require a frequent rebuilding of the neighbor-list, while high values of rs would make the verlet-list

less efficient. Chialvo et al. explore the dependence of rs on the cutoff, temperature and density in

[32].

By combining the link-cell method with verlet-lists one bins the particles into cells and uses the ob-

tained neighborhood information to build the verlet-list. This neighbor-list is then reused throughout
the force calculation (see Listing 2.1) over multiple timesteps. The pseudocode for the neighbor-list
build is outlined in Listing 2.2. We will revisit these algorithms with respect to DCM in Section 3.2.
To account for all the force contributions within the cutoff volume V cut = 4/3πrc3 , we need to check a much bigger
volume V C L = (3rc )3 , with V cut /V C L ≈ 15.5%.
6 cut
V /V verlet with V verlet = 4/3π(rc + rs )3

5

11

2

BACKGROUND AND RELATED WORK

Algorithm 2.1 Force calculation
1:
2:
3:
4:
5:
6:
7:
8:

for all atoms i do

1:

fi ← 0

2:

for all neighbors k of i do

3:

j ← nbrsi [k]

nNbrs[i] ← 0

iBin ← getBin(i)

for all jBin ∈ neighbors(iBin) do
for all atoms j of jBin do

5:

if |rij | ≤ rc then

rij ← ri − rj

6:

f ← forceLJ(|rij |)

if |rij | ≤ rc AND i < j then

7:

fi ← fi − f × rij

nbrsi [nNbrs[i]] ← j

8:

force[j] ← force[j] + f × rij

nNbrs[i] ← nNbrs[i] + 1

9:

end if

end if

10:

11:

end for

11:

12:

force[i] ← force[i] + fi

12:

13:

for all local atoms i do

4:

rij ← ri − rj

9:
10:

Algorithm 2.2 Verlet-list build

end for

13:

end for
end for
end for

The refined controlflow of a MD simulation using verlet-lists is depicted in Figure 7.
Particle positions

Neighborlist build
required?

yes

Build Neighbor-list

no
Compute forces

Update particles

New Particle positions

Figure 7: MD simulation overview with verlet-lists.

2.2.2. Adaptive Cutoff
Two weeks before submitting the thesis we received the paper of Wang et al. [33] which is an excellent
work on error estimates for inhomogeneous systems of short-range force calculations. To our surprise
this paper already introduced the idea of dynamic cutoffs. Despite some commonalities between
their method and ours, there are also substantial differences.
Even though the idea of dynamically changing the cutoffs is the same, the approach is different.
They estimate the error for each point in space and adapt the cutoff accordingly; we detect the in12

2

BACKGROUND AND RELATED WORK

terface — the region which exhibits a large error — and compute the dynamic cutoff based on the
distance to the interface. The fundamental difference between their solution and ours is that their
solution requires Fast Fourier Transforms (FFTs) and therefore requires global communications and
an O(N log N) complexity. Our method, on the other hand, is purely based on local communication

and scales linearly with the number of particles. Moreover, while their work is focused on the error
estimates, our work focuses on large systems and scalability. As such we present performance results

for systems with up to 39.3 million particles on up to 32, 768 cores. Another noticeable difference is
that their method is only applicable to systems with periodic boundary conditions while ours does
not have such a constraint.

2.2.3. Long-Range Tail-Corrections
One approach to incorporate long-range dispersion forces to short-range calculations was introduced
by Mecke et al. [34]. They showed that tail-corrections to the short-range force calculations at runtime
yield a significant accuracy boost. Their approach is somewhat orthogonal to the solution we propose
in this thesis. However, we strongly believe that a combination of these online tail-corrections and
our solution would significantly improve the accuracy of both approaches.
2.2.4. Particle-Mesh Ewald Methods
The underlying idea of the Ewald summation algorithm [8] is to split the potential into two absolute convergent sums of which one is solved in real-space — using the short-range methods of the
previous section — while the other sum is solved in reciprocal space.
Hockney and Eastwood [13] introduced the first mesh-based Ewald method called particle-particle
particle-mesh (PPPM). The particle mesh ewald (PME) [11] and the smooth particle mesh ewald (SPME)
[12] are further examples of mesh-based Ewald methods. Deserno et al. give a detailed comparison
of these different variants [14].
In contrast to the original Ewald summation, mesh-based Ewald summations use FFTs to deal with
the reciprocal sum. This reduces the complexity from O(N 3/2 ) to O(N log N). A detailed discussion
of these algorithms goes well beyond the scope of this thesis; the interested reader is referred to the
literature [11–14] for further details.

Figure 8a shows the bad weak-scaling behaviour of the PPPM algorithm on the BlueGene/Q system at Forschungszentrum Jülich. Despite the fact that in this experiment the number of particles per
core is constant, the time to solution for PPPM (orange line) increases with the size of the systems.
By contrast, our dynamic cutoff method (purple line) on the other hand shows perfect weak scaling
behavior.
Figure 8b identifies the FFTs as the main cause for the limited scaling of PPPM. These findings are
in line with the result of Schulz et al. [35], who also identify the FFTs as the main scalability bottleneck for the mesh-based Ewald summation. The reasons for the bad scaling of the FFTs is manifold:
First, FFTs are memory-bound operations and require a lot of communication; second, FFTs exhibit
an all-to-all communication pattern which prevents them from scaling perfectly; third, FFTs exhibit a
O(N log N) complexity and therefore scale worse than O(N) methods as the total number of particles
N increases.

13

2

BACKGROUND AND RELATED WORK
short-range

FFT (long-range)
other

other (long-range)
100

DCM 3.0/8.0
PPPM

50

80

% of runtime

Time [s]

40

30

20

40

20

10

0
1024

60

0
2048

4096

8192

16384

1024

32768

2048

4096

8192

#cores

#cores

(a) Weak scaling.

(b) Breakdown of runtime.

Figure 8: (a) Weak scaling of PPPM for an interfacial system with 1200 Lennard-Jones paricles per core.
(b) Percentage of PPPM spent on short-range calculations, long-range calculations and other.
The dynamic cutoff method overcomes these bottlenecks and presents an alternative solution for
large interfacial systems using short-range potentials.

14

3

DYNAMIC CUTOFF METHOD

3. Dynamic Cutoff Method
The dynamic cutoff method (DCM) can be seen as an extension to the classical short-range calculations covered in Section 2.2.1. The basic idea of DCM is to change the cutoff of each particle as a
function of the particle’s distance to the interface. This idea is motivated by the fact that particles
close to the interface require a larger cutoff than particles farther away from the interface (see Section
1). This approach both retains the low computational demands of a small cutoff while achieving the
accuracy of a large cutoff (see Section 6.1).
The advantages of DCM being based on the cell-list algorithm are manifold:
• the DCM exhibits a strictly local communication pattern (i.e., neighbor-neighbor communication),

• it has a linear-time complexity with respect to the number of particles, and
• the computational demand is solely based on the number of particles and does not increase
with void spaces in the domain.

However, DCM also has its disadvantages: first, the since different particles can have different
cutoffs, it is possible that particle i might impose a force on particle j but not vice versa (i.e., the

particle-particle interaction is no longer symmetric). This asymmetry leads to the problem that Newton’s third law 7 is no longer applicable. To be precise, DCM has to store j in the verlet-list of i and

also to store i in the verlet-list of j (a.k.a. full verlet-list) and evaluate the force between i and j twice.
Hence, as a result increases the computational demand by a factor of two. While this sounds like
a major disadvantage, we show in Section 6.2 that it is well compensated for by the better scaling
of the DCM. Moreover, due to the preferred memory access patterns and fewer data dependencies,
Newton’s third law is also commonly neglected in GPU implementations [36, 37].
Second, load balancing becomes a challenging task because domain decomposition does not do a
good job for a domain with varying numbers of particles in each subdomain. Moreover, due to the
varying computational demand per particle (depending on the cutoff) an even particle distribution
across the processors does not necessarily result in good load balancing either. Even though the
former obstacle seems only relevant for short-range solvers, and does not have a severe impact on
mesh-based Ewald solvers, it is not. Current LAMMPS [38] implementations of mesh-based Ewald
solvers do not exploit the fact that subdomains with few particles would require less computations
that those subdomains with many particles. Hence, mesh-based Ewald solvers would run into the
same load-balancing problems once they are also able to exploit these low-density regions.
Figure 9 shows a breakdown of the DCM. The DCM, in contrast to classical short-range methods,
has to compute an individual cutoff per particle — based on the particle’s distance to the interface —
before constructing the neighbor-list. This poses the challenge of detecting the interface automatically.
Once the interface has been detected, it can be reused over hundreds of timesteps since the interface,
in a typical MD simulation, only changes very slowly.
The following sections will focus on the computational challenges associated with the DCM and
assume that an interface representation is given. A detailed discussion of the interface detection is
delayed until Section 4.
7

If a body i exerts a force f onto another body j, then j also exerts −f on i.

15

3

DYNAMIC CUTOFF METHOD
Particle positions

Neighborlist build
required?

Interface
detection
required?

yes

no

yes

Detect Interface

no

Compute forces

Compute dynamic cutoff

Update particles

Build Neighbor-list

New Particle positions

Figure 9: MD simulation overview using DCM.
The interface detection method subdivides the computational domain L ∈ R3 uniformly into a set

of Nx × Ny × Nz boxes, with Nx , Ny , Nz ∈ N respectively being the number of boxes in x, y and z

direction. For now, we assume that the 3D interface is given as a set D of box-interface distances
di,j,k :

D = {di,j,k ∈ R | 0 ≤ i < Nx , 0 ≤ j < Ny , 0 ≤ k < Nz }.

(5)

D p = {di ∈ R | 0 ≤ i < N},

(6)

Once the set of box-interface distances D is available, the set of particle-interface distances D p can

be computed:

where N is the number of particles. Algorithm 3.1 shows the pseudocode for a trilinear interpolation

which is used to compute the approximate particle-interface distances D p from the box-interface
distances D.

Algorithm 3.1 Verlet-list build
1:
2:
3:
4:
5:
6:

Figure 10: Trilinear interpolation. Source: [39].

7:
8:
9:
10:

for all atoms i do

(x, y, z) ← getLowerLeftFrontBoxIndex(i)

C0,0 ← linearInterpolation(dx,y,z , dx+1,y,z , r[i]x )

C1,0 ← linearInterpolation(dx,y+1,z , dx+1,y+1,z , r[i]x )
C0,1 ← linearInterpolation(dx,y,z+1 , dx+1,y,z+1 , r[i]x )

C1,1 ← linearInterpolation(dx,y+1,z+1 , dx+1,y+1,z+1 , r[i]x )
C0 ← linearInterpolation(C0,0 , C1,0 , r[i]y )

C1 ← linearInterpolation(C0,1 , C1,1 , r[i]y )

d[i] ← linearInterpolation(C0 , C1 , r[i]z )

end for

Line 2 of Algorithm 3.1 returns the index of the lower, left, front box responsible for particle i. It

is not sufficient to just determine the box of particle i because two particles in different quadrants of
16

3

DYNAMIC CUTOFF METHOD

the same box have to interpolate over different boxes di,j,k (see Figure 11). Based on the computed

particle-interface distance di in Line 9, we need to compute its dynamic cutoff rc [i]. This involves the
use of a cutoff function which will be discussed in the next subsection.

Figure 11: Two particles belonging to the same box might have to interpolate over different neighboring
boxes.

3.1. Design of Cutoff Functions
Linear

rcmax

Exponential

rcmin

9

7

6

5

4

Interface

Out[32]=

Cutoff rc !Σ"

8

0

2

4

6

8

10

Distance from interface d !Σ"

Figure 12: Exemplary cutoff functions for a minimum cutoff of 4.5σ and a maximum cutoff of 9.0σ .

The cutoff function calculates the cutoff of each particle with respect to the particle-interface distance di . Two examples of such functions are given in Figure 12. Measurements of the root-meansquare (RMS) error of the forces (see Section 6.1) suggest that there exists an optimal function such

that the error close to the interface is of the same order as the error within the interior. The choice of
the cutoff function influences both the runtime as well as the accuracy of the method. The effect of
the cutoff function on the RMS error of the forces can be found in Section 6.1.1.
The “exponential” cutoff function, shown in Figure 12 and used in the remainder of this thesis, is

17

3

DYNAMIC CUTOFF METHOD

defined as

where


max

 rc ,
max
rc (d) =
a ∗ 2.0α (rc −d) + b,

 min
rc ,

if d ≤ 0

if 0 ≤ d < rcmax ,

(7)

if d > rcmax

rcmax − rcmin
,
max
2α rc − 1.0
b =rcmin − a

a=

(8)

The parameter α determines the steepness of the function; the more negative α becomes, the more
step-like the function gets (e.g., α = ∞ yields a step function); we use α = −0.5. d represents the

distance to the interface; we distinguish between particles inside the high-density phase (e.g., liquid
phase) d > 0 and particles inside the low-density phase d < 0 (e.g., vapor phase). All particles at the

interface or within the low-density phase are assigned the maximum cutoff rcmax , while all particles

which are at least rcmax away from the interface will be assigned the minimum cutoff rcmin . Moreover,
particles within the liquid phase will be assigned a progressively larger cutoff as they approach the

interface. Although Equation 7 only suggests three cases, we could also imagine a function which
applies the idea of a dynamic cutoff to those particles within the low-density phase as well and
therefore require more cases.

3.2. Neighbor-list Build and Force Calculation

(a) Naive. Bin size (rs + rcmax )3 .

(b) Improved. Bin size (rs + rcmin )3 .

Figure 13: Neighbor-list stencil. The blue and red boundaries mark the volume that needs to be traversed
for the blue or red particle, respectively.
The varying cutoffs pose additional challenges with respect to an efficient implementation of the
neighbor-list build. A naive implementation which uses bins of size rs + rcmax (see Figure 13a) would

result in poor neighbor-list performance because particles with small cutoffs would still have to traverse the whole volume of (3(rs + rcmax ))3 to find all the interacting particles. This approach is es-

pecially bad for interior particles. For instance, if we use rcmin = 12 rcmax , particles within the interior
would exhibit ≈ 98.5% spurious distance calculations which is far worse than the 84.5% which we

have seen in Section 2.2.1. Given that the neighbor-list build is a memory-bound operation, this
18

3

DYNAMIC CUTOFF METHOD

naive implementation would result in the same performance as if we did not use a dynamic cutoff at
all.
The solution to this problem is the following: instead of using bins with edge length rs + rcmax ,

we use an edge length of rs + rcmin or smaller (see Figure 13b). Thus, while we have to traverse a

slightly larger stencil of neighboring boxes, the traversed volume is much smaller and results in fewer
spurious distance calculations. Figure 13 illustrates this change. Notice that even the neighbor-list

volume of interfacial particles decreases. This optimization yields a 4 − 6× speedup of the neighborlist build over the naive implementation.

The force calculation (see Listing 3.2) and 3.3) are very similar to the methods of Section 2.2.1

(compare Listing 2.1 and 2.2). The main difference is that the DCM versions use a dynamic cutoff
per particle instead of a global, static cutoff; moreover, Newton’s third law is no longer applicable.
To be precise, while computing the force between particle i and j the force contribution to particle

j is no longer used; resulting in twice as many calculations. Note that the “magic” of the improved
neighbor-list build is hidden in Line 4 of Listing 3.3.
Algorithm 3.2 DCM. Force calculation
1:
2:
3:
4:
5:

Algorithm 3.3 Verlet-list build

for all atoms i do

1:

fi ← 0

2:

for all neighbors k of i do

3:

j ← nbrsi [k]

4:

rij ← ri − rj

5:

8:

if |rij | < rc [i] then

8:

9:

end if

9:

6:
7:

6:

f ← forceLJ(|rij |)

7:

fi ← fi − f × rij

10:

end for

10:

11:

force[i] ← force[i] + fi

11:

12:

end for

12:
13:

19

for all local atoms i do
nNbrs[i] ← 0

iBin ← getBin(i)

for all jBin ∈ neighbors(iBin) do
for all atoms j of jBin do
rij ← ri − rj

if |rij | < rc [i] AND i 6= j then
nbrsi [nNbrs[i]] ← j

nNbrs[i] ← nNbrs[i] + 1

end if
end for
end for
end for

4

INTERFACE DETECTION

4. Interface Detection
Particle positions

Perform 3D binning

Segment volume

Remove small connected components

Compute distance to interface

Interface position

Figure 14: Interface detection.

Figure 15: Binning.

Figure 16: Create gray-scale image/volume.
Figure 14 shows a schematic of our interface detection method. The first step is to bin all particles
into 3-dimensional boxes bi,j,k of size bx × by × bz (see Figure 15) and treat the set B of these “density
values” as a 3-dimensional image of gray-scale values (see Figure 16),

B = {bi,j,k ∈ N | 0 ≤ i < Nx , 0 ≤ j < Ny , 0 ≤ k < Nz },

where Nx , Ny and Nz are integers computed according to Equation 10:
Nx = Lx /bx ,

Ny = Ly /by ,

Nz = Lz /bz ,

(9)

(10)

such that bx ≈ by ≈ bz ; and 0 < α ≤ bx , by , bz ≤ β. The resolution of the surface is determined by
the parameter β; smaller values of β result in a higher resolution. Since the runtime of this algorithm
is determined by Nx , Ny , and Nz , β influences the runtime as well. For obvious reasons there are
limitations on the choice of both α and β. For instance, if bx , by , or bz become too small (e.g., less
than σ ) it becomes highly unlikely that a box will have more than a single particle; moreover, many

boxes will be empty and the segmentation algorithm (see Section 4.1) will fail. If, on the other hand,

bx , by and bz become too large, the grid becomes too coarse and the interface cannot be precisely

represented.

The optimal values for bx , by , and bz depend on the average distance between two particles and
20

4

INTERFACE DETECTION

thus on the system under investigation. A good choice for these values is given by the following rule
of thumb: make bx , by and bz small enough to capture the interface to a desired accuracy and large
enough that all boxes within the liquid phase have roughly the same number of particles (i.e., they
should be large enough such that density fluctuations within the liquid phase are mitigated). For our
LJ system we used a value of bx ≈ by ≈ bz ≈ 2.6σ .

The current implementation is restricted to the detection of interfaces of systems with just two

different densities (e.g., water in the liquid and vapor phase) but an extension to multicomponent
systems is easily accommodated.
The proceeding subsections discuss the remaining tasks in Figure 14.

4.1. Segmentation
Once the set B (see Equation 9) has been computed we have to segment this image/set into a seg-

mented set S

S = {si,j,k ∈ {0, 1} | 0 ≤ i < Nx , 0 ≤ j < Ny , 0 ≤ k < Nz },

0, if b
i,j,k ≈ ρl ,
such that: si,j,k =
,
1, if b
≈
ρ
v
i,j,k

(11)
(12)

where ρl and ρv are the average liquid and vapor densities, respectively. Boxes with a density closer

to the liquid density should thus be labeled with 0, while boxes with a density closer to the vapor

density should be labeled with 1. (see Figure 17).

Figure 17: Segmented 2D imgage. Black boxes denote liquid phase; white boxes denote vapor phase.
The segmentation can be computed as the minimization of the piecewise constant Mumford-Shah
functional for two-phase segmentation [28, 40]:
L

E[V , ρv , ρl ] =

Z

VL

2

(ρ(x) − ρv ) dx +

Z

VV

(ρ(x) − ρl )2 ds + νPer(V L )

(13)

where V L , V V ∈ V denote the liquid and vapor volumes, respectively, and Per(V L ) is the perimeter

of set V L . The minimization of this functional yields the segmented volumes V L and V V .

The parameter ν ∈ R+ can be seen as a penalty for a large perimeter. It determines the smoothness

of the surface (see Figure 18). If ν is too small, the interfaces will not be properly detected (see Figure

18a). If, on the other hand, ν is too large the penalty will become too large and all boxes will belong

to the same set; hence, no interface will be detected (not shown). For our simulation and accuracy
demands, we empirically determined ν = 0.005. Even though we did not experience any problems
with this value for any of our systems, automatic selection of ν remains an open question.

The implementation of this segmentation method is available in the open-source library QuocMesh
21

4

INTERFACE DETECTION

(a) ν = 0.001.

(b) ν = 0.003.

(c) ν = 0.005.

(d) ν = 0.01.

Figure 18: Effect of ν on the smoothness of the interface. The interfaces shown represent the liquid/vapor interface of a LJ system with 16000 particles after one million timesteps at a rather high
temperature.
[29]. However, the QuocMesh implementation was not suited for periodic boundary conditions
(PBC), so we had to extend QuocMesh’s functionality. Most of the runtime of segmentation is spent
on calculating finite differences. This makes a MPI parallelization straightforward, Using domain
decomposition and the concept of ghost-cells, and will not be discussed in this thesis.
Automatic detection of density values.

To make the interface detection as automatic as possible

we need to find good estimates for ρv and ρl . Naive estimates for ρl and ρl would be ρv = mini,j,k bi,j,k

and ρl = maxi,j,k bi,j,k . However, this does not always give good results since the minimum and

maximum might have large variances. Moreover, poorly chosen ρl and ρv values can result in totally

useless results. Hence, it is desirable to detect these values automatically and with high precision.

One solution to this problem is the K-Means clustering algorithm [41, 42]. K -Means tries to mini-

mize the following objective function:

E=

e
K X
N
X
j=1 i=1

||bi − cj ||2 ,

(14)

e = Nx × Ny × Nz the total numwhere cj denotes a cluster center, K the total number of clusters and N

ber of boxes. In our case K = 2, c0 ≡ ρl and c1 ≡ ρv . Figure 19 shows the histogram over all boxes

bi,j,k of a system with 16000 particles after one million timesteps at a high temperature. Figure 19 can
be interpreted as follows: bars to the left of the plot denote low density areas (i.e., vapor phase); bars
to the right of the plot denote high density areas (i.e., liquid phase) and boxes in between represent
interface densities. The vertical, red lines show that K -means automatically detects estimates for ρl
and ρv .

22

4

INTERFACE DETECTION
350

350

v

ρ

l

ρ

300

300

250

250

200

200

#boxes

#boxes

ρ

150

ρ

l

150

100

100

50

50

0

v

0
0

2

4

6

8

10

12

14

16

18

0

2

4

6

#particle per box

8

10

12

14

16

18

#particle per box

(a) Real data

(b) Artificial data

Figure 19: Histogram and k-means results. Red lines denote the clusters found by k-means.

4.2. Removal of Small Connected Components
Segmented volume S
Connectedcomponent Labeling

Remove small
components of S
Invert S
Connectedcomponent Labeling

Remove small
components of S
Invert S
New segmented volume without small components

Figure 20: Outline for the removal of small connected-components.
Depending on the choice of ν and the accuracy of the segmentation algorithm small bubbles within

the vapor phase may be falsely identified as liquids (or vice versa). Figure 21a shows an artificial
example of such a case. To remove all small components (or bubbles) smaller than a certain volume,

we apply a connected-component labeling algorithm [43] to find these connected components and
their volumes. The outline of this algorithm is depicted in Figure 20. The two boxes labeled with
remove small components of S remove all bubbles within the liquid phase and vapor phase, respectively.
This technique could also be used to determine the total number of connected components or restrict
23

4

INTERFACE DETECTION

this number to a user-specified value.
Our 3D implementation of the connected-component labeling algorithm is based on the existing
2D implementation already available within the QuocMesh library [29]. We extended this algorithm
to 3D, incorporated PBCs, and parallelized it for distributed-memory systems (see Section 5.2). The
existing implementation (see Algorithm 4.1) within QuocMesh is itself based on [44].
Algorithm 4.1 Connected-component Labeling
1:
2:
3:
4:
5:

for 0 ≤ j < Ny do

. Loop over the image in row-major order

for 0 ≤ i < Nx do

if si,j == 0 then

nlNeighbors ← getNumLabeldNeighbors(i, j)
if nlNeighbors == 0 then

labeli,j ← getNewLabel()

6:

else if nlNeighbors == 1 then

7:

10:

labeli,j ← getLabelOfNeighbor()

11:

markNeighboringLabelsAsEquivalent()

8:

else

9:

labeli,j ← getLabelOfNeighbor()

end if

12:
13:
14:

. Only consider pixels of one phase

end if
end for

15:

end for

16:

mergeAllEquivalentLabels()

For better readability we restrict the description of the labeling algorithm to the two dimensional
case. However, a 3D-extension is straightforward. A high-level overview of the connected-component
labeling procedure is shown in Algorithm 4.1. The algorithm loops over all pixels si,j (the output of

the segmentation algorithm) in row-major order and assigns each pixel a label according to the following rules:

• If a pixel has no neighbor with a label, a new label will be assigned.
• If a pixel has exactly one neighbor, the label of the neighbor will be assigned.
• In case multiple neighbors have a label, the smallest label will be assigned and all neighboring
labels will be marked as equivalent.

The details hidden in Line 11 and 16 will not be discussed further because these are technical details not required to understand for the parallel implementation discussed in Section 5.2. Figure 21
illustrates the steps necessary to remove all small bubbles for an artificial segmented image.

4.3. Compute Distance to Interface
Once we have removed all small components in S, we need to compute the distances di,j,k to the
interface of the vapor phase set S v ⊆ S and the liquid phase set S l ⊆ S with S v ∩ S l = ∅. Since our
24

4

INTERFACE DETECTION

(a) Initial Segmentation.

(b) Labeling.

(c) Merge equivalent labels.

(d) Remove small components.

(e) Invert.

(f) Labeling.

(g) Merge equivalent labels.

(h) Remove small components.

(i) Final Segmentation.

Figure 21: Removal of small connected-components under periodic boundary conditions.
current implementation only applies the dynamic cutoff method to particles within the liquid phase
(compare Equation 7) we can set di,j,k = 0 (see Equation 6) when si,j,k = 1 holds. The remaining distances can be computed by solving the signed distance function which satisfies the Eikonal equation
[24]:

|∇d(x)| = 1, x ∈ Rn ,

d(x) = 0, x ∈ Γ ⊂ Rn

(15)

where d(x) and Γ represent the distance function and the interface, respectively.
e log N)
e time with respect to the
The fast marching method (FMM) [45] solves this equation in O(N
e Zhao et al. introduced an alternative solution to the Eikonal equation
total amount of grid points N.
called the fast sweeping method (FSM) [24] as well as provided a detailed analysis of this algorithm,
e and is our preferred choice
including error bounds [46]. The FSM has a linear-time complexity O(N)
for the dynamic cutoff method. The remainder of this section outlines the original FSM as described
in [46] and then introduces our own version which makes some improvements.

25

4 INTERFACE DETECTION
5 PARALLELIZATION
4.3.1. Fast Sweeping Method
� and is our preferred choice for
be found in [48]. It is shown to have a linear-time complexity �(N)
our fast
DC method.
The
remainder
of this based
sectionon
will
outline schemes
the original
as it is described
in [48]
The
sweeping
method
is mainly
upwind
andFSM
a Gauss-Seidel
solver
using
n
and
then
introduce
our
own
version
which
makes
some
improvements
to
it.
different sweeping directions [24]. For n dimensions, Zhao et al. showed that 2 Gauss-Seidel itera-

tions with alternating sweeping orderings are enough to give a good approximation to the distance
4.3.1 Fast Sweeping Method
function for arbitrary data sets [46].
4 INTERFACE DETECTION
The
fast sweeping
method
mainly based on
upwind
schemes
and a Gauss-Seidel
solver using
Algorithm
4.2 outlines
theisimplementation
of FSM
for two
dimensions
(a 3D implementation
with
multiple sweeping directions [46]. Zhao et al. shown that 2� (� being the dimension) Gauss-Seidel
PBC is straightforward). D represents the grid of distances, and�up-right, up-left, down-right and
iterations
with
alternating
sweeping
are enough
to give
a good
the disbe
found in
[48].
It is shown
to haveorderings
a linear-time
complexity
�(N)
and approximation
is our preferredtochoice
for
down-left
denote
the
sweeping
directions.
tance
function
forThe
arbitrary
data sets
[48].
our
DC
method.
remainder
of this
section will outline the original FSM as it is described in [48]
Algorithm
4.2 outlines
theversion
implementation
for two
dimensions
(a 3Dto
implementation
with PBC
The
discretization
of Equation
15 using
Godunov
upwind
scheme
results
in the following
nonand
then
introduce
our
own
whichamakes
some
improvements
it.
is
straight-forward).
D
represents
the
grid
of
distances
up-right,
up-left,
down-right
and
down-left
linear equation [47]:
denote the sweeping directions. A comprehensive explanation of this algorithm can be found in [48];
4.3.1
Fast Sweeping Method
we will content ourself at this point with an+example
which illustrates the concepts of this algorithm.
[(di,j − dxmin ) ]2 + [(di,j − dymin )+ ]2 = h2 ,
Figure
.
The
fast17sweeping
method is mainly
based on upwind schemes and a Gauss-Seidel solver using
(16)
≤ i < that
Nx , 02�≤(�j <
Ny , the dimension) Gauss-Seidel
multiple sweeping directions [46]. Zhao et al. 0shown
being
iterations with alternating sweeping orderings are enough
to give a good approximation
up-left
up-right to the diswhere
h is the grid
spacing,data sets [48].
tance function
for arbitrary
Algorithm 4.2 outlines the implementation of FSM for two dimensions (a 3D implementation with
PBC is straight-forward).dxmin
D represents
the,grid
up-right,
up-left, down-right and down= min(di−1,j
di+1,jof), distances
dymin = min(d
(17)
i,j−1 , di,j+1 ),
down-right
down-left can be found in
left denote the sweeping directions. A comprehensive explanation
of this algorithm
[48]; we will content ourself at this point with an example. Figure 17 illustrates the correctness of
and
(
the fast sweeping method on the example of a circle.
x, Let
x >the
0 grey and white area respectively denote
(x)+
= to have the proper
(18)
the liquid and vapor phase. In oder for all
boxes
distances
we need four sweeps
up-right
up-left
0, x ≤values
0
in total. A sweep in up-right direction yields the correct
for the upper-right vapor phase and
the lower-left
liquid
phase (see
Figure
After grid
threeDadditional
sweeps
in down-left,
up-left
and
The
Gauss-Seidel
iterations
sweep
over 17).
the entire
in alternating
directions
and solve
the Equadown-right direction all distances will be computed correctly.
tion 16 for each grid point. The solution of Equation 16down-left
di,j is then used to update
each grid point
down-right

according to

new
old
up-left
up-right
di,j
= min(di,j
, di,j ).
(19)
Figure
17:
Trilinear
interpolation.
CC
[40].
Algorithm 4.2 2D Fast Sweeping Method
Figure
22 illustrates the correctness of the fast sweeping method using a circle as an example. For all
1: GaussSeidel( D, N , N , up-right)
�

�

down-right

down-left

boxes
have the proper
we need a total of four sweeps. An up-right GS sweep yields the
2: to
GaussSeidel(
D, N� distances,
, N� , up-left)
3: GaussSeidel(
N� , Ngrid
� , down-right)
correct
distances forD,those
points marked with up-right (see Figure 22). After three additional
up-right
up-left
4: GaussSeidel( D, N� , N� , down-left)
sweeps—up-left, down-right, and down-left— all distances will be computed correctly.
down-left

Algorithm
4.2 2D Method
Fast Sweeping Method
4.3.2
Fast Sweeping

down-right

up-left
up-right
1: for all
≤ � an
< arbitrary
N� do set of discrete points Γ = {�� }M , the numerical
Theorem
4.1.0 For
solution �� (�� Γ) by the fast
�=1
Figure 17: 2D FSM for a circle [48].
� sweeps,
2:
for all 0after
≤ �2<
N� do satisfies
sweeping
method
3:
if ���� == 1 then
down-right down-left
�
4:
���� ← 0
¯ � (�� Γ) ≤� �Vapor
�
(�� Γ)phase
≤ �(�� Γ) + O(|� log �|)�
5:
else
6:
���� ← inf
� Liquid phase
up-right
up-left
5 Parallelization
7:
end if
8:
end for
spatial decomposition which already exists within LAMMPS
9: end for
down-left
down-right
10: GaussSeidel( D, N� , N� , up-right)
11: GaussSeidel( D, N� , N� , up-left)
Figure 22: 2D FSM for a circle [46].
12: GaussSeidel( D, N� , N� , down-right)
13: GaussSeidel( D, N� , N� , down-left)

4.3.2 Fast Sweeping Method

20
26

�
Theorem 4.1. For an arbitrary set of discrete points Γ = {�� }M
�=1 , the numerical solution � (�� Γ) by the fast
�
sweeping method after 2 sweeps, satisfies

4

INTERFACE DETECTION

4.3.2. Cutoff-based Fast Sweeping Method
Problems arising from a parallel implementation of FSM (see Section 5.3) and the special needs of our
distance calculations led us to develop a cutoff-based fast sweeping method (CFSM). The fundamental
ideas behind CFSM are the following:
• Restrict distance calculations to those boxes close to the interface (i.e., di,j,k < rcmax + h).
• Start computation of the distances at the interface.

The former idea is specially tailored for the DCM, which does not require particle-interface distances
beyond a certain distance. The latter idea is more general and is also applicable to the original FSM.
Algorithm 4.3 describes the implementation of CFSM.

Algorithm 4.3 Cutoff-based fast sweeping method.
1:
2:
3:
4:
5:
6:

initialize(D p , Q)

for 0 ≤ iter < iter

max

do

for all boxes (x, y, z) ∈ Q do

dnew ← solveEikonal(D p , (x, y, z))
if

dnew
if

7:

≤

grid
rc

|(D p )x,y,z

then

− dnew | > ∆e then

end if

9:

(D p )x,y,z ← dnew

e (x, y, z))
addNeighborsToQueue(Q,

10:
12:

end if
end for

16:

e
swapQueues(Q, Q)
e
emptyQueue(Q)

17:

end if

13:
14:
15:

18:

. local FSM

∆e ← |(D p )x,y,z − dnew |

8:

11:

. Add interfacial boxes to Q

if ∆e < ε then
break

end for

Figure 23 illustrates the propagation of the interface. Red dots mark those grid-points just updated,
while non-updated grid-points remain black.
Performance.

Our numerical experiments show that CFSM and FSM give the same results for

grid points within the cutoff (data not shown). Moreover, as is evident in Figure 24 CFSM is preferable over FSM for various inputs. The first six inputs of Figure 24 denote artificial examples with
randomly placed points (i.e., di,j,k = 0) while the remaining two inputs are much more similar to
a typical input from an MD simulation (i.e., consisting of contiguous regions). The performance
of CFSM increased by roughly an order of magnitude over the original fast sweeping method for
relevant MD inputs.
27

4

INTERFACE DETECTION

(b) 1st iteration.

(a) Initialization.

(c) 2nd iteration.

Speedup

Figure 23: Cutoff-based FSM illustrated. Dashed line denotes the cutoff.

15
10
5
0

30_0.05

30_0.5

30_0.8

50_0.05

50_0.5

50_0.8

slab_30

slab_40

Figure 24: Speedup of CFSM over FSM for various examples using a cutoff of 4h. Examples are encoded as
follows: x_y, x denotes a grid of size x 3 and y denotes the relative amount of randomly placed
“interfacial” grid-points. The last two examples slab_x represent a domain of size x 3 with a
continuous slab of width 2x .

28

5

PARALLELIZATION

5. Parallelization
Our parallel DCM implementation is based on spatial decomposition as this is the main parallelization technique used in LAMMPS [38]. This means that each process is responsible for a certain subdomain and all its local particles. Hence, each process only has a memory requirement of O(N/P),

where N and P are the number of particles and processes, respectively.

Figure 25 outlines the distributed-memory implementation of the dynamic cutoff method; red

boxes mark operations which require communication. DCM is designed to be highly scalable and
to avoid non-local communication wherever possible.

Particle positions

Neighborlist build
required?

Perform 3D binning

Detect initial density values
yes

Migrate particles

Segment volume

no
Interface
detection
required?

Update particle positions

yes

Remove small connected components

no

Compute Forces

Compute distance to
interface for each particle

Compute distance to
interface for each bin

Compute dynamic cutoff

Interface Detection

Advance particles

New Particle positions
Build Neighbor-list

Figure 25: Breakdown of single DCM step.
The parallelization of the binning procedure, K -means and the segmentation algorithm will not be

discussed in this thesis. Instead, we focus on the parallel 3D connected-component algorithm and
the parallel CFSM.

We have also added a shared-memory parallelization based on OpenMP [48] in addition to our
MPI implementation of DCM. This allows us to start multiple threads per MPI rank which reduces
the memory requirements as well as the communication overhead and results in better performance
(see Section 6.2.1).

5.1. Neighbor-list Build and Force calculation
The OpenMP implementation of the neighbor-list build and force evaluation is based on the USEROMP package of LAMMPS. The shared-memory parallelization of the kernels for the force evaluation
29

5

PARALLELIZATION

and the Verlet-list build (see Algorithms 5.1 and 5.2) only requires a single OpenMP directive in front
of the outer-most loops8 .
Algorithm 5.1 OpenMP DCM force calculation.

Algorithm 5.2 OpenMP DCM Verlet-list build.

1:

#pragma omp for schedule(dynamic,20)

1:

#pragma omp for schedule(dynamic,20)

2:

for all atoms i do

2:

for all local atoms i do

3:
4:
5:
6:

fi ← 0

3:

for all neighbors k of i do

4:

j ← nbrsi [k]

5:

rij ← ri − rj

6:

if |rij | < rc [i] then

7:

7:

f ← forceLJ(|rij |)

8:

for all atoms j of jBin do
rij ← ri − rj

10:

end if

11:

end for

11:

12:

force[i] ← force[i] + fi

12:

13:

for all jBin ∈ neighbors(iBin) do
nbrsi [nNbrs[i]] ← j

9:

end if

10:

iBin ← getBin(i)

if |rij | < rc [i] AND i 6= j then

8:

fi ← fi − f × rij

9:

nNbrs[i] ← 0

end for

13:
14:

nNbrs[i] ← nNbrs[i] + 1

end for
end for
end for

We note here a small but important detail. Instead of using the default, static schedule, DCM
requires a dynamic schedule with a chunk-size of about 20 particles. A static schedule (i.e., evenly
distributing the particles among the threads) would result in severe load imbalancing because some
particles have larger cutoffs than others and therefore require much more computational time. The
dynamic schedule removes this problem by always assigning chunks of 20 particles at a time to
each thread. Figure 26 shows the speedups of a dynamic schedule over the static schedule for both
kernels separately. We observe that the force evaluation is more strongly affected by this change than

Speedup

the neighbor-list build.
2

Force
NL Build

1.5
1

0.5
0

2 Threads

4 Threads

8 Threads

Figure 26: Speedup of dynamic schedule over static schedule running 1 MPI process with a 2, 4 or 8 threads.
DCM with rcmin = 4σ and rcmax = 9σ . Host: Intel Xeon CPU E5-2670.

8

The actual implementation requires a few more tweaks but the concept remains the same.

30

5

PARALLELIZATION

5.2. Removal of Small Connected-Components
Our goal for a distributed-memory parallelization for the removal of small connected components
(see Section 4.2) was to design a scalable algorithm with a strictly local communication pattern. This
restriction, as we will see in this section, makes it necessary to content ourselves with an approximate
solution sufficient for our purposes. To be precise, this method gives a lower bound on the actual
label sizes (the size of each connected component) with only very little local communication.
Algorithm 5.3 gives a high-level description of the parallel labeling algorithm in two dimensions
(the 3D implementation is similar). Each process p starts to label its own subset Sp ⊂ S using the
non-parallelized labeling algorithm from section 4.2 without applying PBCs. The next step is to com-

municate the boundary labels to the left, right, top and bottom neighbors, including PBCs. During
this step equivalence classes are built and the labels are updated by summing up their sizes. Due to
the fact that we prescribe an order on the communication—first left to right, then top to bottom— information flows first horizontally, then vertically. Hence, we need an additional local communication
in the reverse direction (top-to-bottom) while using the maximum of the previously computed label
sizes. Once this is done each process can independently remove those labels smaller than a certain
threshold.

Algorithm 5.3 Parallel removal of small components in 2D.
1:
2:
3:
4:

L ← label(Sp )

sendBoundary(L, left-right)
B ← recvBoundary()

7:

L ← updateBoundary(L, B, SUM)

8:

sendBoundary(L, top-bottom)

9:

B ← recvBoundary()

5:
6:

MPI_Barrier()

12:

L ← updateBoundary(L, B, SUM)

13:

sendBoundary(L, left-right)

14:

B ← recvBoundary()

10:
11:

15:
16:
17:

MPI_Barrier()

L ← updateBoundary(L, B, MAX)

. Serial labeling
. Non-blocking

. Blocking

. Global synchronization
. Non-blocking

. Blocking

. Global synchronization
. Non-blocking

. Blocking

Sp ← removeSmallComponents(Sp , L)

Figure 27 illustrates the steps of the algorithm on an artificial example. After the first communication (see Figure 27c) neighboring processes will have detected equivalent labels and updated their
sizes according to their left and right neighbors. Figure 27d shows the labels and label sizes after
the communication with the top and bottom neighbors. Please note that this is not sufficient to compute the label sizes correctly (marked with red circles). Another communication step in the reverse
direction eliminates this problem (see Figure 27e).
31

5

PARALLELIZATION

1

1 3
3
1 3

1

3 3 3
1 2 2
5
1 1
5
5 5 5 3
3 3
1
1

(a) Initial.
3

6 5
5
6 5

12 12 12
5 5 5
12
5 5
12
12 12 12 8
8 8
3
6

1

(d) T/B communication (sum).

4

2

(c) L/R communication (sum).

6 6
6
6 6

6 6
6
6 6

12 12 12
5 5 5
12
5 5
12
12 12 12 12
12 12
3
6

2

4 4 4
3 3 3
8
2 2
8
8 8 8 8
8 8
2
1

(b) Local Labeling.
3

4

1

5 5
5
5 5

12 12 12
5 5 5
12
5 5
12
12 12 12 12
12 12
6

12

3

(e) L/R communication (max).

12

(f) Remove small components.

Figure 27: Parallel removal of small components using PBCs. All components smaller than 5 are removed.
Numbers denote the component size.
We note that there exist pathological cases for which this algorithm does not give the correct solution; Figure 28 illustrates such a case. However, these cases are highly unlikely to occur in MD
simulations. Moreover, this algorithm would still yield good lower bounds of the component sizes.
These lower bounds suffice for our purposes, since we are interested in removing small9 components.
As an example, a label spread among multiple processors (see Figure 28) is typically large enough
that our communication pattern suffices to increase the label size above the selected threshold.

P
0

P
1

P
2

P

P

3

4

P
5

P
6

P
7

Figure 28: Pathological example which illustrates that local communication does not suffice to get the
P0
correct component
size.

P1 29 shows the speedup of the parallel labeling algorithm over its sequential
Performance. Figure
counter part. Section 6.2 gives a detailed performance profile for the entire dynamic cutoff method
and puts the performance of this labeling algorithm into perspective.
9

If we were interested in large components this algorithm would not be nearly as good.

32

PARALLELIZATION

Speedup

5

7
6
5
4
3
2
1
0

40_0.2
40_0.1

2 processes

4 processes

8 processes

Figure 29: Speedup of parallel connected component labeling algorithm over serial implementation. Examples are encoded as follows: x_y, x denotes a grid of size x 3 and y denotes the relative amount of
randomly placed “wholes” (e.g., y = 0 corresponds to a single large component without “wholes”).

5.3. Cutoff-based Fast Sweeping Method

Zhao et al. [47] propose two different parallelization schemes for their fast sweeping method: (1)
Parallelize over the different sweeping directions and (2) subdivide the domain evenly among the
processes and apply multiple sweeps to each subdomain before exchanging boundary information
with neighboring processes. While (1) does not scale beyond the number of sweeping directions,
which is typically much less than the number of processes P, it also does not scale because every
e memory. On the other hand, strategy (2) yields a much more scalable solution
process requires O(N)

but might require many more floating-point operations than the serial version (see Appendix A).
Moreover, their paper shows neither speedups nor scalability results.

The remainder of this section introduces the distributed-memory parallelization of our cutoffbased fast sweeping method. This parallel implementation is based on the principle of spatial decomposition. The underlying idea is simple: each process i applies the serial CFSM to its subdomain
Dp , exchanges the boundaries after every CFSM step and terminates once the maximum change ∆e
is smaller than some threshold ε. The outline for this parallelization is shown in Algorithm 5.4.

33

5

PARALLELIZATION

Algorithm 5.4 Parallel CFSM.
1:
2:
3:
4:

initialize(D p , Q)

for 0 ≤ iter < iter

max

do

for all boxes (x, y, z) ∈ Q do

dnew ← solveEikonal(D p , (x, y, z))
grid
rc

7:

if dnew ≤

8:

end if

9:

(D p )x,y,z ← dnew

5:
6:

if

12:
13:
14:

|(D p )x,y,z

then

− dnew | > ∆e then

e (x, y, z))
addNeighborsToQueue(Q,

end if
end for

e
swapQueues(Q, Q)
e
emptyQueue(Q)

17:

∆e ← MPI_Allreduce(∆e, MAXIMUM)

18:

end if

15:
16:

19:
20:
21:

. local FSM

∆e ← |(D p )x,y,z − dnew |

10:
11:

. Add interfacial boxes to Q

if ∆e < ε then
break

ghostExchange(D p )

. Local communication

addModifiedBoundariesToQueue(D p , Q)

end for

Figure 30 illustrates the successive steps of Algorithm 5.4 on an example.
Performance.

The speedup of the parallel CFSM over the serial implementation is shown in

Figure 31. The example was run on a 16-core Intel Xeon CPU E5-2670.

34

5

PARALLELIZATION

(a) Initial.

(b) Local CFSM step.

(c) Ghost Exchange.

(d) Local CFSM step.

(e) Ghost Exchange.

(f) Local CFSM step.

Speedup

Figure 30: Parallel CFSM. Colored nodes denote active grid-points belonging to a certain processor. Black
nodes denote inactive/complete grid-points. Dashed grid-cells denote ghost cells (i.e., cells
belonging to a different process).

8
6
4
2
0

30_0.05

0.3_0.5

0.3_0.8

50_0.05

50_0.5

50_0.8

slab_30

slab_50

Figure 31: Speedup of parallel CFSM with 8 MPI processes over serial CFSM for various examples using a
cutoff of 4h. Examples are encoded as follows: x_y, x denotes a grid of size x 3 and y denotes
the relative amount of randomly placed “interfacial” grid-points. The last two examples slab_x
represent a domain of size x 3 with a continuous slab of width 2x .
35

6

RESULTS

6. Results

Figure 32: Interfacial system with 16000 Lennard-Jones particles in a domain of size 22σ × 22σ × 88σ .

This section compares the accuracy and performance of DCM with both PPPM and the classi-

cal/static cutoff method. If not otherwise mentioned, we show results for a test system of 16000 LJ
particles randomly placed in a 22σ × 22σ × 44σ box. This box is then centered in a 22σ × 22σ × 88σ

domain (see Figure 32). The parameters are scaled to LJ units with temperature T ∗ = kB T /ε = 0.7

(i.e., well below the critical temperature), timestep τ = 0.005 and distance r ∗ = r/σ . Moreover, all

simulations are run with LAMMPS version 1-Feb-14 using an NVT ensemble (i.e., keeping the number of particles, the volume and the temperature constant). All measurements are taken after an
equilibrating period of 100, 000 timesteps.

6.1. Accuracy

This section compares the accuracy of DCM against that of both PPPM and a static cutoff. Section 6.2
reuses these findings to perform a fair performance comparison.
Figure 33 illustrates that our method detects the interface at the proper location and adjusts the
cutoff accordingly.
7.5

0.9

cutoff
density

0.8

7

Cutoff rc [σ]

0.6
6
0.5
5.5
0.4
5
0.3
4.5

Particle density ρ

0.7

6.5

0.2

4

0.1

3.5

0
-40

-30

-20

-10
0
10
z-coordinate [σ]

20

30

40

Figure 33: Dynamic cutoff (rcmin = 4.0σ and rcmax = 7.0σ ) superimposed on the particle density.
6.1.1. Forces
We follow the same approach as Isele-Holder et al. [15] and measure the relative error in the forces
according to:
∆F rel

v
u
2
N  DC M
u1 X
Fi
− Fiexact
t
=
,
N
Fiexact
i=1

36

(20)

6

RESULTS

where Fiexact and FiDC M respectively represent the exact force and the force calculated by DCM for

particle i after a single timestep. We calculated Fiexact by using the Ewald summation with a large
cutoff and a large number of reciprocal vectors to get a precise estimate of the exact solution.

Figure 34 shows the root-mean-square error in the forces for linear and exponential cutoff functions

(see Section 3.1). We observe that:
1. the error in the forces decreases when increasing the minimal or maximal cutoffs;
2. increasing the minimum cutoff beyond ≈ 4.0σ hardly affects the accuracy;

3. the accuracy gain by increasing the maximum cutoff diminishes as the maximum cutoff increases; and

4. the exponential cutoff function yields a higher accuracy than the linear cutoff function.
Point (1) is somewhat obvious and shows that the algorithm works as expected: an increase in the
computational demand corresponds to an increase in the accuracy. Points (2) and (3) suggest that
there is an optimal combination of minimum and maximum cutoffs for a desired accuracy threshold ∆F rel ; we investigate this later. To explain the contribution to the error for different parameter

choices, a more elaborate analysis is required. Moreover, observation (4) is also expected because the
“exponential” cutoff function assigns larger cutoffs to more particles than the “linear” cutoff function
does and hence results in higher accuracy.
0.0008

0.001

linear
exponential

linear
exponential

0.0009

0.0007
0.0008
0.0007
0.0006

0.0005

∆ Frel

∆ Frel

0.0006

0.0004

0.0005
0.0004
0.0003

0.0003

0.0002
0.0002
0.0001
0.0001

0
4

5

6
7
max cutoff rcmax [σ]

8

9

3

(a) Minimum cutoff fixed at 4.5σ .

4

5

6
7
min cutoff rcmin [σ]

8

9

10

(b) Maximum cutoff fixed at 10σ .

Figure 34: RMS error in the forces.

We measured the error of the forces ∆Fz in the z-direction perpendicular to the interface for dif-

ferent parameter choices (see Figure 35). Figure 35a shows ∆Fz for the classical cutoff method. We

observe that the error decreases as the cutoff increases and that the error is clearly directed [49];
meaning that particles at the left interface, near z ≈ −21, are missing a force “pulling” them to the

right, while particles at the right interface, near z ≈ +21, are missing a force “pulling” them to the
left. It also highlights the main idea of the DCM since increasing the overall cutoff decreases the
error at the interface significantly, while the effect on internal particles is moderate at best. Moreover,

the error of the interfacial particles is roughly two order of magnitude larger than the error of the
internal particles (i.e., −10 ≤ z ≤ 10); this demonstrates that a static cutoff is not the ideal choice for
interfacial systems.

37

6

RESULTS
Figure 35b shows ∆Fz for DCM with varying maximum cutoffs and a linear cutoff function. It

mainly indicates three things. First, as the maximum cutoff increases, the error for interfacial parti-

cles decreases. Second, interfacial particles are still the main source of error. Finally, a linear cutoff
function is not optimal since the error increases with an increasing distance to the interface before it
drops again. The latter observation led us to design of the “exponential” cutoff function in Equation
7 that reduces the error at the interface even further (see Figure 35c). A more step-like exponential

function would reduce the error even further10 . However, the design of an analytical cutoff function
that respects the interaction potential is left for future work. The choice of the cutoff function not
only influences the accuracy but also the performance so there might be an optimal choice for this
function with respect to performance as well.
0.05
0.04

3.0/7.0
3.0/8.0
3.0/9.0

0.04

0.03

0.03

0.02

0.02

0.01

0.01
∆Fz

∆Fz

0.05

rc = 4.5
rc = 5.0
rc = 5.5

0

Dyn.

Static

Dyn.

0

-0.01

-0.01

-0.02

-0.02

-0.03

-0.03

-0.04

-0.04
-0.05

-0.05
-40

-30

-20

-10
0
10
z-coordinate [σ]

20

30

-40

40

(a) Static cutoff.

-30

-20

-10
0
10
z-coordinate [σ]

20

30

40

(b) DCM with linear cutoff function.

0.025

linear
exponential

0.02
Dyn.

Static

Dyn.

0.015
0.01

∆Fz

0.005
0
-0.005
-0.01
-0.015
-0.02
-0.025
-40

-30

-20

-10
0
10
z-coordinate [σ]

20

30

40

(c) Linear vs exponential with rcmin = 3.0σ and rcmax = 7.0σ .

Figure 35: Absolute error of the forces ∆Fz . The yellow boxes denote the interfacial region while the gray
box marks the internal region.
An important questions that still remains to be answered is how to choose the DCM parameters optimally. Figure 36 tries to answer this question. It shows the RMS error of the forces ∆F rel with respect
to both the minimal cutoff as well as the maximum cutoff. The origin of each arrow (see Figure 36)
indicates the optimal parameter choice for a desired accuracy with respect to the best performance.

For instance, to achieve an accuracy of at least ∆F rel = 0.00149, rcmin = 3.0σ and rcmax = 5.0σ would

be optimal, while rcmin = 3.0σ and rcmax = 5.5σ would be the best choice for a slightly higher accu10

This can be achieved by reducing the α value of Equation 8.

38

6

RESULTS

racy of ∆F = 0.00133. This figure validates our previous observation that the error monotonically

decreases with an increasing minimum and maximum cutoff. The more important observation, however, is indicated by the arrows. The sequence in which to pick the optimal rcmin /rcmax combination for

a given accuracy is indicated by the arrow heads. We make the following observations:

1. The minimum cutoff has a stronger effect on performance11 than the maximum cutoff. This
is shown by the fact that the arrows point towards a larger maximum cutoff before pointing
towards a larger minimum cutoff. In other words, with respect to performance it is preferable
to increase the maximum cutoff instead of the minimum cutoff to achieve higher accuracy.

2. The optimal parameter choices are far away from the diagonal (the diagonal represents a static
cutoff); hence, a dynamic cutoff is superior to a static cutoff.
3. Increasing the maximum cutoff beyond some value does not increase the total accuracy further,
because the error of the internal particles eventually dominates. At this point a larger minimum
cutoff is required.

Min cutoff [σ]

3

3

5

7

9

3

5

5

7

7

9

9

3

5

7
Max cutoff [σ]

10-3

10-4

9

Figure 36: RMS error of the forces ∆F rel . The origin of each arrow denotes that this configuration is
preferable w.r.t. performance while meeting a certain accuracy threshold. Arrow direction indicate
the sequence in which to pick optimal parameters for an increasing accuracy threshold.

6.1.2. Density, Surface Tension and Energy
The surface tension, its uncertainty and the particle density have been computed in the same way as
proposed in [15].
As it is evident from Figure 37, the DCM results for the density and surface tension converge to
the results of PPPM as we increase the maximum cutoff. The monotonic convergence of the DCM results toward the PPPM solution once again shows the proper behavior of DCM. Figure 37 reveals
11

This is especially true as the ratio between interfacial and internal particles decreases.

39

6

RESULTS

many noteworthy properties: First, DCM with rcmin = rcmax gives exactly the same results as the

static rc version. Second, the results for DCM and the static cutoff are almost identical for cutoffs

of rc ≥ 8σ , although DCM achieves these results 4 − 5 times faster. Hence, we do not trade accuracy

for performance. For instance, if one compares DCM 3.0/8.0 against static 5.0 much more accurate

results are obtained in just a fraction of the time for the static 5.0 case. Third, the DCM results with
rcmax ≥ 8.0 can be considered identical to those of the PPPM because the relative difference in the

density between DCM and PPPM is less than 1% and because the slightly larger relative difference

in the surface tension can be corrected by a tail correction (not shown). Moreover, PPPM is also an
approximation to the correct solution and itself has an intrinsic error. Fourth, the surface tension

and density monotonically increase as larger maximum cutoffs keep the particles closer together and
preserve the interface more accurately.
1.2

0.845
0.84

1.1
Surface tension γ

Particle density ρ

0.835
0.83
0.825
0.82

1

0.9

0.815
0.8

rcmin= 3.0
rcmin= 4.0
static

0.81

rcmin= 3.0
rcmin= 4.0
static

0.7

0.805
3

4

5

6
rcmax

7

8

3

9

4

5

[σ]

(a) Density.

6
rcmax [σ]

7

8

9

(b) Surface Tension.

Figure 37: Solid line denotes the PPPM results.
Table 1 shows the energy of the DCM in comparison to the results of the static cutoff method and
PPPM. Again, we see that DCM behaves as expected: as the maximum cutoff increases, we approach
the PPPM and static solutions. Moreover, we also observe that the results of DCM rcmin /rcmax are right

in between static rc = rcmin and static rc = rcmax , exactly as we would expect.
static 3.0
static 5.0
static 7.0
DCM 3.0/5.0
DCM 3.0/7.0
DCM 3.0/8.0
DCM 4.5/9.5
PPPM

Total energy
-4.415
-4.772
-4.830
-4.557
-4.636
-4.660
-4.805
-4.859

Std. deviation
0.013
0.016
0.013
0.013
0.017
0.016
0.014
0.013

Table 1: Total energy.

40

6

RESULTS

6.2. Performance
This subsection compares the performance of DCM against the statical cutoff method and PPPM. We
use a Intel Xeon E5-2670 cluster for small test cases and the JUQUEEN supercomputer at Forschungszentrum Jülich for large scalability measurements. If not otherwise mentioned, we use 16 cores (i.e.,
two Intel Xeon E5-2670 CPUs). The minimum and maximum cutoff of DCM are respectively set to
rcmin = 3.0 and rcmax ≥ 8.0, such that the DCM results match the results of PPPM.

Figure 38 shows the speedup of DCM over its static counterpart. It shows the benefit of using the

dynamic cutoff idea with respect to performance. The shown speedups range from 4.8× to 5.7×. The

Speedup

speedup increases as the ratio between interfacial to internal particle decreases.
5
4
3
2
1
0

3.0/7.0

3.0/8.0

3.5/8.0

3.0/9.0

3.5/9.0

Figure 38: Speedup of DCM (using the exponential cutoff function) over a static cutoff with rc = rcmax for
a system with 192, 000 particles. Encoding: rcmin /rcmax .
6.2.1. MPI + OpenMP
Table 2 shows the speedup of our hybrid (MPI + OpenMP) DCM implementation over our pure MPI
implementation. The total number of tasks (MPI ranks times threads) is kept constant for a varying
number of threads. The results strengthen the intuition that a hybrid implementation is at least as
good as a pure MPI implementation. In fact, the hybrid implementation increases the performance
by up to 29%.

#tasks
1024
2048
4096

2 threads
1.056
1.084
1.001

4 threads
1.280
1.121
1.016

8 threads
1.247
1.291
1.019

16 threads
1.199
1.229
1.127

Table 2: Speedup of OpenMP + MPI over a pure MPI implementation. All simulations use the same system
consisting of three million particles. Host: IBM BlueGene/Q.

6.2.2. Scalability
The scalability results of this section are conducted on the IBM BlueGene/Q system JUQUEEN at
Forschungszentrum Jülich. JUQUEEN comprises 28, 672 nodes arranged in 28 racks. Each node
consists of an IBM PowerPC A2 with 16 cores and 4-way simultaneous multi-threading (SMT). This

yields a total of 458, 752 cores (16, 384 cores per rack). Performance measurements indicate that SMT
41

6

RESULTS

yields significant speedups (data not shown). Hence, all of our tests fully exploit the SMT feature of
the IBM PowerPC running 4 tasks per core.
The domain for the weak scaling results is set up such that the volume V p =

Lx
Px

×

Ly
Py

×

Lz
Pz

and

particle density per MPI process p are kept constant with Lx , Ly , Lz , Px , Py and Pz being the com-

putational domain and the number of MPI ranks in the x, y, and z directions, respectively. The
average number of particles per core is thus kept roughly constant at 1200 particles/core. Since

LAMMPS applies domain decomposition we try to keep load balancing issues to a minimum by

avoiding void spaces as much as possible. To be precise, N = #cores · 1, 200 particles are randomly
placed into a box B of volume L × L × 2L. B itself is then centered in a surrounding box of volume

L × L × (2L + d) with d = 25.0σ such that the interface-interface distance d (due to periodic boundary

conditions) is sufficiently large. Strong scaling measurements use a system which is ideal for 32, 768

cores (i.e., 32, 768 · 1, 200 = 39, 321, 600 particles).

Before we start our scalability discussion of DCM and PPPM, we stress that the settings for PPPM

are chosen according to [15] and are considered to be optimal.
The red curve in Figure 39a shows the timings for DCM with sequential interface detection while
the blue and purple curves show those timings with our parallel implementation of the interface
detection. This shows both that parallel interface detection is clearly required and that our parallel
implementation performs well. Moreover, DCM shows perfect weak-scaling and surpasses PPPM
at around 1, 000 to 2, 000 cores. This makes DCM the preferred choice for large MD simulations on

highly parallel systems.

The ideal strong scaling of DCM is shown in Figure 39b: while PPPM slows down significantly

around 4000 to 8000 cores, DCM does not.
70

256

DCM 3.5/9.0 (serial interface)
DCM 3.0/8.0
DCM 3.5/9.0
PPPM

60

DCM 3.5/9.0
DCM 3.0/8.0
PPPM

128

50
40

Time [s]

Time [s]

64

30

32

16
20

7.7x
8

10
0
1024

2048

4096

8192

16384

32768

#cores

4
1024

2048

4096

8192

16384

32768

#cores

(a) Weak scaling.

(b) Strong scaling.

Figure 39: Scalability of DCM vs. PPPM. Results are accumulated over 200 MD steps.
Figure 40a shows the runtime of the neighbor-list build (NL) including interface detection (Interface), short-range Force calculation, communication (Comm) and remainder (Other) separately. Figure 40a highlights several interesting properties of the DCM. First, the communication time remains
more or less constant, as we would expect from an algorithm with strictly local communications. In
addition, the interface-detection contributes less than 1% of the overall runtime.
42

6

RESULTS
Figure 40b breaks the interface detection method down into its subroutines (see Section 4). We

observe that 75 − 80% of the interface detection time is due to volume segmentation. Moreover, the

time it takes to segment the volume increases slightly when going from 1024 cores to 16384 cores.
This is the source of a residual scaling issue that does not affect the overall scalability of DCM. The

remaining subroutines behave nicely and have only a minor effect on the runtime of the interface
detection routine.
NL-Build
Comm

Interface
Other

Force

Binning
CFSM

K-Means
Interpolate

% of interface-detection

100

Runtime [s]

10
5
0

Segment
Remove

80
60
40
20
0

1024 2048 4096 8192 16384 32768
#cores
(a) Breakdown DCM.

1024

2048

4096
#cores

8192

16384

(b) Breakdown Interface detection.

Figure 40: Breakdown of (a) DCM 4.5/9.5 and (b) its interface detection for a weak scaling simulation with
1200 particles per core. Results are accumulated over 200 MD steps.

43

7

FUTURE WORK

7. Future Work
As we have mentioned several times throughout this thesis, there are a number of items left to be explored in the future. This section lists some possible future directions for the dynamic cutoff method
with respect to performance, features, and accuracy.
Performance. Since DCM avoids the use of Newton’s third law (i.e., has a good memory access
pattern) and requires large cutoffs to work (i.e., it is a compute-bound algorithm), it is a perfect
candidate for co-processors such as GPUs [36, 37] or Intel’s Xeon Phi [50]. The work of Nguyen
et al. [51] indicates that the speedup due to GPUs increases as the cutoff becomes larger. Another
optimization along those lines would be to vectorize the two most important kernels —the force
calculation and the neighbor-list build—via C/C++ AVX2 intrinsics. Pennycook et al. [50] show
promising speedups from using explicit vectorization in MD simulations.
So far we have restricted ourselves to periodic boundary conditions. This, however, made it necessary to introduce a void space to separate the interfaces from each other. Hence, we want to extend
DCM to allow non-PBCs as well and apply this idea to the same system that we have studied so far.
Moreover, a careful design of the cutoff function with respect to both accuracy and performance is
left as future work.
Accuracy. Although the interface detection method reproduced the interface well for our purposes,
we would like to carry out more accuracy studies and improve the accuracy of the detected interface
if required.
Features. Moreover, this thesis applied the DCM only to a single-component Lennard-Jones system.
The next step would be to extend DCM’s capabilities to multi-component systems and apply the
dynamic cutoff idea to other short-range potentials.
Although we have applied DCM only to a system with a planar liquid-vapor interface, we believe
that DCM is also applicable to systems with non-planar and liquid-liquid interfaces as well.
As pointed out in Section 3.1 we want to extend the dynamic cutoff principle not only to particles
within the liquid phase but also to those particles within the vapor phase as well.
A possibility to increase the accuracy of DCM quite substantially is to incorporate tail correction
to the potential at runtime [34]. Given that we know the approximate volume V L and V V of the

liquid and vapor phase and their approximate particle density ρL and ρV , we could follow the ideas
of Mecke et al. [34] and compute the force correction ∆Fitail for each particle i according to:

where
∆utail
i

=

Z

rij >rc
rj ∈V L

∆Fitail = −
L

d(∆utail
i )
dr

uij (rij )ρ drij +

Z

rij >rc
rj ∈V V

uij (rij )ρV drij

(21)

(22)

at every timestep. This correction does not require particle-particle interactions and is rather inexpensive to compute. Due to the rapid decay of short-range potentials, one could restrict the evaluation
of the integral to the proximity of each particle, say 2rc , to keep the communication local.
44

8

CONCLUSION

8. Conclusion
We have developed a dynamic cutoff method for computing long-range dispersion interactions in
molecular simulations which is based on computing the distance between particles and the interface.
We also present a scalable, linear-time interface detection method for non-planar interfaces. This
interface detection method enabled us to develop the DCM for short-range potentials. DCM is specially tailored for massively parallel supercomputers simulating interfacial systems with millions of
particles.
We have implemented DCM as part of LAMMPS and showed that it exhibits desired properties
such as (1) linear-time complexity, (2) local communication, and (3) ideal weak- and strong-scaling.
Moreover, our accuracy results show that DCM is able to achieve the same accuracy as state-of-theart algorithms for interfacial Lennard-Jones systems, while outperforming them for large systems.
For instance, looking at a system with 39.3 million particles, DCM achieved the same accuracy as

PPPM but was 7.7 times faster. While more case studies are required and more DCM features are
desirable, our preliminary results indicate that DCM is a promising algorithm for massively parallel
supercomputers.

We also introduced an interface detection method which is highly scalable and fast enough to be
applicable in real time throughout the course of a MD simulation. As such, it might open the door to
a wide variety of new MD applications.

45

References

References
[1] TOP500 Supercomputing site. http://www.top500.org/, 2014.
[2] Gustavo A Chapela, Graham Saville, Stephen M Thompson, and John S Rowlinson. Computer
simulation of a gas–liquid surface. part 1. Journal of the Chemical Society, Faraday Transactions 2:
Molecular and Chemical Physics, 73(7):1133–1144, 1977.
[3] EM Blokhuis, D Bedeaux, CD Holcomb, and JA Zollweg. Tail corrections to the surface tension
of a lennard-jones liquid-vapour interface. Molecular Physics, 85(3):665–669, 1995.
[4] Ahmed E Ismail, Mesfin Tsige, Pieter J Veld In’t, and Gary S Grest. Surface tension of normal
and branched alkanes. Molecular Physics, 105(23-24):3155–3163, 2007.
[5] Mingxue Guo, D-Y Peng, and Benjamin C-Y Lu. On the long-range corrections to computer
simulation results for the lennard-jones vapor-liquid interface. Fluid phase equilibria, 130(1):19–
30, 1997.
[6] P. P. Ewald. Die berechnung optischer und elektrostatischer gitterpotentiale. Annalen der Physik,
369(3):253–287, 1921.
[7] DE Williams. Accelerated convergence of crystal-lattice potential sums. Acta Crystallographica
Section A: Crystal Physics, Diffraction, Theoretical and General Crystallography, 27(5):452–455, 1971.
[8] John W Perram, Henrik G Petersen, and Simon W De Leeuw. An algorithm for the simulation of
condensed matter which grows as the 3/2 power of the number of particles. Molecular Physics,
65(4):875–893, 1988.
[9] Naoki Karasawa and William A Goddard III. Acceleration of convergence for lattice sums. The
Journal of Physical Chemistry, 93(21):7320–7327, 1989.
[10] Josh Barnes and Piet Hut. A hierarchical o (n log n) force-calculation algorithm. 1986.
[11] Tom Darden, Darrin York, and Lee Pedersen. Particle mesh ewald: An n log(n) method for
ewald sums in large systems. The Journal of chemical physics, 98(12):10089–10092, 1993.
[12] Ulrich Essmann, Lalith Perera, Max L Berkowitz, Tom Darden, Hsing Lee, and Lee G Pedersen.
A smooth particle mesh ewald method. The Journal of Chemical Physics, 103(19):8577–8593, 1995.
[13] R W Hockney and J W Eastwood. Computer Simulation Using Particles - R.W Hockney, J.W
Eastwood - Google Books. 1988.
[14] Markus Deserno and Christian Holm. How to mesh up ewald sums. i. a theoretical and numerical comparison of various particle mesh routines. The Journal of chemical physics, 109(18):7678–
7693, 1998.
[15] Rolf E Isele-Holder, Wayne Mitchell, and Ahmed E Ismail. Development and application of a
particle-particle particle-mesh ewald method for dispersion interactions. The Journal of chemical
physics, 137(17):174107, 2012.
46

References
[16] Leslie Greengard and Vladimir Rokhlin. A fast algorithm for particle simulations. Journal of
computational physics, 73(2):325–348, 1987.
[17] David J Hardy, John E Stone, and Klaus Schulten. Multilevel summation of electrostatic potentials using graphics processing units. Parallel Computing, 35(3):164–177, 2009.
[18] Daniel Tameling, Paul Springer, Paolo Bientinesi, and Ahmed E Ismail. Multilevel summation for dispersion: A linear-time algorithm for r- 6 potentials. The Journal of chemical physics,
140(2):024105, 2014.
[19] RW Hockney, SP Goel, and JW Eastwood. Quiet high-resolution computer models of a plasma.
Journal of Computational Physics, 14(2):148–158, 1974.
[20] Loup Verlet. Computer" experiments" on classical fluids. i. thermodynamical properties of
lennard-jones molecules. Physical review, 159(1):98, 1967.
[21] Lívia B Pártay, György Hantal, Pál Jedlovszky, Árpád Vincze, and George Horvai. A new
method for determining the interfacial molecules and characterizing the surface roughness in
computer simulations. application to the liquid–vapor interface of water. Journal of computational chemistry, 29(6):945–956, 2008.
[22] Fernando Bresme, Enrique Chacón, and Pedro Tarazona. Molecular dynamics investigation
of the intrinsic structure of water–fluid interfaces via the intrinsic sampling method. Physical
Chemistry Chemical Physics, 10(32):4704–4715, 2008.
[23] Marcello Sega, Sofia S Kantorovich, Pál Jedlovszky, and Miguel Jorge. The generalized identification of truly interfacial molecules (itim) algorithm for nonplanar interfaces. The Journal of
chemical physics, 138(4):044110, 2013.
[24] Hong-Kai Zhao, Stanley Osher, Barry Merriman, and Myungjoo Kang. Implicit and Nonparametric Shape Reconstruction from Unorganized Data Using a Variational Level Set Method.
Computer Vision and Image Understanding, 80(3):295–314, December 2000.
[25] Herbert Edelsbrunner and Ernst P Mücke. Three-dimensional alpha shapes. ACM Transactions
on Graphics (TOG), 13(1):43–72, 1994.
[26] Adam P Willard and David Chandler. Instantaneous liquid interfaces. The Journal of Physical
Chemistry B, 114(5):1954–1958, 2010.
[27] Miguel Jorge, Pál Jedlovszky, and M Natália DS Cordeiro. A critical assessment of methods
for the intrinsic analysis of liquid interfaces. 1. surface site distributions. The Journal of Physical
Chemistry C, 114(25):11169–11179, 2010.
[28] Benjamin Berkels. An unconstrained multiphase thresholding approach for image segmentation. In Proceedings of the Second International Conference on Scale Space Methods and Variational
Methods in Computer Vision (SSVM 2009), volume 5567 of Lecture Notes in Computer Science, pages
26–37. Springer, 2009.
47

References
[29] AG Rumpf, Institute for Numerical Simulation, University of Bonn. Quocmesh software library.
http://numod.ins.uni-bonn.de/software/quocmesh/.
[30] Martin Karplus and J Andrew McCammon. Molecular dynamics simulations of biomolecules.
Nature Structural & Molecular Biology, 9(9):646–652, 2002.
[31] Pedro Gonnet. A simple algorithm to accelerate the computation of non-bonded interactions
in cell-based molecular dynamics simulations. Journal of Computational Chemistry, 28(2):570–573,
2006.
[32] Ariel A Chialvo and Pablo G Debenedetti. On the use of the verlet neighbor list in molecular
dynamics. Computer Physics Communications, 60(2):215–224, 1990.
[33] Han Wang, Christof Schütte, and Pingwen Zhang. Error estimate of short-range force calculation
in inhomogeneous molecular systems. Physical Review E, 86(2):026704, August 2012.
[34] Matthias Mecke, Jochen Winkelmann, and Johann Fischer. Molecular dynamics simulation of
the liquid–vapor interface: The lennard-jones fluid. The Journal of chemical physics, 107(21):9264–
9270, 1997.
[35] Roland Schulz, Benjamin Lindner, Loukas Petridis, and Jeremy C Smith. Scaling of multimillionatom biological molecular dynamics simulation on a petascale supercomputer. Journal of Chemical Theory and Computation, 5(10):2798–2808, 2009.
[36] Joshua A Anderson, Chris D Lorenz, and Alex Travesset. General purpose molecular dynamics
simulations fully implemented on graphics processing units. Journal of Computational Physics,
227(10):5342–5359, 2008.
[37] W Michael Brown, Peng Wang, Steven J Plimpton, and Arnold N Tharrington. Implementing molecular dynamics on hybrid high performance computers–short range forces. Computer
Physics Communications, 182(4):898–911, 2011.
[38] Steve Plimpton. Fast parallel algorithms for short-range molecular dynamics. Journal of computational physics, 117(1):1–19, 1995.
[39] Marmelad. https://en.wikipedia.org/wiki/File:3D_interpolation2.svg.
[40] David Mumford and Jayant Shah. Optimal approximations by piecewise smooth functions and
associated variational problems. Communications on pure and applied mathematics, 42(5):577–685,
1989.
[41] John A Hartigan. Clustering algorithms. John Wiley & Sons, Inc., 1975.
[42] John A Hartigan and Manchek A Wong. Algorithm as 136: A k-means clustering algorithm.
Applied statistics, pages 100–108, 1979.
[43] Hanan Samet and Markku Tamminen. Efficient component labeling of images of arbitrary dimension represented by linear bintrees. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 10(4):579–586, 1988.
48

References
[44] Connected Components Labeling. http://homepages.inf.ed.ac.uk/rbf/HIPR2/label.htm, 2014.
[45] James A Sethian. A fast marching level set method for monotonically advancing fronts. Proceedings of the National Academy of Sciences, 93(4):1591–1595, 1996.
[46] Hongkai Zhao. A fast sweeping method for eikonal equations. Mathematics of computation,
74(250):603–627, 2005.
[47] Hongkai Zhao. Parallel implementations of the fast sweeping method. Journal of Computational
Mathematics, 25(4), 2007.
[48] OpenMP Architecture Review Board. OpenMP application program interface version 4.0, May
2014.
[49] Rolf E Isele-Holder, Wayne Mitchell, Jeff R Hammond, Axel Kohlmeyer, and Ahmed E Ismail.
Reconsidering Dispersion Potentials: Reduced Cutoffs in Mesh-Based Ewald Solvers Can Be
Faster Than Truncation. Journal of Chemical Theory and Computation, 9(12):5412–5420, December
2013.
[50] Simon J Pennycook, Chris J Hughes, M Smelyanskiy, and SA Jarvis. Exploring simd for molecular dynamics, using intel R xeon R processors and intel R xeon phi coprocessors. In Parallel
& Distributed Processing (IPDPS), 2013 IEEE 27th International Symposium on, pages 1085–1097.
IEEE, 2013.
[51] Trung Dac Nguyen, Jan-Michael Y Carrillo, Andrey V Dobrynin, and W Michael Brown. A case
study of truncated electrostatics for simulation of polyelectrolyte brushes on gpu accelerators.
Journal of Chemical Theory and Computation, 9(1):73–83, 2012.

49

A PARALLEL FAST SWEEPING METHOD

A. Parallel Fast Sweeping Method
We like to point out that the parallel implementation of the fast sweeping method by Zhao et al. [47]
potentially results in many more floating-point operations than the serial implementation would
require.
As an example, let D ∈ Rn×n be the two dimensional domain of interest, Di,j ∈ Rn/2×n/2 , Di,j ⊂ D

are the subdomains of four processes pi,j , 0 ≤ i, j < 2. Let N = n × n be the total number of grid

points. Imagine that we place a single point into the center of D0,0 . All processes would start to
process their subdomain requiring ≈ N/4 operations each (i.e., N operations in total). However, all
processes except for p0,0 would do useless work. At the end of the first computational step only

D0,0 has been computed correctly. The next step is to exchange the solution at the boundary of

each subdomain with its neighbors. Once the communication is done, each process again starts to
process its subdomain again requiring ≈ N/4 operations each. At the end of this second computation

phase D0,1 D1,0 have also been computed correctly. However, to compute D1,1 correctly another

computation step is required. All in all this algorithm took roughly three times more operations than
the serial implementation and as a consequence results in poor scaling.

50

