Raising the Performance of the Tinker-HP Molecular Modeling Package [Article v1.0]

TL;DR

This paper reviews and demonstrates how to optimize the Tinker-HP molecular dynamics package for modern Intel CPU architectures, achieving significant performance improvements through vectorization and parallelization strategies.

Contribution

The paper introduces a comprehensive optimization strategy for Tinker-HP on Intel CPUs, focusing on AVX-512 vectorization and parallel scaling, with practical performance results.

Findings

Performance gains from AVX-512 optimization

Effective vectorization of key code hotspots

Enhanced scalability on petascale architectures

Abstract

This living paper reviews the present High Performance Computing (HPC) capabilities of the Tinker-HP molecular modeling package. We focus here on the reference, double precision, massively parallel molecular dynamics engine present in Tinker-HP and dedicated to perform large scale simulations. We show how it can be adapted to recent Intel Central Processing Unit (CPU) petascale architectures. First, we discuss the new set of Intel Advanced Vector Extensions 512 (Intel AVX-512) instructions present in recent Intel processors (e.g., the Intel Xeon Scalable and Intel Xeon Phi 2nd generation processors) allowing for larger vectorization enhancements. These instructions constitute the central source of potential computational gains when using the latest processors, justifying important vectorization efforts for developers. We then briefly review the organization of the Tinker-HP code and identify the computational hotspots which require Intel AVX-512 optimization and we propose a general and optimal strategy to vectorize those particular parts of the code. We intended to present our optimization strategy in a pedagogical way so it could benefit to other researchers and students interested in gaining performances in their own software. Finally we present the performance enhancements obtained compared to the unoptimized code both sequentially and at the scaling limit in parallel for classical non-polarizable (CHARMM) and polarizable force fields (AMOEBA). This paper never ceases to be updated as we accumulate new data on the associated Github repository between new versions of this living paper.

Tables

Table 1. Table 1: MS used for the performance measurements. The numbers of cores are taken from  [ 1 ] for comparison. The CPU2 raw gives the number of cores which produced the best performance (See tables  5 , 6 and 7 ). For the sequential performance measures, only one core was used.

MS	Ubiquitin	DHFR	Puddle	COX-2
Atoms	$\mathrm{9\; 737}$	$\mathrm{23\; 558}$	$\mathrm{96\; 000}$	$\mathrm{174\; 219}$
CPU	$480$	$680$	$\mathrm{1\; 440}$	$\mathrm{2\; 400}$
CPU2		$960$		$\mathrm{3\; 000}$
MS	Pond	Lake	STMV	Ribosome
Atoms	$\mathrm{288\; 000}$	$\mathrm{864\; 000}$	$\mathrm{1\; 066\; 628}$	$\mathrm{3\; 484\; 755}$
CPU	$\mathrm{2\; 400}$	$\mathrm{7\; 200}$	$\mathrm{10\; 800}$	$\mathrm{10\; 800}$
CPU2			$\mathrm{16\; 200}$	$\mathrm{16\; 200}$

Table 2. Table 2: Profiling of Rel using Intel VTune Amplifier. Simulations ran on one core and 100 steps. MS is DHFR with AMOEBA polarizable force field and with CHARMM force field (no polarization). Most important NUC and computational hostspots are shown in separate frames. vmlinux is the system kernel, performing memory operations and system calls. For CHARMM calculation, image is splitted in two parts. The vectorized routines will use image(2). So, only the starred lines are counted in the total CPU time for comparison with Vec . The † † {\dagger} on some lines indicate routines that have not been vectorized in Vec . Thus, they don’t count in the total CPU time for comparison.

Module	CPU Time	Vector
	(s)	usage %
NUC hotspots
Total CPU time : 36.0896 s
vmlinux	$27.5005$	$100.00$
libmkl_avx512.so	$5.7625$	$100.00$
libmpi.so.12.0	$2.7144$	$0.00$
libc-2.17.so	$0.0862$	$0.00$
libmkl_intel_lp64.so	$0.0120$	$0.00$
libiompstubs5.so	$0.0090$	$0.00$
libmpifort.so.12.0	$0.0050$	$0.00$
DHFR (AMOEBA, polarizable model)
Computational hotspots
Total CPU time : 278.9512 s (100 steps)
tmatxb_pme2	$100.9210$	$0.00$
epolar1	$52.7085$	$0.00$
ehal1	$52.4679$	$0.00$
empole1	$28.9127$	$0.00$
image	$25.2141$	$0.00$
efld0_direct2	$17.4910$	$0.00$
mlistcell†	$8.2543$	$0.00$
vlistcell†	$7.2870$	$0.00$
torque	$2.1355$	$0.00$
DHFR (CHARMM, non-polarizable model)
Computational hotspots
Total CPU time : 24.7982∗ s (100 steps)
elj1∗	$15.3259$	$0.00$
echarge1∗	$6.7309$	$0.00$
image (1)	$3.4130$	$0.00$
clistcell†	$2.8517$	$0.00$
image (2)∗	$2.7414$	$0.00$

Table 3. Table 3: Profiling of Vec using Intel VTune Amplifier. Simulations ran on one core and 100 steps. MS is DHFR with AMOEBA polarizable force field and with CHARMM force field (no polarization). For the vectorized routines, the Vector Usage percentages go from 63.9 to 100%. Only the starred lines are counted in the total CPU time for DHFR with CHARMM.

Module	CPU Time	Vector	Boost
	(s)	Usage %	factor
NUC hotspots
Total CPU time : 37.3438 s
vmlinux	$25.4116$	$100.00$
libmkl_avx512.so	$6.1404$	$100.00$
libmpi.so.12.0	$2.7094$	$0.00$
libmkl_vml_avx512.so	$2.6733$	$100.00$
libc-2.17.so	$0.0703$	$0.00$
libmkl_intel_lp64.so	$0.3208$	$0.00$
libmpifort.so.12.0	$0.0110$	$0.00$
libiompstubs5.so	$0.0070$	$0.00$
DHFR (AMOEBA, polarizable)
Computational hotspots			$2.0410$
Total CPU time : 136.4569 s (100 steps)
tmatxb_pme2vec	$62.9675$	$63.90$	$1.6027$
epolar1vec	$29.2285$	$94.90$	$1.8033$
ehal1vec	$19.9918$	$67.90$	$2.6245$
empole1vec	$11.7175$	$90.20$	$2.4675$
efld0_direct2vec	$6.9914$	$82.60$	$2.5018$
imagevec	$4.9416$	$100.00$	$5.1024$
torquevec2	$0.6186$	$85.70$	$3.4521$
DHFR (CHARMM, no polarization)
Computational hotspots			$1.9535$
Total CPU time : 13.0355∗ s (100 steps)
elj1vec∗	$8.2493$	$64.60$	$1.8578$
echarge1vec∗	$3.6786$	$90.90$	$1.8297$
image (1)	$3.4130$	$0.00$	$1.0000$
imagevec∗	$1.1076$	$100.00$	$2.4751$

Table 4. Table 4: 1 step measured times of execution and boost factors for different test MS using Rel or Vec . Simulations ran on 1 core. Values are averaged over 10 steps.

MS	Ubiquitin
	Time Rel (s)	Time Vec (s)	Boost
time_vdw	$0.0964$	$0.0567$	$1.7002$
time_elec	$0.1352$	$0.0967$	$1.3981$
time_polar	$1.2326$	$0.8758$	$1.4062$
MS	DHFR
	Time Rel (s)	Time Vec (s)	Boost
time_vdw	$0.2359$	$0.1453$	$1.6235$
time_elec	$0.2823$	$0.2012$	$1.4031$
time_polar	$2.6051$	$1.8181$	$1.4329$
MS	COX-2
	Time Rel (s)	Time Vec (s)	Boost
time_vdw	$1.8906$	$1.1362$	$1.6639$
time_elec	$2.3816$	$1.7398$	$1.3689$
time_polar	$22.2782$	$15.6152$	$1.4267$
MS	STMV
	Time Rel (s)	Time Vec (s)	Boost
time_vdw	$1.9878$	$1.2260$	$1.6214$
time_elec	$3.8826$	$2.9314$	$1.3245$
time_polar	$64.2167$	$45.4406$	$1.4132$

Table 5. Table 5: Best production performances and boost factors for the different MS using Rel or Vec . For DHFR, COX-2, STMV and Ribosome, optimal results with CPU2 setup are also shown (see table  1 ).

AMOEBA polarizable Force Field
MS	$P_{\mathrm{𝐑𝐞𝐥}}$	$P_{\mathrm{𝐕𝐞𝐜}}$	Boost
	(ns/day)	(ns/day)	factors
Ubiquitin	$11.6875$	$16.9779$	$1.4526$
DHFR (CPU)	$9.1725$	$13.3312$	$1.4533$
DHFR (CPU2)	$9.4761$	$14.6054$	$1.5413$
Puddle	$3.5421$	$5.2417$	$1.4798$
COX-2 (CPU)	$1.9608$	$2.9343$	$1.4965$
COX-2 (CPU2)	$2.0137$	$3.1514$	$1.5650$
Pond	$1.7921$	$2.7620$	$1.5412$
Lake	$0.7066$	$1.1025$	$1.5602$
STMV (CPU)	$0.4995$	$0.7921$	$1.5858$
STMV (CPU2)	$0.5152$	$0.8140$	$1.5799$
Ribosome (CPU)	$0.2295$	$0.3420$	$1.4901$
Ribosome (CPU2)	$0.2368$	$0.3527$	$1.4894$

Table 6. Table 6: Best production performances and boost factors for different MS using Rel or Vec with CHARMM force field. For DHFR, optimal results with CPU2 setup are also shown (see table  1 ).

CHARMM non polarizable Force Field
MS	$P_{\mathrm{𝐑𝐞𝐥}}$	$P_{\mathrm{𝐕𝐞𝐜}}$	Boost
	(ns/day)	(ns/day)	factors
Ubiquitin	$39.3068$	$48.8269$	$1.2422$
DHFR (CPU)	$24.2333$	$31.7408$	$1.3098$
DHFR (CPU2)	$26.4805$	$34.8272$	$1.3152$
Puddle	$9.4749$	$12.8026$	$1.3512$
COX-2	$8.1411$	$11.3459$	$1.3936$
Pond	$5.1579$	$6.8394$	$1.3260$

Table 7. Table 7: Best production performances for the different MS using Rel2 , Rel2-multi (multi-timestep) and Vec2-multi (multi-timestep). For DHFR, COX-2, STMV and Ribosome, optimal results with CPU2 setup are also shown (see table  1 ).

AMOEBA polarizable Force Field
MS	$P_{\mathrm{𝐑𝐞𝐥𝟐}}$	$P_{\mathrm{𝐑𝐞𝐥𝟐}-\mathrm{𝐦𝐮𝐥𝐭𝐢}}$	$P_{\mathrm{𝐕𝐞𝐜𝟐}-\mathrm{𝐦𝐮𝐥𝐭𝐢}}$
	(ns/day)	(ns/day)	(ns/day)
Ubiquitin	$11.6875$	$28.28$	$40.32$
DHFR (CPU)	$9.1725$	$22.20$	$32.26$
DHFR (CPU2)	$9.4761$	$22.93$	$35.33$
Puddle	$3.5421$	$8.57$	$12.68$
COX-2 (CPU)	$1.9608$	$4.74$	$7.09$
COX-2 (CPU2)	$2.0137$	$4.87$	$7.65$
Pond	$1.7921$	$4.34$	$6.69$
Lake	$0.7066$	$1.70$	$2.65$
STMV (CPU)	$0.4995$	$1.21$	$1.92$
STMV (CPU2)	$0.5152$	$1.25$	$1.97$
Ribosome (CPU)	$0.2295$	$0.55$	$0.82$
Ribosome (CPU2)	$0.2368$	$0.57$	$0.85$

Table 8. Table 8: Profiling of Vec using Intel VTune Amplifier. Simulations ran on one core and 100 steps. MS is DHFR with AMOEBA polarizable force field and with CHARMM force field (no polarization). For the vectorized routines, the Vector Usage percentages go from 78.0 to 100%. As imagevec has been fully vectorized, there is no more separation for the image CPU time in the CHARMM part of the table. Neighbor list building routines have been fully vectorized. So, the lines with † † {\dagger} in the table  2 have been reintegrated in the total CPU time to compute the general boost factor.

Module	CPU Time	Vector	Boost
	(s)	Usage %	factor
DHFR (AMOEBA, polarizable)
Computational hotspots			$2.4566$
Total CPU time : 119.8797 s (100 steps)
tmatxb_pme2vec	$57.4947$	$100.00$	$1.7553$
epolar1vec	$21.0042$	$100.00$	$2.5094$
ehal1vec	$16.6440$	$78.10$	$3.1524$
empole1vec	$7.8784$	$90.20$	$3.6699$
efld0_direct2vec	$6.4751$	$100.00$	$2.7013$
vlistcellvec	$4.7310$	$100.00$	$1.5403$
mlistcellvec	$3.0972$	$100.00$	$2.6651$
imagevec	$2.5108$	$100.00$	$10.0422$
torquevec2	$0.4432$	$100.00$	$4.8183$
DHFR (CHARMM, no polarization)
Computational hotspots			$2.9201$
Total CPU time : 10.6051 s (100 steps)
elj1vec	$6.3047$	$75.00$	$2.4309$
echarge1vec	$2.6411$	$100.00$	$2.5486$
clistcellvec	$1.0424$	$100.00$	$2.7357$
imagevec	$0.6169$	$100.00$	$9.9763$

Equations

$$B=\frac{T_{\mathbf{Rel}}}{T_{\mathbf{Vec}}}$$

$$B=\frac{P_{\mathbf{Vec}}}{P_{\mathbf{Rel}}}$$

Full text

1 Introduction

Tinker-HP is a massively MPI parallel package dedicated to classical molecular dynamics (MD) and to multiscale simulations, especially using advanced polarizable force fields (PFF) encompassing distributed multipoles electrostatics[1]. It is an evolution of the popular Tinker package code [2] which conserves Tinker’s simplicity of use and its developer-friendliness, allowing for the rapid development of new algorithms. Tinker-HP offers the possibility to perform large scale simulations while retaining the Tinker reference double precision implementation dedicated to CPU-based (Central Processing Unit) petascale architectures. The parallel scalability of the software is demonstrated via benchmark calculations. Overall, a several thousand-fold acceleration over a single-core computation is observed for the largest molecular systems, thus allowing long reference polarizable MD simulations on large molecular systems containing up to millions of atoms.

Despite this strong acceleration, and due to the development model of the Tinker software suite (now version 8 [2]), which favours new scientific development over optimization, no attempt has previoustly been made to adapt the Tinker-HP code to a particular CPU architecture. Each execution of the existing code took little or no advantage of the underlying capabilities of the CPU it was running on. The deployment of numerous processors with SIMD (Single Instruction/Multiple Data) capabilities that in principle allow substantial speedups leads us to address this issue. The existence of strong vector capabilities on modern Intel®architectures, and particularly the Intel®Advanced Vector Extensions 512 (Intel AVX-512) on Intel® Xeon® Scalable and Intel Xeon® Phi™ processors, combined with the excellent support from the Intel® Fortran Compiler, motivates us to change the overall design of the code, while trying to keep its simplicity and readability. The goal of this paper is two-fold. First, it provides a comprehensive living review dedicated to the capabilities and performance of Tinker-HP’s main production methods (i.e. force fields) on Intel’s architectures. Second, it is organized to present in a pedagogical fashion our code optimization efforts, outlining an optimal strategy for vectorization of Tinker-HP. Although there have been theoretical discussions of AVX-512 vectorization[3], and a description of the effective vectorization of a Quantum Monte-Carlo production code[4], such practical examples are rarely documented and we think our experience could prove useful to other software developers.

The present version of the paper is organized as follows. After reviewing the specificities of the latest Intel Xeon Scalable Processors (code-named Skylake) and particularly their Intel AVX-512 vector instruction set, we present the general structure of the most CPU intensive Fortran subroutines in Tinker-HP, show their performance hotspots and propose a general strategy to vectorize the code within the Intel AVX-512 instruction set. We then give concrete examples of the vectorization process. Performance comparisons between the current Tinker-HP release and the vectorized version is then made, first for execution on a single core, to show brute force acceleration of the code, and then in the context of a realistic parallel execution (with up to 16 000 cores). Finally, extended benchmarks on meaningful systems are provided with the AMOEBA polarizable force field [5, 6, 7] and also using an initial implementation of classical force fields such as CHARMM [8], to illustrate what can be achieved with a typical non-polarizable force field (AMBER,[9] OPLS-AA[10] etc…).

2 Intel Xeon Scalable processors

We are using in our study a system with a CPU from the Intel Xeon Scalable processor family (code-named Skylake). These processors feature up to 28 cores per processor with two hyperthreads per core for a total of up to 56 threads per processor. A new mesh interconnect reduces the latency of communication within cores and controllers in the processor. Each core is capable of issuing up to four instructions per cycle, out-of-order. Up to two of these can be Intel AVX-512 instructions[11]. The Intel AVX-512 instruction set is a 512-bit SIMD instruction set that allows performing computations with a single instruction using SIMD registers that contain eight double-precision (DP) or sixteen single-precision (SP) floating-point values, as well as a variety of integer data sizes. The Intel AVX-512 instruction set also supports unaligned loads, fused-multiply and add, vector masking, shuffle and permutation instructions, histogram support, and hardware-accelerated transcendentals. Using all available cores and SIMD units is key to unlocking the full potential performance of these processors.

A significant change from its predecessor, the Intel Xeon processor v4, is the reorganization of the cache hierarchy to better balance the cache usage for server workloads. To achieve this, the L2 cache has increased in size to 1MB. The last level cache (LLC) has been reduced in size (up to 38.5 MBs), but it is now a non-inclusive cache (meaning that data is not evicted from caches closer to the cores when evicted from the LLC).

The Intel Xeon Scalable processors provide two memory controllers with three memory channels each that support DDR4 memory at up to 2600 MHz. This provides up to 123 GB/s of bandwidth to main memory for each socket. Three Intel®Ultra Path Interconnect links, each providing 10.4 GT/s, allow multiple processors to communicate to create bigger systems (e.g, dual-socket systems).

3 Considerations on vectorization

Efficient use of the SIMD units available in the processor is very important to achieve the best performance on modern (and probably future) processors. Most vector parallelism is extracted from loops. In this section, we outline some ideas that can help to achieve good performance from loop vectorization.

Modern compilers are able to auto-vectorize loops but need to be able to determine that the vectorization does not break any possible cross-iteration dependencies. This is not always possible due, for example, to variable aliasing or loop complexity[12]. Programmers can assist the compilers by rewriting their loops with idioms that compilers recognize or by using directives that guide the vectorization process.

Once a loop is vectorized, it is useful to consider if vector code generated for the loop is the most efficient possible. There are tools that can help with this assessment as well as provide feedback on how to improve efficiency[13]. Common issues that are worth looking into are:

• •

Unaligned loads or stores. When data is not aligned to cache boundaries, some penalty will be incurred in load and store instructions compared to well aligned data. Therefore, it is recommended to align data to cache boundaries and also use the proper mechanisms to inform the compiler.

• •

Loop prologues and remainders. To provide better data alignment, the compiler might generate different code for the first iterations of the loop to ensure the main vectorized portion of the loop runs on aligned data. Similarly, if the compiler cannot deduce the number of loop iterations, it will generate a remainder loop to compute the final iterations of a loop that do not fill a vector register completely.

• •

Unnecessary masking. While the Intel AVX-512 instruction set supports vectorization of loops with conditional statements, this feature requires the use of mask instructions and masked vector operations which reduces the vectorization efficiency. Therefore it is recommended to move, as much as possible, conditional statements out of vectorized loops. This might require spliting the loop (with different code for those iterations where the branch was taken vs. was not taken) or duplicating the loop (with different code for each branch of the conditional).

• •

Non-unit strides. Use of non-unit strides will often force the compiler to generate gather and/or scatter instructions which are less efficient than regular loads and stores. This also includes accessing a field in a structure, as a non-unit stride will be needed to access the same field between consecutive elements of an array. This is why a Struct-of-Arrays (SoA) layout is preferred over a more conventional Array-of-Structs (AoS) layout for efficient vectorization. Additionally, the Array-of-Structs-of-Arrays (AoSoA) layout can be used where a given number of elements of a field are placed consecutively in memory like in SoA, usually a small multiple of the cache line size, but element fields are still interleaved in memory as in AoS. This layout can improve locality when multiple fields are used closely together in a code section at the expense of increasing code complexity.

• •

Indirect accesses. Indexing one array with the values from another array will also result in generation of gather and/or scatter instructions. Therefore, this pattern should be avoided as much as possible.

• •

Register pressure. The number of vector registers is limited (e.g., Intel AVX-512 provides 32 vector registers). If the number of arrays that are used in a given loop, plus temporary values that might be required for the operation of the loop, exceeds the number of registers, then the compiler will need to spill some of the arrays to the stack and restore them afterwards which significantly reduces vectorization efficiency. To avoid this situation, it is better to use different loops for independent array operations rather a single big loop containing all arrays.

4 Working environment and definitions

4.1 Definitions

In this paper, we will use two versions of Tinker-HP :

the Release Version 1.1, referred to as Rel. 2. 2.

the Vectorized Version 1.1v, referred to as Vec.

A third version, Release Version 1.2 (referred to as Rel2), is mentioned in the Perspective subsection 8.3 to give a feeling for the performance gains anticipated with newly-developed algorithms. Vectorization of Rel2 is in progress. It will be referred to as Vec2.

We ran Tinker-HP exclusively on supercomputers under Unix/Linux Operating System (OS). These machines aggregate hundreds or even thousands of interconnected systems called computing nodes, or simply nodes, each typically having tens of CPU cores and hundreds of gigabytes of memory. On Unix/Linux computers, the code is executed by a process, which uses memory and CPU resources managed by the OS.

What we called the code can be split in two parts :

the User Code (UC), which comprises all the Fortran code. Here, it is the code for Tinker-HP and the 2Decomp library. 2. 2.

the Non–User Code (NUC), which comprises all the code executed by the process because of library calls from the UC, system calls done on behalf of the UC or code introduced implicitly by the compiler.

Of course, the way we write the UC (what library we used, how we setup our data,…) has an influence on the NUC executed by the process. As we want to raise the performance, we have to take into account what NUC gets executed because of the UC code that we write.

We use the term Molecular System (MS) to denote all the physical systems we have simulated for this paper.

Note that Fortran code listings shown in this paper have been taken as is, while all compilation reports and assembly code listings have been edited for publication purposes.

4.2 Compilation setup

We worked with the Intel® Parallel Studio XE 2018 development suite[14], containing the Intel Fortran Compiler, the Intel® MPI Library, the Intel® Math Kernel Library (Intel MKL) with implementations of Blas[15], Lapack[16] and Fftw3[17] routines, and the Intel®VTune™Amplifier for profiling and analysis.

The Tinker-HP sources are compiled with the flags shown in listing 1 where :

• •

-xCORE-AVX512 flag forces the generation of binaries for the Intel Xeon Scalable processors.

• •

-qopt-zmm-usage=high flag instructs the compiler to use zmm (512 bits) registers as much as possible.

• •

-align array64byte instructs the compiler to align all static arrays to 64 bits memory address boundaries

• •

-falign-functions=64 tells the compiler to align functions on 64 bits boundaries

• •

-qopt-report-phase and -qopt-report=5 flags produce vectorization reports.

• •

-S flag produces assembly code listings.

The vectorization reports are of major interest, because they give precise indications on how the compiler understands the code and what decisions it makes, allowing the programmer to modify the code or give the compiler pertinent directives. Some indications are of particular interest:

• •

the vector length, which gives the number of elements (real, integer or logical) affected by one operation

• •

the scalar cost, which gives the number of scalar operations necessary to handle one iteration of the loop

• •

the vector cost, which gives the same information as the scalar cost, but for vector operations

• •

the estimated potential speedup. Most of the time, this is the ratio scalar cost/vector cost. This speedup is not in time of execution, but in number of operations. A speedup of 10.0 does not mean that a loop will execute 10 times faster, but rather that there will be 10 times fewer vector operations.

Even if we do not use OpenMP*[18] in the dynamic simulation engine, other parts of the Tinkertools suite (Tinker 8 [2]) use OpenMP directives. So, object files are linked with the flags shown in listing 2 where :

• •

-mkl=sequential flag tells the linker to use the sequential Intel MKL library, which is lighter and faster than the multi-threaded ones (e.g., OpenMP).

• •

-qopenmp-stubs flag enables compilation of OpenMP programs in sequential mode.

4.3 Execution setup

For the performance tests, calculations were performed on nodes running under the RedHat* Enterprise Linux* Server Release 7.4 operating system, with the following configuration :

• •

2 $\times$ Intel Scalable Xeon 8168 processor – 2,7 GHz – 24 cores/processor

• •

192 GB of DDR4 memory,

• •

InfiniBand* EDR interconnect.

We chose 8 MS from those studied in [1] with increasing sizes ranging from 9 737 to 3 484 755 atoms : the Ubiquitin protein, the prototypic Dihydrofolate Reductase (DHFR), the COX-2 dimer, the Satellite Tobacco Mosaic Virus (STMV), the Ribosome full structure in polarizable water, and three water boxes (Puddle, Pond and Lake).

The table 1 gives for each MS the name, the number of atoms and the number of cores used for the parallel calculations.

All calculations have been made starting from an equilibrated MS with a timestep of $2fs$ and a RESPA integrator[19]. For the one core tests, we define the time $T$ (in seconds) as the time of execution. For the parallel tests, we define the performance $P$ (in ns/day) as the duration of simulation that can be achieved in one day. Our goal is to optimize throughput via vectorization. That means lowering the time $T$ and increasing the performance $P$. Thus, we define the boost factors $B$ as :

[TABLE]

where $T_{\mathbf{Rel}}$ and $T_{\mathbf{Vec}}$ are times for Rel and Vec respectively, or

[TABLE]

where $P_{\mathbf{Rel}}$ and $P_{\mathbf{Vec}}$ are performances for Rel and Vec respectively.

To get the profiles for Rel and Vec, we used DHFR MS with the AMOEBA polarizable force field and with the CHARMM classical force field (no polarization). Both simulations ran on one core and for 100 steps.

5 Release version 1.1 of Tinker-HP

5.1 Polarizable AMOEBA force field

We focus here on the part of the code dedicated to the AMOEBA polarizable force field which is the most computationally challenging and gives a lower bound for Tinker-HP performance[5, 20]. AMOEBA has been shown to have wide applicability for physical systems ranging from liquids to metals ions,[21, 22] including heavy metals,[23, 24] in gas and solution phase, to proteins [6, 7] and to DNA/RNA[7]. It uses distributed atomic multipoles through quadrupole moments and a Thole damped point dipole polarization interaction model. Van der Waals interactions use the Halgren buffered 14–7 function. In this paper, we used the AMOEBA 2013 protein parametrization [6, 7] coupled with the 2003 water model[5].

5.2 General Structure

Tinker-HP uses a 3D spatial decomposition to distribute atoms on compute cores. Every process is assigned to a subsection of the simulation box and is responsible of updating the positions, velocities and accelerations of the atoms present in this subdomain at each timestep[1]. The most computationally intensive part of Tinker-HP is devoted to forces and electric fields calculations.

All the compute routines follow the same organizational scheme :

• •

an external loop over all the atoms held by the process

• •

the selection of the neighbour sites using various criteria (e.g. cutoff distances, …)

• •

a long internal loop over the selected sites, where all quantities are computed.

In Rel, this internal loop computes all quantities for each atom-neighbour pair on the fly, with no attempt to pre-calculate or store intermediate quantities that can eventually be re-used. This results in a serious issue regarding cache registers and processing units, and in an extensive use of memory-core transfer instructions. By contrast, there’s almost no use of arrays, besides indexing. This means data are often not contiguous in memory, and therefore accesses to memory are irregular. Thus, the possibility to take advantage of the vector extension capabilities of the Intel AVX-512 instructions is very low.

5.3 Hotspots

The table 2 shows the profiling of Rel, running on one core for DHFR with AMOEBA (polarizable model) and with CHARMM (non-polarizable model). We give the name of the module or routine, the real CPU time spent executing it, and the vector usage percentage. All routines are sorted with higher time-consumption listed first. We can see that Rel has two kinds of hotspots :

NUC hotspots : these are mainly due to libraries calls, system calls and memory management operations (initialization, copy, allocation and deallocation). 2. 2.

Computational hotspots : these are mainly due to the computation of :

• •

the matrix-vector product operation applied at each iteration of the polarization solver (tmatxb_pme2), which can be called up to 12 times at each step, depending upon the convergence method

• •

the dipole polarization energy and forces (epolar1)

• •

the van der Waals energy and forces for the Halgren buffered 14–7 function (ehal1)

• •

the multipolar permanent electrostatic energy and forces (empole1)

• •

the right hand size of the polarization equation (efld0_direct2)

• •

the van der Waals energy and associated forces for the Lennard-Jones 6-12 function (elj1)

• •

the charge-charge interaction energy and associated forces (echarge1)

The routines used to build neighbor lists (vlistcell, mlistcell and clistcell ) appear. Other widely used utility routines (image and torque) also appear.

In order to raise the performances of Rel, we need to address these hotspots. Two observations guide us :

vmlinux is taking almost as much CPU time as the multipole polarization energy and forces computation routines. This means the process makes many system calls and performs many memory operations. 2. 2.

the vector usage percentage is strictly $0.00$ for all the computation subroutines. This confirms that these routines only use scalar operations.

The first observation led us to investigate the library and system calls and, first and foremost, to work on the memory management of a process running Tinker-HP in order to reduce memory operations.

The second observation forced us to rewrite the computation routines. Since using vector operations generally means using loops on arrays, the on the fly method of computation in Rel must no longer be used.

6 Optimization strategy

6.1 Non-User Code hotspots

NUC hotspots come from libraries calls, system calls (file open, read or close, MPI function calls,…) and memset, memcpy and calls to malloc() or its derivatives that each process makes during its life.

6.1.1 Libraries and System calls

Libraries calls

The vast majority of libraries calls comes from the Intel MKL library, which actually does computing work. As using it wisely can give a significant speedup, we have to provide the right Fortran code (UC) to access the right vectorized functions of Intel MKL.

System calls

When running, Tinker-HP reads a few files at the very beginning and outputs a log file which contains the simulation results and, periodically, a file containing atoms coordinates. As these input/output operations are done by the MPI–rank–0 process only, the open, read and close system calls do not really stress the Unix/Linux system, even if the simulation runs on millions of atoms.

So, most of the system calls come from the MPI library and are due to two design choices in Tinker-HP :

As explained before, each Tinker-HP process hold a portion of the space (a domain) and maintains MPI communications with MPI processes that hold other domains nearby it, so that each process can exchange information with the others and track atoms coming in or out of its own domain. As the other processes can run on other nodes, there can be even more time spent in system calls because of network transmissions. 2. 2.

A memory region is shared between all MPI processes on a computing node, using the sharing capabilities of the MPI library. The access control of this region is implemented through system calls (semaphores, flock()…) to synchronize processes and guarantee non overlapping when writing data.

Minimizing the time spent in system calls is not so easy. We tried different distributions of the MPI processes over the nodes to favour local MPI communications, but that did not give convincing results. We also tried to improve the use of MPI framework by masking communications with computations. We used non blocking versions of MPI_SEND and MPI_RECEIVE functions, and did some calculations before calling the corresponding MPI_WAIT. The performance gain is not really noticeable for now. But improving this part of the code is a work in progress.

6.1.2 Memory management

The figure 1 gives a simple picture of the memory layout of a Unix/Linux process. The text zone contains the code, the data zone contains initialized data. The bss zone contains uninitialized fixed size data and has itself a fixed and limited size, set at compile time. The heap zone contains allocated data and can grow or shrink, subject to allocations or de-allocations. The stack zone contains a Last-In-First-Out structure, where values are pushed and pulled at each subroutine or function call or if the processor runs out of free registers.

Process memory setup

Historically, Tinker dynamically allocates and de-allocates all arrays it uses, because it was originally built to run on workstations with limited amount of memory and cores. These allocations are made by calls to the system malloc() group of functions. As a consequence, data are put in the heap section of the process, whose size is managed by the OS, allowing it to grow or shrink as needed.

As Tinker-HP is a pure MPI program which distributes data and can potentially run on hundreds of different nodes, each of them with gigabytes of memory, the problem of the memory consumption is not that important. On a computing node, each core (so, each MPI process) can easily have 2 or even 4 gigabytes of memory for its own use.

Still the size of some arrays are proportional to the size of the systems and therefore, the MS-size dependent data, declared when entering a subroutine, can be very large. In a normal run, each process maintains hundreds of numbers for each atom it holds. And we can have thousands of atoms held by each process and tens of MPI processes on one node. So, allocation and de-allocation of data for each process and in each subroutine constitutes a big stress for the OS.

Considering that the overall size of data held by one process is often under 2Gb, and that we maintain size derived constants throughout the program, we decided to remove all dynamic allocations in the vectorized routines, and declare all arrays with fixed sizes. That moves the data in the bss section of the process, which is lighter than heap for the OS to handle, lowering the stress on it.

Memset and memcpy

The execution cost of the memset and memcpy operations cannot be easily evaluated, as they come from the compiler libraries and are built upon the C library. But, because of their potential (big) effect on the performances (see table 2 and discussion on page 2), these operations have been extensively tracked.

Many of memset come from unnecessary zeroing and have been easily removed. Some of them come from the use of intrinsic Fortran90 functions, where the compiler creates temporary storage and introduces memcpy operations to work with it (for example, PACK). We tried to remove the use of intrinsic functions as much as possible. Some of the memset or memcpy operations also come from the way Fortran passes arrays to subroutines or functions (as a whole array, or as a slice). We avoided these operations wherever possible.

After this optimization, the real CPU time on one core for NUC hotspots can be shorter by up to 10%. But this depends a lot on the MS simulated and the activity of the Unix/Linux system outside of Tinker-HP.

6.2 Computational hotspots

The strategy we used can be developed in five guidelines :

Rewrite all big internal loops. As using vector operations means using arrays, big loops can be split in numerous short ones, each loop computing only one or a few quantities for all the involved atoms. This way, the quantities calculated can be kept in arrays and vector operations can be executed on them. 2. 2.

Cope with the way the compiler works on loops. As stated in section 3, when the compiler tries to vectorize a loop, it can generate 3 execution loops :

• •

a Peeled loop (P–loop), to treat array elements up to the first aligned one.

• •

a Kernel loop (K–loop), to treat the biggest possible number of array elements with vector operations.

• •

a Remainder loop (R–loop), to treat elements that remain untreated by the previous loops.

As the K–loops are the most effective and the fastest loops, we must eliminate as much P–loops and R–loops as possible. We’ll show below what we did to achieve this goal. 3. 3.

Use vectorized mathematical operations as much as possible. This can be difficult sometimes, because each compiler or library implements them in its own way. For example, using the Intel Compiler, the sqrt(X) function is not vectorized. But the power function ** is. So loops with X**0.5 have a better vectorization score than loops with sqrt(X). As the two functions can give slightly different numerical results, care must be taken to be sure to always have the same numerical results than the Rel. 4. 4.

Limit the use of Fortran intrinsics The use of Fortran intrinsics (SUM, PACK, etc…) is a common approach for vectorization in Molecular Dynamics package. But we found that, in many cases, the presence of intrinsics can prevent the compiler from finding the best optimization (see discussion about PACK on page 7). Intrinsics tend to protect data as much as possible. In doing so, they make frequent use of memset or memcpy operations which tends to result in unnecessary memory operations. 5. 5.

Have no dependency between arrays in the loops, because the compiler will refuse to vectorize any loop where it cannot be sure that there is no dependency.

With that in mind, knowing that Intel AVX-512 vector registers can hold eight 8-bytes reals or sixteen 4-bytes integers, we should have a significant improvement of the speed if the number of neighbouring atoms is big enough to fill them. That is generally the case in Tinker-HP calculations, except for very in-homogeneous MS.

To summarize, filling in the 512 bits registers in an efficient way and using as much vector operations as possible in a loop need :

• •

No dependency, to be vectorized

• •

Low number of arrays used, to reduce the register pressure

• •

Arrays as close as possible in memory, to reduce cache miss and cost of memory operations

• •

Data aligned on a suitable boundary, to eliminate the P–loop

• •

No subroutine calls, no un-vectorized math operations, to get the best of the K–loop.

• •

Loop count carefully chosen, to eliminate the R–loop

• •

No if-test. If tests are mandatory (as in the selection process), they should be built in logical arrays before being used in the loop.

6.2.1 Dependency

Short loops calculate only one or a few unrelated quantities. They use the lower possible number of arrays. Thus, dependencies do not often occur. Where the non dependency cannot be automatically determined, we can easily see it and give directives to the compiler, or at worst rewrite the loop.

6.2.2 Data alignment

Historically, in Tinker, data were located in commons, that were themselves organized with scientific development in mind. Some compilers have options to align commons. But they may be inefficient if data are not correctly organized, with memory representation in mind.

We decided to replace commons with modules, which have many advantages :

• •

Arrays can be efficiently aligned using directives when declared

• •

Overall readability is better, due to modularity

• •

Code can be introduced in modules, so we can group operations and write them once and for all.

In all the modules, we used an ATTRIBUTE ALIGN::64 directive for each array declaration. At the very beginning of this work, we used arrays like pos(n,3) to represent, for example, the three spatial coordinates of an atom. But, sometimes, we saw that the initial alignment of the first row of pos was not kept by the compiler for the following ones, preventing it from fully optimizing the code and forcing it to generate extra P–loops. All arrays are now mono-dimensional. The coordinates are represented by arrays like Xpos(n), Ypos(n) and Zpos(n). The three coordinates are treated in the same loop, with no dependency, allowing for vectorization.

6.2.3 Data layouts in memory

The figure 2 shows 3 different data layouts for arrays in memory.

• •

In the setup ①, no ATTRIBUTE ALIGN::64 directive has been given. There is no memory loss, but the real8 array is not on a 64bits boundary. During execution, elements in this array will not be aligned. If no ASSUME_ALIGNED::64 directive is given, the compiler will generate P-loops. If an ASSUME_ALIGNED::64 directive is given, no P-loop will be generated. But the process will pick up wrong real8 numbers in the K–loop, and give wrong results, or even crash.

• •

In the setup ②, all arrays are aligned on a 64 bits boundary with an ATTRIBUTE ALIGN::64 directive. No P-loop will be generated. But there can be a memory hole, if the number of elements in the integer arrays is odd. When running, the process could have to do some jumps in memory to get the real*8 numbers, loosing time.

• •

In the setup ③, all arrays are aligned on 64bits with an ATTRIBUTE ALIGN::64 directive. No P-loop will be generated. There’s no memory hole, because the number of elements in the integer arrays is kept even.

We decided to implement the setup ③, which represents the best trade-off between speed and memory consumption. So, we ended up with the typical array declarations in a module (here, MOD_vec_vdw) shown in listing 3, where :

• •

All 4-bytes integer arrays are listed before 8-bytes real ones, and each parametric dimension is a multiple of 64, eliminating the holes in the data layout and ensuring correct alignment in any cases. For example, in the listing 3, the parameter maxvlst, which represents the maximum number of neighbours in van der Waals pair list, is set to 2560.

• •

Arrays used in a loop are listed close to each others in modules and in the same order of utilization as much as possible, to reduce the time for memory operations.

For arrays declared in the subroutines (outside of modules), with shapes depending on the size of the MS, we calculated the next multiple of 16 bigger than the size, and used it as the dimension of the array.

Although the setup ③ implies an over-consumption of memory for arrays with MS–dependent sizes, this is almost not noticeable, because we add at worst 15 elements to the size, which goes usually from thousands to millions in Tinker-HP.

Almost all P–loops have been removed this way.

6.2.4 Loop counts

As we have many short loops in each subroutine, we really must carefully choose the loop counts. The number of sites being dependent of the size of the MS we work on, we cannot impose a fixed value. To be sure to completely fill the 512 bits registers at each loop, we decided to maintain two working loop counts :

a Real loop count, multiple of 8 for loops on real*8 arrays ($8*8=64$ bytes, so 512 bits). 2. 2.

an Integer loop count, multiple of 16 for loops on integer arrays ($16*4=64$ bytes, so 512 bits)

As an example, if nnvlst is the number of sites we work on, these working loop counts are computed as in listing 4.

Since nnvlst can already be a multiple of 8 or 16, we should use the mod and merge constructs to get the smallest loop count. We used nvloop8 for real8 operations, and nvloop16 for integer operations. real8 arrays are loaded in registers by chunks of 8 elements, and integer arrays by chunks of 16 elements. We eliminated almost all the R–loops this way.

The flaw of this method is that, when nnvlst is very low, we do an overwork. However, the number of neighbours is generally big enough (typically between 200 and 2000), so these extra iterations are not too expensive. Furthermore, they are executed as K–loops, where maximum vectorization is in effect. We just have to remember that we have now extra (useless) calculated values and find a way to drop them in our calculations.

6.2.5 Design of loops

The use of short loops allows the programmer to better understand the code. That is exactly the same for the compiler ! So, loops in Vec are :

• •

simple. We use a small number of arrays (a maximum of 8 seems to be a good start). Due to the intrinsic complexity of the mathematical model, we use many temporary arrays.

• •

short. The loops contains 10 instructions at most, and 3 or 4 most of the time.

• •

mostly if-test free. Most of the time, an if-test in a loop prevents the compiler from vectorizing. For the loops in the neighbours selection process, which cannot be if-test free by nature, tests are built as logical array masks before being used. This way, loops using these masks can easily be vectorized.

6.2.6 Editing considerations

The large refactoring effort on the Rel code may seem to be of limited efficacy at first glance. But we found that this allowed us to better understand what the code does, which was a crucial step in replacing the internal big loop with multiple short loops. It was then far easier to reason about vectorization, and modifications and debugging were simplified.

7 Vectorized loops in Tinker-HP

We have 2 kinds of vectorized loops in Tinker-HP :

• •

selection loops that select sites using various cutoffs and work on integers or logicals

• •

compute loops that compute quantities and work on reals.

We will show each of them in details, and give some insights on how we have improved the vectorization. The typical loops have been extracted from ehal1vec, which use the module MOD_vec_vdw shown in listing 3.

7.1 Typical selection loop

We built a selection mask mask1 with the appropriate test, using nvloop16 here, as we work on integers (listing 5). We first tell the compiler to assume the loop count value is a multiple of 16, so that it does not generate any R–loop. We really have to be sure of that, otherwise the process will pick up numbers from other arrays in the memory and wrong results or bad crashes will occur.

The vectorization report on listing 6 shows that the loop is perfectly vectorized.

We then applied the mask on the set of atoms we worked on to select those of interest. In Fortran, we have an intrinsic function PACK that does exactly this. So, in a first attempt, we wrote the selection loop as shown in the listing 7.

Unfortunately, because of the PACK function, which does an implicit loop over all the elements of the array it works on, each line was seen as an independent loop by the compiler, and optimization was made on that line only.

The corresponding Fortran compiler report on listing 8 shows that a vectorized loop is generated for the count line, with a vector cost of $0.810$. Then, the first PACK line generates 3 loops :

• •

one over the 2-bytes logical array mask1 (vector length 32) with a vector cost of $1.250$.

• •

one over the 4-bytes integer array kglobvec (vector length 16) with a vector cost of $1.000$.

• •

one reduced to a memset or memcpy for the assignment of kglobvec1

The total vector cost is $2.250$ for one PACK operation. We also have 3 loads and 1 store for each PACK line.

For this selection loop, we obtained a total vector cost of $0.810+4*2.250=9.810$, plus 4 memset or memcpy, and a total of 13 loads and 4 stores. We cannot easily know the final cost of this selection loop, because, as stated before, the implementation of the 4 memory operations depends on the compiler.

To get a controlled and constant vector cost, whichever compiler we use, we decided to get rid of the PACK function.

After all, packing data is just a matter of selecting array elements and putting them contiguously in a new array. So, we ended up with a functionally equivalent loop depicted in listing 9.

Although there is a test in this loop, the corresponding Fortran compiler report (listing 10) clearly shows that :

• •

The loop is vectorized

• •

Every reference is aligned, so are the loads. The stores cannot be aligned, because of the packing.

• •

The vector length is 16 which means 16 integers will be picked up in each operation

• •

The potential speedup is more than 7. This is very good in the presence of a test.

• •

4 vector compress instructions are generated, which correspond to the 4 assignments.

The last remark is very interesting : the Intel compiler was able to recognize this construct as a packing loop, and implemented it directly with vpcompressd instructions, which are Intel AVX-512 pack instructions at the machine code level.

A look to the assembly code in listing 11 confirms that the vpcompressd instructions operate on zmm pure 512-bit vector registers.

We obtained a vector cost of only $2.500$ and no memset or memcpy. That is 4 times smaller than the PACK version, and much more if we count the memset or memcpy operations. The number of loads is also reduced by a factor of 3. This version is really faster than the first one.

7.2 Typical compute loop

The calculation loops follow the scheme we have described above. They are short, simple and easy to read. The listing 12 shows a typical compute loop.

We first tell the compiler to assume the loop count value is a multiple of 8 (we work on reals here). All arrays are independent and used in the order they were declared in the module (see listing 3).

We can see from the corresponding Fortran compiler report in the listing 13 that :

• •

The loop is vectorized (no dependency).

• •

Every reference is aligned, so are the loads and stores.

• •

The vector length is 8 which means 8 numbers will be picked up in each operation

• •

The potential speedup is around $10.5$.

• •

2 vectorized math library calls are made for the 2 ** function.

A look to the assembly code in listing 15 shows that all multiplications are done with vector operations vmulpd and vfmadd213pd and vfmadd231pd, which are fused multiply-add operations. These vector instructions operate on zmm registers. We can also see the two calls to the vectorized version of the ** function so we are fully using Intel AVX-512 capabilities.

If ever we had used a division, instead of **- one, we would have got :

The estimated potential speedup in this case is less than half the previous one. And the utilization of the vector units is not so optimal.

So, a careful reading of the vectorization report is always necessary to ensure the best choices have been made.

7.2.1 Final profile for Vec.

The tables 3 shows the profile and the boost factors between Rel and Vec for the final vectorized routines.

NUC hotspots

The Real CPU Time is almost the same as for the Rel version. We can see a reduction of about 10% for the real CPU time of vmlinux. The libmkl_vml_avx512.so shared library has appeared, because we use vectorized mathematical functions and replaced all calls to the complementary error function erfc, which was in the sources of Rel, by calls to vderfc, which is a vectorized implementation in Intel MKL library.

Computational hotspots

The vector usage percentage varies between 64% and 100%, and the boost factors are between 1.60 and more than 5.10. The real CPU time has shrunk from $278.95s$ to $136.45s$ for AMOEBA calculation, and from $24.80s$ to $13.03s$ for CHARMM calculation, giving an overall boost factor of roughly 2.

Naive users may expect a higher boost, since Intel AVX-512 registers are designed to hold 8 reals and vector operations are built to apply on the 8 reals at the same time. But other phenomena, like indirect indexing of arrays, memory operations, vectorization coverage, and processor frequency changes are also in effect which limits the boost that can be achieved in practice. Also, even if the reported speedups for each loop can be between 6 and 15 or even more, some masked load and store operations also tend to lower the boost factor.

In a recent study, Watanabe and Nakagawa[25] have obtained a boost factor between 1.4 and 1.8, depending on the data setup, for the vectorization of the Lennard-Jones potential on AVX2 and AVX-512 architectures. So, achieving comparable, and often superior, boost factors on all the vectorized routines of Tinker-HP seems quite satisfactory. One interesting way to see what could be achieved is also to look at what has been done in another community code developed at Argonne and devoted to Quantum Monte-Carlo[4]. By definition, Monte-Carlo is the opposite of MD, as it is a highly parallel approach with few communications. In that context, gains up to 4.5 have been obtained, highlighting some kind of upper limit of what is possible to achieve.

Watanabe and Nakagawa also showed that the boost factors are very dependent of the structure of data. This means a code written with AVX-512 capabilities in mind will not be so efficient on AVX2-capable processors. Thus, developers need to adapt their code to a specific processor in order to obtain significant boosts. But the general strategy used here should remain useful on all kind of architectures.

8 Performance on Intel Scalable Processors

8.1 Sequential performance

We have evaluated the overall performance boost due to Intel AVX-512 by running calculations on only one core from a dedicated node. In this situation, we can easily measure the execution time of each interesting part of the code, with limited perturbation from the MPI portions of the code or the presence of other processes.

We chose to measure three execution times :

time_vdw, which is the time taken by Van der Waals calculations. Depending on the setup, we used ehal1(vec) or elj1(vec). 2. 2.

time_elec, which is the time taken by electrostatic calculations. Depending on the setup, we used echarge1(vec) (direct charges) + reciprocal charges or empole1(vec) (direct multipoles) + reciprocal multipoles + torque(vec2). 3. 3.

time_polar, which is the time taken by polarization calculation and by far the biggest. We used epolar1(vec) (direct polarization) + efld0_direct(vec) + reciprocal polarization + 9 calls to tmatxb_pme2(vec) + 2 calls to torque(vec2)

Notice that, as the times we measured are for the execution of a combination of vectorized subroutines and non-vectorized subroutines, they cannot be directly related to the CPU times reported in table 3.

In this case, the boost is always a tradeoff between feeding the CPU with enough numbers, which goes better with the size of the MS, and minimizing the exchanges between memory and cores, which goes worse with the size.

For every MS, we ran 10 calculation steps using the polarizable AMOEBA force field. We took the average value of each time, removing the smallest and the biggest. Results are given in the table 4.

The performance boost factors are always very good, even when the size of the MS is quite large (more than 1 million atoms on one core for STMV !). The boosts we obtained are significant and justify the important vectorization efforts we made to get them.

The real performance gain should be estimated in a more realistic situation, where Tinker-HP is running in parallel. In this case, there could be 8 to 48 processes running on one node, each competing for resources, and up to 340 nodes involved, multiplying MPI communications.

8.2 Parallel performance

8.2.1 Polarizable force field : AMOEBA

We focus here on the absolute performance improvements over previous published results. As vectorization did not affect the scaling of the methods, interested readers can refer to the earlier Tinker-HP software publication for a detailed analysis of scalability [1]. Here, we ran calculations of 2 000 steps (4ps), with the core setups shown in table 1. The best performance was taken as the average of the 20 performance evaluations made by Tinker-HP during the run, after removing the first, middle and last values, which are lower because of the writing of intermediate files at these steps.

Results given in table 5 show a boost factor between 1.45 and 1.59 in parallel mode. The boost increases with the size of the MS, indicating a better overall utilization of the vector registers. When the MS is large, other phenomena (MPI memory contention, network communications, …) result in lower boost factors. We are still able to obtain small gains with CPU2 sets, because most of the supplementary cores use vectorized routines. The results are very encouraging, especially given that not all the code has been optimized.

We pushed forward and tried simulations on STMV and Ribosome with up to 16200 cores (CPU2 set). Figures 3 and 4 show the performance obtained for Rel and Vec upon increasing the number of cores.

The boost factors remain relatively constant for these two MS. With very large number of cores (and very large number of nodes), both Rel and Vec speeds are bounded by MPI communication and memory operations.

8.2.2 Non-polarizable force field : CHARMM

Tinker-HP is not yet optimized for traditional simple partial charge force fields as no specific "modern algorithmics" are present. Indeed, our current implementation is essentially a massively parallel version of the initial Tinker code that was aimed toward performance comparisons for Steered Molecular dynamics between polarizable and non-polarizable approaches.[26] Also, we have tested a conservative molecular dynamics setup where bonds lengths are not rigid, reciprocal space computations are done at each timestep, etc. Such a setup was chosen in order to provide reference numbers, but the future performances can likely be accelerated substantially. At present, we need much more cores to get results comparable to those of other prominent codes[27, 28, 29, 30]. Nevertheless, we decided to make performances measurements, firstly to get an idea of the boost that vectorization can provide in this case and, secondly, to know if we can still benefit from the scalability of the code, which is one of its greatest strengths. We used the same MS and the same CPU sets, limited to a maximum of 2 400 cores (i.e. as they were chosen for AMOEBA).

Vectorization boost

The table 6 shows the performances we obtained for Rel and Vec.

Overall, the speedup factor in using non-polarizable force fields is found to be between 3 and 4. The boost factors are lower than for AMOEBA, mainly because the vectorized part of the code which is actually executed is itself smaller. The results show the same behaviour as for AMOEBA as the size of the MS increases, with a peak value reached for smaller systems (around 200 000 atoms). Beyond this size, the non-vectorized code become the limiting speed factor.

Scalability

We tested the scalability of the code with three MS : Ubiquitin, DHFR and COX-2. As for AMOEBA, we ran for 2000 steps with increasing number of cores, and took the average performance given by the code. Figures 5, 6 and 7 show the performance obtained for Rel and Vec when increasing the number of cores.

For all MS simulated, the scalability is still very good. The boost factor remains almost constant for the two smaller MS (Ubiquitin and DHFR). For COX-2, the boost factor decreases from 1.41 to 1.39 when increasing the number of cores because, with 2 400 cores, communications tends to lower the benefit of the vectorization. In practice, this version of the code is a first step towards an efficient engine for non-polarizable MD but work is still required and is in progress to obtain better performance with updated code.

8.3 Perspectives on Tinker-HP 1.2 performance

This section gives an indication of the performance gains that will appear in the forthcoming Tinker-HP Release 1.2 version (Rel2). Indeed, despite being not fully vectorized, this major update proposes significant algorithmic speedups. For now, we can point out that a strong performance gain without accuracy loss is observed in using Tinker-HP with the new multi-timestep BAOAB-RESPA1 integrator[31], and with hydrogen mass repartitioning. This newly introduced integrator splits the energy terms in three levels evaluated at different timesteps: the bonded terms are evaluated every 1 fs, the non-bonded terms (including polarization) are split into short and long range, the short-range being evaluated every 10/3 fs and the long range every 10 fs. Furthermore, short-range polarization is treated with the non-iterative TCG-1 (Truncated Conjugate Gradient) solver[32, 33] and the outer-level uses the Divide-and-Conquer Jacobi Iterations (DC-JI) [34] approach, offering a net global acceleration a factor of 4.91 compared to standard 1 fs/Beeman/ASPC (7 without ASPC) simulations without loss of accuracy, enabling an accurate evaluation of properties such as free energies[31].

Preliminary results (where not all routines are yet vectorized) are reported in Table 7. We intend to review the full 1.2 vectorized version Tinker-HP in a future update of this living review.

As of July 2019, we have vectorized the neighbor list building routines and made improvements to the vectorization of other routines. The table 8 shows the new boost factors for the computational hotspots. The individual speedups for neighbor list subroutines range from $1.54$ to $2.73$, raising the overall vectorization boost factor to $2.45$ for the polarizable force field and $2.92$ for classical force fields. The overall performance boost (parallel gain) for classical forces field thus increases from 1.4 to 2.0 bringing interesting perspectives towards a future highly optimized classical force field MD engine.

Finally, beside the focus on the AMOEBA polarizable force field, performances will be given for other polarizable models as well as on classical force fields (CHARMM, AMBER, OPLS etc…). For now, despite the non-optimization and the absence of use of lower precision of this part of the code, more than a 4-time speedup of the values reported in Table 7 give an initial idea of the reasonable code performances for non-polarizable simulations.

9 Conclusion

In many ways this work represents a fundamental step in the evolution of the Tinker-HP software.

First, it demonstrates that new HPC architectures can offer significant acceleration to an existing massively parallel code like Tinker-HP. A brute performance boost factor between 1.32 and 1.70 can be achieved on computationally intensive parts of the code, leading to an overall acceleration factor between 1.45 and 1.59 for AMOEBA (1.24 and 1.40 for CHARMM) under realistic conditions, including the simulation of molecular systems with millions of atoms. Considering that the many current calculations require a total simulation time of a several microseconds, such a speed gains represent major progress.

Second, it shows that improved speed is not just available from raising the frequency of the CPU or buying more powerful computers. Large accelerations such as those reported here involve a close cooperation between the computational chemists, who write the code, and HPC specialists, who know how the CPU and the system software work. To get these gains, we had to dig into the pieces of code that were the most CPU consuming and to rewrite them almost completely, with simplicity and efficiency in mind. But it was worth the effort. Furthermore, considering the trends observed with prior CPUs, we anticipate that vectorization will also play an important role in future architectures.

Third, this work gives us a strategy and some methods to further improve the code. It can serve as a solid starting point for the future. We are now able to more easily adapt Tinker-HP to new underlying hardware or software advances. That will allow us to make the best of new technologies.

Of course, optimization is far from finished as some parts of the code are not yet vectorized (for example the reciprocal space computations involved in permanent electrostatics and polarization), and other sources of possible speedups exist and will be investigated. In particular, we have to review how we can improve the creation of neighbour lists, implement faster indexing of all the atoms (sorting indexes could be a solution) and achieve better masked MPI communications within computations. Decreasing precision is also possible in specific cases to gain performances while retaining sufficient accuracy. This paper will continue to be updated as we accumulate new data on Github, until a new version of this living document is pushed to review. The next iteration of the paper will also incorporate results on next-generation of Intel Xeon Scalable processors (codenamed Cascade Lake), and attempt to evolve towards an adaptation of the Tinker-HP code to any future architectures proposed by Intel. Future work will focus on the algorithmic boosting of our initial implementation of classical non-polarizable force fields. In addition, the next iteration of the paper will propose more detailed benchmarks for new polarizable approaches, including SIBFA [35, 36] and ongoing modifications of AMOEBA such as AMOEBA+[37] and HIPPO[38, 39].

10 Acknowledgments

We want to give special thanks to David Guibert and Cedric Bourasset (Center for Excellence in Parallel Programming, Atos-Bull) for initial discussions.

LHJ wishes to thank Prof Esmaïl Alikhani (Sorbonne Université) for sharing his 2x18 cores Skylake computer in the early phase of this work.

LHJ also wants to thank Bernard Voillequin for fruitful discussions at the very beginning of this work. B. Voillequin is not a scientist, and not even an engineer. But he’s a man of good sense and good will. That is invaluable.

11 Funding Information

This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 810367), project EMC2 (JPP). JWP and PR acknowledge support from the U.S. National Institutes of Health (R01GM106137 and R01GM114237). The authors acknowledge the generous support provided by GENCI (Grand Equipement National de Calcul Intensif, France) through the Grand Challenge and grant no A0050707671, which provided computing allocations on the Irène Joliot-Curie (TGCC, Bruyère le Châtel, France), Occigen and Frioul machines (CINES, Montpellier, France). The authors also thank the TACC (Texas Advanced Computer Center) for computational resources provided on the Stampede-2 machine.

Disclaimers

Software and workloads used in performance tests may have been optimized for performance only on Intel microprocessors.

Performance tests, such as SYSmark and MobileMark, are measured using specific computer systems, components, software, operations and functions. Any change to any of those factors may cause the results to vary. You should consult other information and performance tests to assist you in fully evaluating any contemplated purchases, including the performance of that product when combined with other products.

For more information go to www.intel.com/benchmarks.

Performance results are based on testing as of March 14th 2019 and may not reflect all publicly available security updates. See configurations described in the paper for details. No product or component can be absolutely secure.

Intel technologies’ features and benefits depend on system configuration and may require enabled hardware, software or service activation. Performance varies depending on system configuration. Check with your system manufacturer or retailer or learn more at intel.com.

Intel, the Intel logo, Xeon, VTune and Xeon Phi are trademarks of Intel Corporation or its subsidiaries in the U.S. and/or other countries.

*Other names and brands may be claimed as the property of others.

Author Information

\makeorcid

[TABLE]