Addressing the challenges of standalone multi-core simulations in molecular dynamics

1 Introduction

The last few decades have seen a rapid rise of computational modelling as a powerful tool to develop materials by postulation, inference and tuning in what can only be described as a semi-closed loop approach. The advantage of this approach is that properties that may not be directly ascertained, such as homogeneous nucleation [1], or are costly to investigate empirically can be evaluated [2, 3, 4, 5, 6]. A molecular simulation begins with choosing the best force field that describes the molecular configuration and ends with the iterative calculation of the phase space under specified boundary conditions that denote the external influences. Such conditions convey the applied forces, thermostats [7], pressure, etc. The choice of mathematical model, which in effect denotes the force field, affects the simulation data structures and therefore the design of the software. This in turn affects what can ultimately be parallelized. In our own work, we focus on the Sutton–Chen implementation of the Finnis–Sinclair embedded atom model potential [3, 7, 8, 9]. Force fields describe the interactions of particles in the configuration in terms of normalized dimensionless perturbations. The general assumption is that the internal force field in the molecular configuration is usually much stronger in comparison with the external perturbation. The required phase can be determined deterministically by molecular dynamics (MD), or stochastically through Monte-Carlo approaches [7, 10]. Ultimately, the results of such computation must be compared with the experiment to ascertain its accuracy or even to gauge the efficacy of the model. This leads either to a tuning of the model if at the model development stage or to a tuning of the material properties by parameter adjustment if the model has been tested and is accepted previously. There are generally two simulation approaches used, namely MD and ab-initio methods [11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. The basic computational modelling approach takes the aggregated, macroscopic properties arising as the consequence of their isolated, pairwise and interactive atomistic dynamics and considering the cumulative effects of these atomistic interactions as aggregations by some method that incrementally integrates their equations of motion. The recurrent idea in both approaches is to first derive a mathematical model that best describes the system, through which a state trajectory can be found. The mathematical model and its effectiveness and range depend on a number of factors that ultimately decide the extent of success of the simulation. Therefore, for a given system, it is imperative to derive a model that is as realistic as possible. However, this particular goal must take into consideration the level of accuracy required and the computational resources available. For atomic and molecular systems of macroscopic significance, the number of particles in the system can be extremely large, and studying their pairwise interactions will generate vast amounts of data and consume equally vast amounts of system time and resources. Using a faster processor does not necessarily alleviate the problem immediately, in spite of any mathematical simplifications such as specification of a cut-off distance for the interactions, simplification of the force field and so on that may have been devised. As interest in the computational study of materials grew, the issue of computational speed and accuracy quickly took centre stage and highlighted the shortcomings of the typical computer for these purposes. Compounding the situation is the fact that conformity to Moore’s law [21, 22] is no longer assured using the device interconnect-size concept on which the law is based. Moore’s law predicts a doubling of computational power approximately every 2 years as a consequence of the falling distances of the component interconnections on a semiconductor wafer. In 1971, the interconnect distance was around 10 μm and currently stands around 10 nm. It is expected to reach 5 nm around 2020. This distance affects the fastest speed that a charge conveying bit information can travel between any two devices on the same chip. Operating speed, however, is not simply a consequence of the interconnect distance. Other engineering considerations must be made, such as how heat dissipation will be handled. Also, as the physics shows, this distance is the realm of the lattice parameter (around 10 times the 0.361 nm lattice parameter of copper, for instance) where quantum mechanical effects manifest significantly. Consequently, processor speeds have more or less stabilized over the last 5 years at clocking frequencies of around 3 GHz. Clearly, to keep up with the demands of increasing computational power, there had to be a shift in the thinking. Today, vast increases in computational power are therefore not so much due to unilateral improvements in processor technology, but due to ingeniously interconnected multi-core, multiprocessor arrays. The development of computational frameworks and efficient algorithms appears to be driving innovations in atomistic MD simulations [7, 8, 9]. This is unwittingly a consequence of advancements in computer gaming, where the demand for high-performance computer graphics for more realistic rendering of ever larger and often networked games is ever greater and far outweighs other applications in the public domain. Thus, better models and computational algorithms are increasingly possible that exploit these features and are supported by ongoing developments in internet and computer technology (ICT). The rise of high-performance, supercomputing clusters [6] as low-cost supercomputers may perhaps be considered a secondary revolution in computation. This is because micron-scale, microprocessor interconnects [23] are now falling as fast as the demands of processing power require. However, in spite of these trends, standalone computation is still of considerable importance because it underscores the role of the individual node computer which is an integral component of any cluster, particularly if its inherent power could be unlocked through careful software design. The latter remains the basic computational element, and several such distributed nodes collectively execute chunks of tasks assigned by a task manager, which is also responsible for collating the results of the processing and making it available to the requester. Unfortunately, apart from the nature of interconnects between processors or cores themselves, the design of the node can easily become the weakest link in a computational arrangement and can lead to overheads whose cumulative effects substantially diminish the overall system performance. By standalone, one does not mean that the system is necessarily single processor or single core, but only that the computing node does not feature as part of a distributed system. Standalone can therefore mean any arrangement on a physically isolated, single machine that can support anywhere from single-core to multiple single-core processors to multiple processor multi-core computers. Such arrangements are today quite ubiquitous as standard features on desktop computers. For instance, the i3, i5 and i7 classes of x86-based machines common on desktops today are single processors with two to four cores and can support hyper-threading, such as Turbo Boost and K-technology [24]. The K-technology allows the overclocking of the CPU but for computational purposes illustrates the law of diminishing returns since it can limit the CPU power to only one core. Only applications written with one core in mind will likely reap any speedup benefit from such overclocking. The i7-Extreme can have up to eight cores and offer a unique processing opportunity if properly coded for the extra capability. Typically, most desktop applications do not use that extra feature and therefore do not observe the extra speedup. Additionally, on the same main board supporting the processor, there can be a high-definition hardware graphics processing unit (GPU). This multi-core CPU/GPU combination is the “typical” standalone configuration being referred to here. In our own work, the simulations are done on a Dell Optiplex 3010, which hosts a third-generation i5-3470 (3.2 GHz) with four cores and L3 cache (6 MB), has no hyper-threading but has K-technology (not activated to the maximum 3.6 GHz in this work). The system comes preinstalled with 8 GB (DDR3–1.6 GHz) of SDRAM. The GPU is AMD Radeon HD 7470 with 1 GB DDR3 SDRAM and interfaces the CPU with a PCI Express, 16-bit bus. Our results presented below are based on this hardware specification. The external storage, which is not part of the benchmark performance and mentioned here for completeness, is a 500 GB hard drive. Its function is merely to hold initial particle data and boundary conditions and to store the ultimate output of a simulation. Once the simulations start, no intermediate file access is done, all intermediate calculations being done in a SDRAM scratch pad. In the typical computer, the random access memory (RAM) is shared between cores or processors through a common, managed memory area called a “heap.” There have also been dramatic developments in physical storage devices associated with the modern desktop to address the demands of traditional multimedia data files, if not for the more demanding storage of large gaming applications and data. For computational applications, physical storage media may be used as swap space and intermediate storage, although in the interest of speed they are generally relegated to input and ultimate data output. These considerations put a heavy demand, and therefore stipulation of type and layout, on the RAM requirements of such a system for successful computational modelling application. As a benefit of falling costs, computing clusters are becoming more common in medium-sized research facilities and universities. However, they are still generally still beyond the reach of many smaller organizations and individual researchers for a number of reasons. These reasons may revolve around an inability to justify their costs by any sort of sustainable research output, external third-party usage, running costs such as steady electrical power requirements which are likely to be proportional to the size of the system, availability of skilled manpower for configuration and maintenance, programming expertise for task scheduling, interpretation and result integration and so on. The starting inertia in computational research therefore tends to be high for the same reasons. As the popularity of computation continues to rise, the bulk of research output remains largely experimental, leaving computation largely unexplored.

As mentioned above, large MD simulations are today being run regularly for large particle systems with ever-increasing spatiotemporal scales. New innovations or extensions to existing systems are constantly being devised for increased computational performance. Many of these developments revolve around identification of code aspects that can be further parallelized or optimized. Parallelization starts typically with a decomposition, where the problem is divided into smaller “chunks” which are then assigned to a specific processor or core. There are three decomposition methods that are widely used to parallelize MD simulations. These are particle, force and domain decomposition using a spatial basis [25, 26, 27]. The particle number tends to be large in practically meaningful simulations, and none of these decompositions is yet very effective, due to memory consumption, a large number of processors and communication partners. Spatial decomposition approaches, based on physical subdomains, are generally superior [28]. In recent years, under the technological limitations imposed by Moore’s law, increases in computational power were achieved in two ways, namely increasing the number of cores in a processor and the vector size in single-instruction, multiple data instruction set (SIMD) computers. It has become almost de facto in this quest that the only avenue available to increase computational power is by an upward scaling of the number of processors and cores per processor in a chip. For a given configuration of processors, other performance enhancements can be devised though improved algorithms. It is clear that tapping into the full power of such a system will require a highly scalable parallel code. Unfortunately, the secondary issue of power dissipation and age degradation failure (e. g. due to any of the billions of transistors failing because of heat subjection) has to be considered carefully [29, 30]. Domain decomposition can present problems for inhomogeneous systems. The cell-task method [31] can overcome some problems, particularly if used as part of a hybrid decomposition scheme. For instance, by the use of a hybrid algorithm, Mangiardi and Meyer [29] devise a domain decomposition and thread-based hybrid parallelization technique based on large vectorization SIMD processors such as the AVX, AVX-2 and Xeon-Phi processors, comprising several thousands of cores. Their hybrid technique involves modelling with short-ranged forces for both homogeneous and inhomogeneous collections consisting of up to tens of millions of atoms. The hybrid algorithm achieves parallelization by using a task-based approach on smaller sub-domains obtained from domain decomposition. Each sub-domain is handled by a team or “pool” of threads running on multiple cores. The flow of execution of tasks cannot access the same particle simultaneously, thereby averting the need for lockout code to prevent such access of a particle by different tasks. At present, several parallelized MD programs are in existence, for instance [32, 33, 34, 35, 36] and others. However, such programs employ rigid parallelization techniques, and completely new versions will be needed to implement improved parallel algorithms that better exploit advancements in computer architectures [29]. In light of this proliferation, a smaller code that is easier to revise with algorithmic and architectural advancements has clear advantages, and is arguably necessary if advances in architectures are to be exploited fully in future [29, 37]. The effort required for such revisions can be small. There are some issues with hybrid methods, such as synchronization of data structure access which can lead to false sharing of these records because all processors constantly amend the data structures of neighbour cells [38]. These irregular memory access patterns pose significant challenges for parallelization and performance optimization. Additional challenges are associated with the hardware limitations that effectively limit active core performance due to chip configuration, such as memory speed, volume, interchip communication bandwidths, power management, etc. In 2010, Peng et al. [25, 26, 27] proposed three incremental optimizations for an emergent platform, the Godson-T multi-core architecture [39], which was formally released in late 2010. They presented behavioural, instruction-level simulations of a chip that was still in development and aimed at petaflops supercomputing. Their work addressed the foregoing issues of power efficiency, performance, programming and parallelism and platform-operating system portability (particularly to Linux and BSD; portability to WindowsCE is reported) by exploiting polymorphic parallelism. Its main contributions to high-performance computing (HPC) is improved interconnection data communications (both internal and external), fine-grained division of threads, improved thread synchronization and locality awareness. Its unique architecture achieves transfer bandwidths of up to 512 Gb/s for internal register transfers and 51.2 Gb/s for off-chip transfers. It implements 64, 64-bit cores and several levels of memory on chip, specifically 16 address-interleaved 128 kB L2 cache, 32 kB scratch pad memory and four external DDR3 memory controllers. The actual chip was shipped in the optimization strategy employed, termed cell-centred addressing that sequentially maps neighbour 3D cells to 2D cores by a scalar transformation vector which are then searched for migrating atoms in O (1) time. Such mappings are readily solved by classical algorithms but what makes the Godson-T and its progenies unique is the presence of on-chip hardware that implements a locality-aware parallel algorithm for data reuse and latency and congestion avoidance. Such issues typically arise from core to core communications when shared data structures are accessed in L2 cache or off-chip memory during the two-body force calculations. In this contribution, we present a short review of the current state of MD simulations with the aim of achieving speedups on a standalone computer. The discussion is expected to permit a reasonably easy starting point for standalone computational MD for the reader, with the view to readily scale such computers to larger clusters by virtue of interconnection. The work is intended to address single machines that support thread and message passing by virtue of the processor and memory architecture, by using any selection of proliferated tools which are freely available for such purposes. We discuss the rationales for the chosen approach of our own ongoing work in the hope that the generally high initial inertia could be overcome by the interested reader. Finally, we present some of our results thus far and point to an online repository where the evolving code is being maintained. We begin by outlining some standalone architectural configurations that are decidedly parallel or can be made so with relative ease.

2 Standalone architectures

The processor memory architectures of the traditional computer fall into one of two broad categories. In the first category, instruction fetch-execution cycles are repeated from the reset vector up until the last instruction is executed. The instructions and data are stored contiguously in the same, indistinguishable memory blocks. The instruction and data fetches happen along the same paths. In the second category, the instructions are generally stored in different memory locations and the paths to them are different. This sort of architecture is generally faster since simultaneous fetches and pre-fetches are possible. Several enhancements have been developed to support this behaviour, such as pipelining or super pipelining, which involve look-ahead methods and ways to avoid the problems that arise when branches are encountered in the linearly pre-fetched sequence of instructions and data. Such methods include branch prediction, cache memory for the frequently used data thereby saving on fetch times, and so on. A complete discussion of these techniques is beyond the present scope.

3 Getting started with standalone computation

4 Parallel processing paradigms in the C-language

4.1 Threads and message passing

The C-language is native to Unix-based systems, but its high portability has made it standardized and widely available across several platforms. Although there are newer contenders to parallelization, but the C-language is still widely used for parallelization for three main reasons. First, complex data structures that closely resemble collections of structured particles in a system can be defined, referenced and manipulated within memory with considerable ease. The basic data structure typically consists of a particle type object with directly accessible attributes such as mass, position and velocity. Entire collections of such particles, which represent particle configurations within a domain, can be manipulated just as readily as its individual particles through structure access functions. Second, the use of pointers within the language introduces a flexibility and enabling modularity. In this way, individual particles or entire collections of particles can be passed on between different sections the code with ease. An instance where this flexibility is called for is during particle migration. This requires that a particle is added or deleted from a collection. Third, the language is fast and has a number of excellent heap memory management functions. It is native and highly popular on Unix-based systems. The main idea behind thread processing as exploited in parallel programming involves a single problem process giving rise or “spawning” multiple, dynamically managed sub-programs called “threads.” The threads run concurrently and independently and are structured such that they can access shared memory. Shared memory access can involve some interthread communication, requiring careful management of all spawned threads, from the moment each is spawned to the moment it terminates. The process of killing a thread releases its resources for reassignment by the manager. A critical role of the manager is to implement careful synchronization of threads to avoid thread collisions and racing conditions by assuring a mutual exclusivity of memory access during a focus by a particular thread. The concept of threads is well-established in the C-programming language and is conveyed by a number of paradigms. For instance, the POSIX threads standard, or p-threads, makes use of a standardized library header “pthread.h,” which contains the necessary definitions of constants, functions and types. The code in Figure 1 spawns three threads using this approach. Other popular parallelization techniques available to C-programs include the multiprocessing application programming interfaces (APIs), of which the open multiprocessing API (OpenMP) is a popular implementation, and various message passing interfaces (MPIs). The MPIs are libraries that define code and routines which run on a diverse cross-section of platforms including Linux and Windows. The most commonly referenced and freely available MPIs are OpenMPI (not to be confused with OpenMP), MPI, MPICH2 and LAM-MPI. Each of these MPIs may be considered a cross between p-threads and OpenMP.

Figure 1:

An example of p-thread code spawning three threads.

4.2 Open multiprocessing programming

A significant amount of time in complex, repetitive calculations is spent executing loop statements. These iterations lend themselves readily to parallelization using OpenMP by the inclusion of the compiler directive #pragma followed by the code to be parallelized. An example is shown in Figure 2. The same directive can be used to specify the aspects of critical code, i. e. code that is common to all the threads, as shown in Figure 3. The critical code can, for instance, manage a variable that is global relative to the other parallelized code sections. This approach is implicitly serial and has the virtue of eliminating the problems caused by collisions where code sections try to access the same variable simultaneously.

Figure 2:

An example of the OpenMP implementation for parallelization.

Figure 3:

An example of ordinary OpenMP and critical OpenMP processes.

The process of combining the outputs of several smaller parallel sections is known as reduction. In OpenMP, there is a separate keyword that is used with the #pragma omp directive to carry out more reduction efficiently because of its importance. The syntax is shown in Figure 4.

Figure 4:

Syntax of the compiler directive that invokes output reduction of parallel processes.

There is a multitude of other keywords that are associated with the #pragma omp directive, for the streamlining of various multiprocessing tasks.

Figure 5 shows the specific application of the authors to which parallelization was applied. The code snippet is only a small part of a bigger program and shows the calculation of the total potential energy using the Sutton–Chen implementation of the Finnis–Sinclair potential.

Figure 5:

An example of a program parallelized using MPI.

4.5 Cloud virtualization

In very recent years, the concept of the computing Cloud has arose initially as a means of getting large, safe and storage facilities for organization files. This has gradually been progressed to include outright distributed computation and large-data crunching. For instance, financial markets generate vast amounts of time-series data during a trade. The needs of forecasting, fraud detection, internet trade, etc. all require an almost real-time computational element. Depending on the intended functionality and application of the remote user, these powerful new features may be supplied in a limited way free of charge, or at a nominal fee. In any case, this current trend absolves the user from making any major financial outlay in terms of the starting hardware while reaping the benefits of powerful compute facilities and storage. However, security and data privacy and a full grasp of the trappings of this relatively young concept remains a concern for large corporations with obvious interest to safeguard their operations. One recurrent feature of such remote compute facilities is the ability to define virtual environments familiar to the user through virtualization. In essence, several instances of a virtual high-performance computer can be created and loaded by a simple script or graphical interface from a simple less powerful one serving as a terminal from across the planet. The instances can be either homogenous or heterogeneous in terms of the loaded operating system. Examples of organizations with such offerings are Yellow Circle [53], Microsoft Cloud (Azure) [54], Google Cloud platform [55] and others that support free virtualization for individuals in a research university through a secure account. However, the free accounts suffer the limitation in the maximum number of possible virtual instances, a drawback not suffered by paying users. In addition, the technical expertise required to have tasks up and running on such systems with carefully enumerated set {P, C} is essentially nil since server side support exists to tailor applications, particularly for the paying users. In the case of Yellow Circle, for instance, the user can remotely configure the network topology by designing routers, internet protocol addresses, subnets and so on, all through a simple graphical user interface. An image capture of the configuration screen on Yellow Circle with a free account is shown in Figure 6. The tasks of assigning compute nodes, assimilating final output data gathering and task scheduling and general network management are equally straightforward. Consequently, clouds may indeed be the ultimate future of cluster computing. Cloud virtualization can improve the utilization of compute resources while reducing implementations costs. It also introduces a desirable ability to a user to customize their system without many of the technical challenges experienced in traditional grids and cluster configuration. However, at present, virtualization-based cloud has limited performance in its infancy and is yet limited for HPC for a number of reasons. One such limitation can be caused by system virtual machine monitors for supervisory control for the cores and virtual memory peripherals. Methods therefore have to be found to minimize such delays. Ren et al. [56] propose a unique, lightweight supervisory cloud infrastructure called nOSV that serves both HPC and normal applications concurrently. The environment constructs high-performance virtual machines with dedicated resources which operate without supervisory interference, thereby keeping the performance of the traditional Cloud. With their implementations, they report significant performance improvements over other similar virtualization environments.

Figure 6:

A remote virtualization cloud implemented in Yellow Circle for computational purposes.

5 Summary of results

Our own work pertains to face-centred cubic (fcc) atomic structures of metallic atoms, such as copper atoms. The force fields were derived from the embedded atom model adaptation by Sutton–Chen of the Finnis–Sinclair potential. To illustrate the speedup, we calculated the equilibrium lattice parameter (a), the vacancy formation energy of copper, where an atom is removed for just below the surface of the bulk structure and taken to a point several hundred atomic radii away from the surface (to simulate a point at “infinity”), with a cut-off distance of 10 lattice parameters. Similarly, we calculate the extraction energy [7, 9]. The bulk array consists of 30,000 atoms, and computation time is measured using system time functions that are outside the particle setup, the equilibration steps and the final data output steps. The same timer functions are used for the non-parallelized iterations as well as the parallelized code. The configured hardware is the Dell Optiplex 3010 mentioned above. Without parallelization, the run times in a calculation of Ev on an Ubuntu Linux terminal were 6.70 hours and 5.21 hours for the deterministic and Monte-Carlo simulations, respectively. The maximum parallelized speedups are shown in Table 1. The results show a speedup of 1.8 on the wall-clock time for 5 % accuracy on the lattice parameter. The Monte-Carlo method exhibited higher speedup, with the lower acceptance rate, much higher than observed for the deterministic, MD-based simulation. A version of the software can be obtained from [57].

6 Conclusions

This article discusses the considerations to embark on computational research in MD using optimized multi-core computers that can also be operated as a standalone computer, or a node in cluster-based computing. At the onset, a fundamental goal was in the establishment of a computing facility that exploits the features of the single machine, with the secondary aim to extend the functionality to a cluster of such machines, making low-cost HPC possible for the group. The article presents a survey of the current state of the art, on which we base our own ongoing work. We consider what it takes to achieve parallelism and efficiency, providing arguments and pointers to aspects to consider if the need to write own as opposed to using third-party software should arise. Furthermore, we discuss the growing trends in high-performance computation, from the use of GPUs to internet-based virtualization on computing clouds. Finally, we present the summarized results for a sizable array of copper atoms. Both MD and Monte-Carlo methods showed speedup when the iterative steps were parallelized using threads on four cores. A link to a version of the evolving code in a public repository, called the Verlet–Stormer velocity program code [57], is also given for the interested reader. Since threads and MPI are concepts that are standardized and supported in libraries across many platforms, the C-program on which the present work is based could readily be ported to other systems.