Antonio Monari
email: antonio.monari@irsamc.ups-tlse.fr

Anthony Scemama

Michel Caffarel

Large-scale quantum Monte Carlo electronic structure calculations on the EGEE grid

A grid implementation of a massively parallel quantum Monte Carlo (QMC) code on the EGEE grid architecture is discussed. Technical details allowing an efficient implementation are presented and the grid performance (number of queued, running, and executed tasks as a function of time) is discussed. Finally, we present a very accurate Li2 potential energy curve obtained by running simultaneously several hundreds tasks on the grid.

Introduction

Quantum Monte Carlo (QMC) methods are known to be powerful stochastic approaches for solving the Schrödinger equation. [START_REF] Foulkes | Quantum Monte Carlo simulations of solids[END_REF] Although they have been widely used in computational physics during the last twenty years, they are still of marginal use in computational chemistry. [START_REF] Hammond | Monte Carlo Methods in Ab Initio Quantum Chemistry[END_REF] Two major reasons can be invoked for that: i) the N -body problem encountered in chemistry is particularly challenging (a set of strongly interacting electrons in the field of highly-attractive nuclei) ii.) the level of numerical accuracy required is very high (the so-called "chemical accuracy"). In computational chemistry, the two standard approaches used presently are the Density Functional Theory (DFT) approaches and the various post-Hartree-Fock wavefunction-based methods (Configuration Interaction, Coupled Cluster, etc.) In practice, DFT methods are the most popular approaches, essentially because they combine both a reasonable accuracy and a favorable scaling of the computational effort as a function of the number of electrons. On the other hand, post-HF methods are also employed since they lead to a greater and much controlled accuracy than DFT. Unfortunately, the price to pay for such an accuracy is too high to be of practical use for large molecular systems.

QMC appears as a third promising alternative method essentially because it combines the advantages of both approaches: a favorable scaling together with a very good accuracy. In addition to this, and it is the central point of the present note, the QMC approaches -in sharp contrast with DFT and post-HF methodsare ideally suited to High-Performance-Computing (HPC) and, more specifically, to massive parallel computations either on homogeneous multiprocessor platforms or on heterogeneous grid infrastructures. As most "classical" or "quantum" Monte Carlo approaches, the algorithm is essentially of the number crunching type, the central memory requirements remain small and bounded and the I/O flows are essentially marginal. Due to these extremely favorable computational aspects plus the rapid evolution of computational infrastructures towards more and more numerous and efficient processors, it is likely that QMC will play in the next years a growing role in computational chemistry.

In the present study, the first implementation of our quantum Monte Carlo program on a large scale grid -the European EGEE-III grid[3]-is presented. As a scientific application we have chosen to compute with a very high accuracy the potential energy curve (PEC) of the Li 2 molecule (total energy of the system as a function of the Li-Li distance). To the best of our knowlegde, the curve presented here is the most accurate PEC ever published for this system. In order to reach such an accuracy two conditions need to be fulfilled. First, a large enough Monte Carlo statistics has to be realized to reduce the final statistical error down to the precision desired. Second, accurate enough trial wave functions must be employed to reduce as much as possible the so-called "Fixed-Node" error (the only systematic error left in a QMC calculation, see Ref. [START_REF] Hammond | Monte Carlo Methods in Ab Initio Quantum Chemistry[END_REF]). The first condition is easy to fulfill since the system is small (only, six electrons) and accumulating statistics is just a matter of making enough Monte Carlo steps and using enough processors ("brute force" approach). The second condition is much more challenging since we need to introduce trial wavefunctions with a controlled nodal quality and which can be improved in a systematic way. Here, we have realized this latter aspect by considering wavefunctions issued from Full Configuration Interaction (FCI) calculations in a large basis-set (technically, the cc-pVQZ basis set, Ref. [START_REF] Dunning | Gaussian basis sets for use in correlated molecular calculations. I) The atoms boron through neon and hydrogen[END_REF]). Such a FCI trial wavefunction is expected to have a very good nodal structure. However, there is a price to pay: To handle such a function is quite expensive. More precisely, the FCI trial wavefunction used here is expressed a sum of 16,138 products of two 3×3 determinants (three α-electrons and three β-electrons) and, at each Monte Carlo step, this wavefunction and its first-and second-derivatives have to be computed. Note that the computational cost in terms of CPU time is directly proportional to the number of products in the trial wavefunction expansion. To the best of our knowledge, it is the first time that such a high number of determinants in a QMC calculation has been used. Let us emphasize that it has been possible here only because of the use of the grid infrastructure.

In Sec. 2 some technical details related to the implementation of a quantum Monte Carlo simulation and the use of our QMC=Chem [START_REF]QMC=Chem is a general-purpose quantum Monte Carlo code for electronic structure calculations[END_REF] program are presented. Section 3 presents the computational strategy employed in our application to the Li 2 molecule. Section 4 gives the results and discusses the performance. Finally, some conclusions are presented in Sec.5.

Technical details

A walker is a vector X of the 3N -dimensional space containing the entire set of the three-dimensional coordinates of the N electrons of the system. During the simulation, a walker (or a population of walkers) samples via a Monte Carlo Markov Chain process the 3N -dimensional space according to some target probability density (the precise density may vary from one QMC method to another). From a practical point of view, the averages of the quantities of interest (energy, densities, etc.) are calculated over a set as large as possible of independent random walks. Random walks differ from each other only in the initial electron positions X 0 , and in the initial random seed S 0 determining the entire series of random numbers used.

In the QMC=Chem code used here, the main computational object is a block. In a block, N walk independent walkers realize random walks of length N step , and the quantities of interest are averaged over all the steps of each random walk. If N step is significantly larger than the auto-correlation time (which is usually rather small), the positions of the walkers at the end of the block can be considered as independent of their initial positions and a new block can be sampled using these configurations as X 0 and using the current random seed as S 0 .

The final Monte Carlo result is obtained by averaging all the results obtained for each block. If the data associated with each block are saved on disk, the averages can be calculated as a post-processing of the data and the calculation can be easily restarted using the last positions of the walkers and the last random seed.

Note that the computation of the averages does not require any time ordering. If the user provides a set of N proc different initial conditions (walker positions and random seed), the blocks can be computed in parallel. In figure 1, we give a pictorial representation of four independent processors computing blocks sequentially, each block having different initial conditions.

Design of the QMC=Chem program

The QMC=Chem program was designed specifically to run on heterogeneous clusters via the Message Passing Interface (MPI) API [START_REF] Gropp | A high-performance, portable implementation of the MPI message passing interface standard Parallel Computing[END_REF] and also in grid environments via Python[8] scripts. The memory requirements, disk input/outputs and network communications were minimized as much as possible, and the code was written in order to allow asynchronous processes. This section presents the general design of the program.

The behavior of the program is the following. A main Python program spawns three processes: an observer, a computation engine, and a data server (see figure 2).

The observer The observer keeps a global knowledge of the whole calculation (current values of the computed averages, total CPU time, wall time, etc). It updates the results using the calculated blocks at regular intervals of time and checks if the code should continue or stop by informing the data server. It also checks if a stopping condition is reached. The stopping condition can be a maximum value of the total CPU time, the wall time, the number of blocks, or a threshold on the statistical error bar of any Monte Carlo estimate. The computation engine The computation engine starts the Fortran MPI executable. The master MPI process broadcasts the common input data to the slaves, and enters the main loop. In the main loop, the program computes one block and sends the results to the data server via a trivial python XML-RPC client. The reply of the data server determines if the main loop should be exited or if it should compute another block. When the main loop is exited, there is an MPI synchronization barrier where all the slave processes send the last walker positions and their last random seed to the master, which writes them to disk.

A linear feedback shift register (LFSR) pseudo-random number generator [START_REF] Ecuyer | Tables of maximally equidistributed combined LFSR generators[END_REF] is implemented in the code. A pool of 7000 initial random seeds was previously prepared, each random seed being separated from the previous one by 6.10 11 seeds. Every time a random number is drawn, a counter is incremented. If the counter reaches 6.10 11 , the next free random seed is used. This mechanism guarantees that the parallel processes will never use the same sequence of random numbers.

The data server The data server is a Python XML-RPC server whose role is to receive the computed data during a simulation and save it into files. Each file is a few kilobytes large and contains the averages of interest computed over a block. The data server also computes an MD5 key [START_REF] Rivest | Technical report. Internet Activities Board[END_REF] related to the critical input values to guarantee that the computed blocks belong to the simulation, and that the input data has not been corrupted.

Generation of the output

The output file of the program is not created during the run which only produces block files via the data server. A separate script analyzes the blocks written to disk to produce an output. This script can be executed at any time: while a calculation is running to check the current values, or when the simulation has finished. The consistence between the input data with the blocks is checked with using the previously mentioned MD5 key.

Adaptation to grid environments A script was written to prepare as much files as there are different requested tasks on the grid, each file name containing the index of the task (for example 12.tar.gz). Each file is a gzipped tar file of a directory containing all the needed input data. The only difference between the input files of two distinct tasks is their vector X 0 and their random seed S 0 . The script also generates a gzipped tar file qmcchem.tar.gz which contains a statically linked non-MPI Fortran i686-linux executable, and all the needed Python files (the XML-RPC data server, the observer, etc). A third generated file is a shell script qmcchem grid.sh which will run the job.

Another script reverses this work. First, it unpacks all the tar.gz files containing the input files and the block files. Then it collects all the block files into a common directory for the future production of the output file.

Termination system signals (SIGKILL, SIGTERM, . . . ) are intercepted by the QMC=Chem program. If any of this signals is caught, the program tries to finish the current block and terminates in a clean way.

Advantages of such a design

Asynchronous processes Note that the transmission of the computed data is not realized through MPI in the main loop of the computation engine, but with a Python script instead. This latter point avoids the need for MPI synchronization barriers inside the main loop, and allows the code to run on heterogeneous clusters with a minimal use of MPI statements. The main advantage of this design is that if the MPI processes are sent on machines with different types of processors, each process will always use 100% of the CPU (except for the synchronization barrier at the end) and fast processors will send more blocks to the data server than the slower processors.

Analysis of the results

The analysis of the blocks as a post-processing step (see section 2.1) has major advantages. First, the analysis of the data with graphical user interfaces is trivial since all the raw data is present in the block files, and they are easy to read by programs, as opposed to traditional output files which are written for the users. The degree of verbosity of the output can be changed upon request by the user even after the end of the calculation, and this avoids the user to read a large file to find only one small piece of information, while it is still possible to have access to the verbose output. The last and most important feature is that the production of the output does not impose any synchronization of the parallel processes, and they can run naturally in grid environments.

Failure of a process As all the processes are independent, if one task dies it does not affect the other tasks. For instance, if a large job is sent on the grid and one machine of the grid has a power failure, the user may not even remark that part of the work has not been computed. Moreover, as the process signals are intercepted, if a batch queuing system tries to kill a task (because it has exceeded the maximum wall time of the queue, for example), the job is likely to end gracefully and the computed blocks will be saved. This fail-safe feature is essential in a grid environment where it is almost impredictable that all the requested tasks will end as expected.

Flexibility The duration of a block can be tuned by the user since it is proportional to the number of steps per trajectory. In the present work, the stopping condition was chosen to be wall time limit, which is convenient in grid environments with various types of processors. When the stopping condition is reached, if the same job is sent again to the queuing system, it will automatically continue using the last walker positions and random seeds, and use the previously computed blocks to calculate the running averages. Note that between two subsequent simulations, there is no constraint to use the same number or the same types of processors.

Computational Strategy and Details

A quantum Monte Carlo study has been performed on the Li 2 molecule. The choice of such a system allowed us to consider a Full Configuration Interaction (FCI) wave-function as the trial wave function, i.e. a virtually exact solution of the Schrödinger equation in the subspace spanned by the gaussian orbital basis. We remind that the computational cost of the FCI problem scales combinatorially with the number of basis functions and electrons and, therefore high quality FCI wave functions are practically impossible to be obtained for much larger systems.

The Q5Cost common data format

Due to the inherent heterogeneity of grid architectures, and due to the necessity of using different codes, a common format for data interchange and interoperability is mandatory in the context of distributed computation. For this reason we have previously developped a specific data format and library for quantum chemistry [START_REF] Scemama | Common Format for Quantum Chemistry Interoperability: Q5Cost format and library[END_REF], and its use for single processor and distributed calculations has already been reported [START_REF] Vetere | A Theoretical Study of Linear Beryllium Chains: Full Configuration Interaction[END_REF]. The Q5Cost is based on the HDF5 format, a characteristic that makes the binary files portable on a multiple platform environment. Moreover the compression features of the HDF5 format are exploited to reduce significantly the file size while keeping all the relevant information and metadata. Q5Cost contains chemical objects related data organized in a hierarchical structure within a logical containment relationship. Moreover a library to write and access Q5Cost files has been released [START_REF] Scemama | Common Format for Quantum Chemistry Interoperability: Q5Cost format and library[END_REF]. The library, built on top of the HDF5 API, makes use of chemical concepts to access the different file objects. This feature makes the inclusion on quantum chemistry codes rather simple and straightforward, leaving the HDF5 low level technical details absolutely transparent to the chemical software developer. Q5Cost has emerged as an efficient tool to facilitate communication and interoperability and seems to be particularly useful in the case of distributed environments, and therefore well adapted to the grid.

Computational Details

All the preliminary FCI calculations have been realized on a single processor by using the Bologna FCI code. [START_REF] Bendazzoli | A vector and parallel full configuration interaction algorithm[END_REF], [START_REF] Gagliardi | Direct-list algorithm for configuration interaction calculations[END_REF] The code has been interfaced with Molcas [START_REF] Karlström | MOLCAS: a program package for computational chemistry[END_REF] to get the necessary one-and two-electron molecular integrals. The FCI computations considered in the present work involved up to 16,138 symmetry adapted and partially spin adapted determinants. All the communications between the different codes has been assured by using the Q5Cost format and library [START_REF] Scemama | Common Format for Quantum Chemistry Interoperability: Q5Cost format and library[END_REF]. In particular a module has been added to the Molcas code to produce a Q5Cost file containing the information on the molecular system and the atomic and molecular (self consistent field level) integrals. The Q5Cost file has been directly read by the FCI code, and the final FCI wave function has been added to the same file in a proper and standardized way. The actual QMC=Chem input has been prepared by a Python script reading the Q5Cost file content. Before running the QMC calculation on the grid an equilibration step was performed (i.e., building "good" starting configurations, X 0 , for walkers) by doing a quick variational QMC run (see, Refs. [START_REF] Foulkes | Quantum Monte Carlo simulations of solids[END_REF] and [START_REF] Hammond | Monte Carlo Methods in Ab Initio Quantum Chemistry[END_REF]) in single processor mode.

QMC computations have been run on the EGEE grid over different computing elements and in a massively parallel way. Typically for each potential energy curve point we requested to use 1000 nodes, obtaining at least about 500 tasks running concurrently. Once the job on each node was completed the results were retrieved and the output file was produced by the post-processing script to obtain the averaged QMC energy. Due to the inherent flexibility of the QMC implementation the fact of having different tasks on different nodes terminating after different number of blocks did not cause any difficulty as the output file was produced independently from the computation phase. Moreover, the failure or the abortion of some tasks did not impact significantly on the quality of the results.

Results and Grid Performance

Our results for the potential energy curve of the Li 2 molecule are presented in Figure 3 (graphical form) and Table 1 (raw data and error bars). Results are given for a set of 31 inter-nuclear distances. Let us emphasize that the data are of a very high-quality and the energy curve presented is, to the best of our knowledge, the most accurate energy curve ever published. To illustrate this point, let us mention that the dissociation energy defined as D e ≡ E(R eq ) -E(R = ∞) is found to be here D e = -0.0395(2) a.u. in excellent agreement with the experimental result of D e = -0.03928 a.u., Ref. [START_REF] Filipp | Multiconfiguration wave functions for quantum Monte Carlo calculations of firstrow diatomic molecules[END_REF]. A much more detailed discussion of these results and additional data will be presented elsewhere [START_REF] Monari | A quantum Monte Carlo study of the fixednode Li2 potential energy curve using Full Configuration Interaction nodes[END_REF].

For this application only 11Mb of memory per Fortran process was needed to allow the computation of the blocks. This feature is twofold. First, it allowed our jobs to be selected rapidly by the batch queuing systems as very few resources were requested. Second, as the memory requirements are very low, the code was able to run on any kind of machine. This allowed our jobs to enter both the 32-bit (x86) and 64-bit (x86 64) queues.

Grid Performance

A typical analysis of the grid performance can be seen on Figure 4. Here we report the number of queued, running and executed tasks as a function of time for the submission of a parametric job of 1000 tasks.

One can see that after a very limited amount of time (less than one hour) almost half of the tasks were in a running status. Correspondingly the number of queued tasks undergoes a very rapid decay, indicating that a consistent percentage of the submitted tasks did not spend a significant amount of time in the queue. This feature is a consequence of the very limited amount of resources requested by the job, and by the fact that virtually all the queues could be used, and should therefore be ascribed to the high flexibility of our approach. Correspondingly the number of running tasks experiences a maximum at about one hour. It is also important to notice the high asymetry of the peak, indicating that while the maximum amount of running tasks is achieved quite rapidly, the high efficiency (number of running tasks) is mantained for a considerable amount of time before a degradation. The number of completed tasks too experiences a quite rapid increase. After about four hours the numbers of running tasks remains constant and quite low, and consequently a small variation is observed on the number of completed tasks. The decayng of the number of running job is due to the fact that a vast majority of tasks have been completed achieving the completion of the desired number of blocks. The remaining jobs (about 50) can be abscribed to tasks running on very slow processors and so requiring a greater amount of time to complete the blocks, moreover in some cases some of the jobs could be considered as "ghost" i.e. labelled as running although being stalled. After 9 hours the job was cancelled since the desired precision level had been achieved, even if only about 700-800 tasks were actually completed. The overall process required on average about 1000 CPU hours for each geometry equivalent to about 40 days of CPU time. It is also noteworthy to comment the distribution of the CPU time per block for each task (Figure 5), this difference reflects the fact that some tasks have been performed on slower nodes. Again this asymetry did not influence the final result since the final statistics were performed on the total amount of blocks.

Conclusions

An efficient grid implementation of a massively parallel quantum Monte Carlo code has been presented. The strategy employed has been presented in detail. Some test applications have been performed on the EGEE-III grid architecture showing the efficiency and flexibility of our implementation. As shown, our approach enables to exploit the computational power and resources of the grid for the solution of a non-trivial N -body problem of chemistry. It must be emphasized that the very important computational gain obtained here has concerned the simulation of the Shrödinger equation for a single fixed nuclear geometry (single point energy calculation). This is in sharp contrast with the common situation in computational chemistry where parallelism is not used (or partially used) for solving the problem at hand (algorithms are in general poorly parallelized) but rather for making independent simulations at different nuclear geometries (trivial parallelization based on different inputs). We believe that the quantum Monte Carlo approach which is based on Markov chain processes and on the accumulation of statistics for independent events can represent an ideal test bed for the use of grid environments in computational chemistry. Finally, we note that the combination of grid computing power and of the QMC ability to treat chemical problems at a high-level of accuracy can open the way to the possibility of studying fascinating problems (from the domain of nano-sciences

Fig. 1 .

 1 Fig. 1. Graphical representation of a QMC simulation. Each process generates blocks, each block being composed of N walk walkers realizing Nstep Monte Carlo steps.

Fig. 2 .

 2 Fig. 2. Inter-process communication of the QMC=Chem program.

Fig. 4 .

 4 Fig. 3. The QMC Li2 potential energy curve

Fig. 5 .

 5 Fig. 5. Histogram of the CPU time per block. This figure shows the heterogeneity of the used processors.

Table 1 .

 1 Quantum Monte Carlo energies (atomic units) as a function of the Li-Li distance (atomic units). Values in parenthesis correspond to the statistical error on the two last digits.

	Distance (a.u.) Energy (a.u.) Distance (a.u.) Energy (a.u.)
	2.2	-14.81854(44)	5.2	-14.99491(29)
	2.4	-14.85192(36)	5.4	-14.99431(30)
	2.6	-14.87861(36)	5.6	-14.99279(30)
	2.8	-14.90269(37)	5.8	-14.99090(30)
	3.0	-14.92195(35)	6.0	-14.99018(32)
	3.2	-14.93861(38)	6.4	-14.98541(26)
	3.4	-14.95324(37)	6.8	-14.98088(29)
	3.6	-14.96428(38)	7.2	-14.97804(29)
	3.8	-14.97320(38)	7.6	-14.97281(28)
	4.0	-14.98250(29)	8.0	-14.96984(28)
	4.2	-14.98629(27)	10.0	-14.95951(27)
	4.4	-14.99016(30)	12.0	-14.95747(15)
	4.6	-14.99285(29)	14.0	-14.95624(12)
	4.8	-14.99358(29)	16.0	-14.95606(14)
	5.0	-14.99479(29)	18.0	-14.95569(17)
	5.051	-14.99492(17)	. . .	. . .
			100.0	-14.95539(11)

Acknowledgements

Support from the French CNRS and from University of Toulouse is gratefully acknowledged. The authors also wish to thank the French PICS action 4263 and European COST in Chemistry D37 "GridChem". Acknowledgments are also due to the EGEE-III grid organization and to the CompChem virtual organization.

to biological systems) which are presently out of reach.


Anthony Scemama

Michel Caffarel

Emmanuel Oseret

William Jalby

QMC=Chem: A Quantum Monte Carlo Program for Large-Scale Simulations in Chemistry at the Petascale Level and beyond

In this work we discuss several key aspects for an efficient implementation and deployment of large-scale quantum Monte Carlo (QMC) simulations for chemical applications on petaflops infrastructures. Such aspects have been implemented in the QMC=Chem code developed at Toulouse (France). First, a simple, general, and fault-tolerant simulation environment adapted to QMC algorithms is presented. Second, we present a study of the parallel efficiency of the QMC=Chem code on the Curie machine (TGCC-GENCI, CEA France) showing that a very good scalability can be maintained up to the more than 70 000 cores. Third, it is shown that a great enhancement in performance with the single-core optimization tools developed at Versailles (France) can be obtained.

Introduction

Quantum Monte Carlo (QMC) is a generic name for a large class of stochastic approaches solving the Schrödinger equation by using random walks. In the last forty years they have been extensively used in several fields of computational physics and are in most cases considered as state-of-the-art approaches. However, this is not yet the case in the important field of computational chemistry where the two "classical" computational methods are the Density Functional Theory (DFT) and post-Hartree-Fock methods. For a review of QMC and its status with respect to the standard approaches, see e.g. [START_REF] Caffarel | Quantum monte carlo methods in chemistry[END_REF]. In the recent years several applications for realistic chemical problems have clearly demonstrated that QMC has a high potential in terms of accuracy and in ability of dealing with (very) large systems. However, and this is probably the major present bottleneck of QMC, simulations turn out to very CPU-expensive. The basic reason is that chemical applications are particularly demanding in terms of precision: the energy variations involved in a chemical process are typically several orders of magnitude smaller than the total energy of the system which is the quantity computed with QMC. Accordingly, the target relative errors are typically 10 -7 or less and the Monte Carlo statistics needed to reach such a chemical accuracy can be tremendously large. Now, the key point for the future is that this difficulty is expected to be largely overcome by taking advantage of the remarkable property of QMC methods (not valid for standard methods of chemistry) of being ideally suited to HPC and, particularly, to massive parallel computations. In view of the formidable development of computational platforms this unique property could become in the near future a definite advantage for QMC over standard approaches.

The stochastic nature of the algorithms involved in QMC enables to take advantage of today and tomorrow's computer architectures through the following aspects: i) Data structures are small inducing a fairly small memory footprint (less than 300 MiB per core for very large systems) and data accesses are organized to maximize cache usage: spatial locality (stride-one access) and temporal locality (data reuse), ii) Most of the computation can be efficiently vectorized making full use of the vector capabilities of recent processors, iii) Network communications can be made non-blocking, iv) Access to persistent storage is negligible and can be made non-blocking, v) Different parallel tasks can be made independent of each other so as to run asynchronously, vi) Fault-tolerance can be naturally implemented. All these features which have been implemented in the QMC=Chem code developed at Toulouse [START_REF]Quantum Monte Carlo for Chemistry@Toulouse[END_REF] are required to take advantage of large-scale computing grids [START_REF] Monari | Large-scale quantum monte carlo electronic structure calculations on the egee grid[END_REF] and to achieve a very good parallel efficiency on petascale machines.

In section 2 a short overview of the mathematical foundations of the QMC method employed here is presented. For a more detailed presentation the reader is referred to [START_REF] Caffarel | Quantum monte carlo methods in chemistry[END_REF] and references therein. Section 3 is devoted to the presentation of the general structure of the simulation environment employed for running QMC=Chem on an arbitrary computational platform. A preliminary version of such an environment has been presented in [START_REF] Monari | Large-scale quantum monte carlo electronic structure calculations on the egee grid[END_REF]. Here, we present an improved implementation where the network communications are now fully handled by a client-server implementation and the computational part is isolated in multiple instances of a single-core Fortran program. These modifications allow the program to survive failures of some compute nodes. In section 4 the results of our study of the parallel efficiency of the QMC=Chem code performed on the Curie machine (TGCC-GENCI, France) thanks to a PRACE preparatory access [START_REF] Caffarel | Large-scale quantum monte carlo simulations for chemistry[END_REF] are presented. In section 5 the results of the optimization of the single-core performance are discussed. The optimization was performed after a static assembly analysis of the program and a decremental analysis.

Overview of a QMC Simulation

In the simulations discussed here, the basic idea is to define in the 3N-dimensional electronic configuration space a suitable Monte Carlo Markov chain combined with a birth-death (branching) process to sample a target probability density from which exact (or high-quality) quantum averages can be evaluated. In our simulations we employ a variation of the Fixed-Node Diffusion Monte Carlo (FN-DMC) method, one of the most popular versions of QMC. In short, we aim at solving the electronic Schrödinger equation written as

H Ψ 0 (r 1 , . . . , r N ) = E 0 Ψ 0 (r 1 , . . . , r N ) ( 1
)
where H is the molecular Hamiltonian operator, Ψ 0 (r 1 , . . . , r N ) the N-electron wave function, and E 0 the total electronic energy.

To do that, the basic idea is to construct a stochastic process having the density π(r 1 , . . . , r N ) = Ψ 0 (r 1 , . . . , r N )Ψ T (r 1 , . . . , r N ) . . . dr 1 . . . dr N Ψ 0 (r 1 , . . . , r N )Ψ T (r 1 , . . . , r N )

(2) as stationary density. Here, Ψ T -called the trial wavefunction -is some good known (computable) approximation of Ψ 0 . The role played by the trial wavefunction is central since it is used to implement the "importance sampling" idea, a fundamental point at the heart of any efficient Monte Carlo sampling of a high-dimensional space. In the important case of the total energy, it can be shown that the exact ground-state energy E 0 may be expressed as the average of the so-called local energy defined as E L (r 1 , . . . , r N ) ≡ H Ψ T Ψ T over the density π. In brief, the stochastic rules employed are as follows:

N step N proc N walk

CPU time

1. Use of a standard Markov chain Monte carlo chain based on a drifted Brownian motion. The role of the drift term is to push the configurations (or "walkers" in the QMC terminology) towards the regions where the trial wavefunction takes its largest values (importance sampling technique). 2. Use of a birth-death (branching) process: the walkers are killed or duplicated a certain number of times according to the magnitude of the local energy (low values of the local energy are privileged).

It can be shown that by iterating these two rules for a population of walkers the stationary density π (Eq. 2) is obtained. Note that in the actual implementation in QMC=Chem, a variation of the FN-DMC method working with a fixed number of walkers is used [START_REF] Assaraf | Diffusion monte carlo methods with a fixed number of walkers[END_REF].

Denoting as X = (r 1 , . . . , r N ) a walker in the 3N-dimensional space, the random trajectories of the walkers differ from each other only in the initial electron positions X 0 , and in the initial random seed S 0 determining the entire series of pseudo-random numbers.

The main computational object is a block. In a block, N walk independent walkers realize random walks of length N step , and the energy is averaged over all the steps of each random walk. N step is set by the user, but has to be taken large enough such that the positions of the walkers at the end of the block can be considered independent from their initial positions. A new block can be sampled using the final walker positions as X 0 and using the current random seed as S 0 . The block-averaged energies are Gaussian distributed and the statistical properties can be easily computed. The final Monte Carlo result is obtained by super-averaging all the block-averages. If the block-averages are saved to disk, the final average can be calculated by post-processing the data and the calculation can be easily restarted at any time. As all blocks are completely independent, different blocks can be computed asynchronously on different CPU cores, different compute nodes and even in different data centers. Figure 1 shows a pictorial representation of three independent CPU cores computing blocks sequentially, each block having different initial conditions.

The core of QMC=Chem is a single-core Fortran executable that computes blocks as long as a termination event has not been received. At the end of each block the results are sent in a non-blocking way to a central server, as described in the next section.

The Client/Server Layer

In the usual MPI implementations, the whole run is killed when one parallel task will not be able to reach the MPI Finalize statement. This situation occurs when a parallel task is killed, often due to a system failure. For deterministic calculations where the result of every parallel task is required, this mechanism is convenient since it immediately stops a calculation that will not give the correct result. In our case, as the result of the calculation of a block is a Gaussian random variable, removing the result of a block from the simulation is not a problem since doing that does not introduce any bias in the final result. Therefore, if one compute node fails, the rest of the simulation should survive.

A second disadvantage of using MPI libraries is that all resources need to be available for a simulation to start. In our implementation, as the blocks can be computed asynchronously we prefer to be able to use a variable number of cores during the simulation in order to reduce the waiting time in the batch queue.

These two main drawbacks lead us to write a lightweight TCP client/server layer in the Python language to handle all the network communications of the program and the storage of the results in a database. The Python program is divided into three distinct tasks shown in figure 2, the first and second tasks running only on the master compute node.

The first task is a manager that watches periodically the database associated with the running simulation. The manager computes the running averages and error bars, and checks if the stopping condition of the calculation is reached. The stopping condition can be for instance a threshold on the error bar of the energy, a condition on the total execution time, etc. If the stopping condition is reached, a stopping flag is set in the database.

The second task is a data server. This task contains two threads: one network thread and one I/O thread. The network thread periodically sends a UDP packet to the connected clients to update their stopping condition and to check that they are still running. The network thread also receives the block averages and puts them in a queue. Simulta- neously, the I/O thread empties this queue by storing in the database the block averages. At any time, a new client can connect to the data server to request the input files and participate to a running calculation. Another way to increase the number of running cores is to submit another run reading and writing to the same database. As the managers read periodically the content of the database, each simulation is aware of the results obtained by the other simulations. This allows the use of multiple managers and data servers that can be submitted to the batch scheduler as independent jobs in order to gather more and more computing resources as they become available.

The third task is a forwarder. Each compute node (as well as the master node) has only one instance of the forwarder. The forwarder has different goals. The first goal is to spawn the computing processes (single-core Fortran executables) and to collect the results via Unix pipes after a block has been computed. Then, it sends the results to the data server while the computing processes are already computing the next block. If every compute node sends directly the results to the data server, the master node is flooded by small packets coming from numerous sources. Instead the forwarders are organized in a binary tree structure, and the second goal of a forwarder is to collect results from other forwarders and transfer them in a larger packet to its parent in the binary tree. Using this structure, the data server has much fewer connected clients, and receives much larger packets of data. All the nodes of the forwarder tree can possibly connect to all their ancestors if the parent forwarder does not respond.

On massively parallel machines an MPI launcher is used to facilitate the initialization step. The launcher sends, via the MPI library, the input files and Python scripts to all the slave nodes allocated by the batch scheduler. The files are written in a RAM disk on every node (the /dev/shm location). The reason for this copy is to avoid too many simultaneous I/O on the shared file system, and also to avoid I/O errors if the shared file system fails, if a local disk fails or is full. The MPI launcher then forks to an instance of a forwarder that connects to the data server.

Parallel Efficiency

Using the design described in the previous section, the parallel section of the program is expected to display a parallel efficiency of about 100% since there is no blocking statement. In this section, we investigate in some detail how the initialization and finalization steps impact the parallel efficiency of QMC=Chem.

For that, a small benchmark was set up during a PRACE preparatory access on the Curie machine (TGCC-GENCI, France). [START_REF] Caffarel | Large-scale quantum monte carlo simulations for chemistry[END_REF] On this machine, all the compute nodes were equipped with the same processors, namely four Intel Nehalem 8-core sockets. Two important parameters are to be kept in mind. First, we chose to run a short simulation for which the average CPU time required to compute one block is 82 seconds. 3Second, we have chosen to send the stopping signal after 300 seconds (during the fourth block). When the forwarders receive the stopping signal, they wait until all the working CPU cores finish their current block. Hence, in this study the ideal wall time for perfect scalability should be 328 seconds.

The additional time T (N core ) with respect to the ideal time can be expressed as

T (N core ) = W (N core ) - C(N core ) N core ( 3
)
where N core is the number of cores, W (N core ) is the wall time and C(N core ) is the CPU time, both measured for N core cores. With 10 000 cores, 149 seconds are needed for the initialization and finalization steps. For this 7 minutes benchmark, a parallel efficiency of 69% was obtained. However, as the parallel section has an ideal scaling, one can extrapolate the parallel efficiency one would obtain for a one hour run. If the stopping signal occurs after one hour, each core would have computed 44 blocks. The total CPU time would be C(N core ) = 11C(N core ). As the additional time T (N core ) does not depend on the number of computed blocks, the wall time would be W (N core ) = T (N core ) + C(N core )/N core . A parallel efficiency of 96% would be obtained for a one-hour run on 10 000 CPU cores (figure 3).

More recently, we were given the opportunity to test QMC=Chem on the thin nodes of the Curie machine (80 000 Sandy Bridge cores), for a real application on a biological molecule made of 122 atoms and 434 electrons (the largest application ever realized using all-electron Diffusion Monte Carlo). Using 51 200 cores for 3 hours, the parallel efficiency was 79%. After the runs were finished, we realized that for such a large molecular system, the CPU time needed to compute one block had quite large fluctuations due to the implementation the dense-sparse matrix product presented in the next section. This implementation considerably reduces the total wall time (which is what the end user wants), but slightly reduce the parallel efficiency. This problem has been solved by making the number of Monte Carlo steps per block non-constant. Nevertheless, these runs confirm that a good scaling can still be obtained for a real simulation.

Single-Core Efficiency

Our choice in the implementation of the QMC algorithms was to minimize the memory footprint. This choice is justified first by the fact that today the amount of memory per CPU core tends to decrease and second by the fact that small memory footprints allows in general a more efficient usage of caches. Today, the standard size of the molecular systems studied by QMC methods and published in the literature usually comprise less than 150 electrons. For a 158 electron simulation, the binary memory footprint (including code and data) per core is only 9.8 MiB. To check the memory footprint of much larger systems, a few Monte Carlo steps were performed successfully on a molecular system containing 1731 electrons; such a large system only required 313 MiB of memory per core. For a system beyond the largest systems ever computed with all-electron QMC methods, the key limiting factor is only the available CPU time and neither the memory nor disk space requirements. This feature is well aligned with the current trends in computer architecture for large HPC systems.

As the parallel scaling is very good, single-core optimization is of primary importance: the gain in execution time obtained on the single-core executable will also be effective in the total parallel simulation. The Fortran binary was profiled using standard profiling tools (namely gprof [START_REF]The GNU profiler[END_REF] and Scalasca [START_REF] Wolf | Scalasca. In Encyclopedia of Parallel Computing[END_REF]). Both tools exhibit two major hot spots in the calculation of a Monte Carlo step. The first hot spot is a matrix inversion, and the second hot spot is the product of a constant matrix A with five different matrices B 1 . . . B 5 . These two bottlenecks have been carefully optimized for the x86 micro-architectures, especially the for the AVX instruction set of the Sandy Bridge processors.

To measure the performance of the matrix inversion and the matrix products, small codelets were written. The final results are given in table 1, compared to the performance of the single core executable, with different molecular system sizes. Computational complexity (with respect to FP operations) of the matrix inversion is O(N 3 e ) where N e is the number of electrons. Exploiting the sparse structure of the right matrices in the matrix matrix products, computational complexity of such products is reduced to O(N 2 e ) with a prefactor depending on N basis , the size of the basis set used to describe the wave function.

The matrix inversion is performed in double precision (DP) using the Intel MKL library, an implementation of the LAPACK [START_REF] Anderson | Lapack: a portable linear algebra library for highperformance computers[END_REF] and BLAS [START_REF] Dongarra | An extended set of fortran basic linear algebra subprograms[END_REF] APIs. To maximize MKL efficiency, arrays were padded to optimize array alignment (lined up on 256 bit boundaries) and the leading dimension of the array is chosen to be a multiple of 256 bits to ensure that all the column accesses are in turn properly aligned.

For the matrix matrix products, similar alignment/padding techniques were used. Loops were rearranged in order to use full vector length stride-one access on the left dense matrix and then blocked to optimize temporal locality (i.e. cache usage).

An x86 64 version of the MAQAO framework [START_REF] Djoudi | MAQAO: Modular assembler quality Analyzer and Optimizer for Itanium 2[END_REF] was used to analyze the binary code and to generate best possible static performance estimates. This technique was used not only on the matrix matrix products but also on all of the hottest loops (i.e. accounting at least for more than 1% of the total execution time). This allowed us to detect a few compiler inefficiences and to fix them by hard coding loop bounds and adding up pragmas (essentially for allowing use of vector aligned instructions).

Then, the DECAN tool [START_REF] Koliai | A balanced approach to application performance tuning[END_REF] was used to analyze performance impact of data access. For that purpose, for each loop, two modified binaries were automatically generated: i) FPISTREAM: all of AVX load instructions are replaced by PXOR instructions (to avoid introduction of extra dependencies) and all of the AVX store instructions were replaced by NOP instructions (issuing no operation but preserving the binary size). FPISTREAM corresponds to the ideal case where all of the data access are suppressed. ii) MIS-TREAM: all of the AVX arithmetic instructions were replaced by NOP instructions. By comparing cycle counts of FPISTREAM, MISTREAM binaries with the original binary, potential performance problems due to data access (essentially cache misses) Table 1. Single core performance (GFlops/s) of the two hot routines: inversion (DP), matrix products (SP), and of the entire single-core executable. Measurements were performed on an Intel Xeon E31240, 3.30GHz, with a theoretical peak performance 52.8 GFlops/s (SP) and 26.4 GFlops/s (DP). The values in parenthesis are the percentages with respect to the peak. The turbo feature was turned off, and the results were obtained using Likwid performance tools. [START_REF] Treibig | Likwid: A lightweight performance-oriented tool suite for x86 multicore environments[END_REF]. 1 As the matrix to invert is block-diagonal with two N e /2 × N e /2 blocks, the inversion runs on the two sub-matrices.

System sizes

Matrix inversion could be easily detected. Such an analysis revealed that for most of the hot loops, data access was accounting for less than 30% of the original time, indicating an excellent usage of the caches. Measurement of the binaries were performed directly with the whole application running allowing to take into account runtime context for the loops measured.

General Conclusion

In December 2011, GENCI gave us the opportunity to test our program on Curie (TGCC-GENCI, France) while the engineers were still installing the machine. At that time, up to 4 800 nodes (76 800 cores) were available to us for two sessions of 12 hours. As the engineers were still running a few benchmarks, our runs were divided into 3-hour jobs using 400 nodes (6 400 cores). In this way the engineers were able to acquire resources while our job was running. This aspect points out the importance of our flexible parallel model, but makes it impossible to evaluate rigorously the parallel efficiency. At some point, all the available nodes were running for our calculation during several hours. As we had previously measured a sustained performance of 200 GFlops/s per node for this run, we can safely extrapolate to a sustained value of ∼960 TFlops/s (mixed single/double-precision) corresponding to about 38 % of the peak performance of the whole machine for a few hours. As the machine was still in the test phase, we experienced a few hardware problems and maintenance shutdowns of some nodes during the runs. Quite interestingly, it turns out to be an opportunity for us to test the robustness of our program: it gave us the confirmation that our fault-tolerant scheme is indeed fully functional.

In this work we have presented a number of important improvements implemented in the QMC=Chem program that beautifully illustrate the extremely favorable computational aspects of the QMC algorithms. In view of the rapid evolution of computational infrastructures towards more and more numerous and efficient processors it is clear that such aspects could be essential in giving a definite advantage to QMC with respect to other approaches based on deterministic linear algebra-type algorithms.

Fig. 1 .

 1 Fig. 1. Graphical representation of a QMC simulation. Each process generates blocks, each block being composed of N walk walkers realizing N step Monte Carlo steps.

Fig. 2 .

 2 Fig. 2. The communication architecture of QMC=Chem.

Fig. 3 .

 3 Fig. 3. Number of computed blocks as a function of the number of CPU cores for a fixed computation time.

This is not representative of a real simulation since it is far too small: the time spent in network communications will be over-estimated compared to real simulations where the time to compute one block is typically much greater (10 minutes and more).

Acknowledgments. AS and MC would like to thank ANR for support through Grant No ANR 2011 BS08 004 01. This work was possible thanks to the generous computational support from CALMIP (Toulouse) under the allocation 2011-0510, GENCI, CCRT (CEA), and PRACE through a preparatory access (No PA0356).