N Melab

J Gmys

M Mezmaz

D Tuyttens
email: daniel.tuyttens@umons.ac.be

Many-core Branch-and-Bound for GPU accelerators and MIC coprocessors

Coprocessors are increasingly becoming key building blocks of High Performance Computing platforms. These many-core energy-efficient devices boost the performance of traditional processors. On the other hand, Branch-and-Bound (B&B) algorithms are tree-based exact methods for solving to optimality combinatorial optimization problems (COPs). Solving large COPs results in the generation of a very large pool of subproblems and the evaluation of their associated lower bounds. Generating and evaluating those subproblems on coprocessors raises several issues including processor-coprocessor data transfer optimization, vectorization, thread divergence, and so on. In this paper, we investigate the offload-based parallel design and implementation of B&B algorithms for coprocessors addressing these issues. Two major many-core architectures are considered and compared: Nvidia GPU and Intel MIC. The proposed approaches have been experimented using the Flow-Shop scheduling problem and two hardware configurations equivalent in terms of energy consumption: Nvidia Tesla K40 and Intel Xeon Phi 5110P. The reported results show that the GPU-accelerated approach outperforms the MIC offload-based one even in its vectorized version. Moreover, vectorization improves the efficiency of the MIC offload-based approach with a factor of two.

Introduction

Coprocessors are increasingly becoming key building blocks of High Performance Computing platforms. In addition to their energy efficiency, they boost the performance of traditional processors through the combination of a larger number of processing cores, vector-SIMD processing and multi-threading. Currently, the most used coprocessors (Top500 ranking [START_REF]Top500 HPC international ranking[END_REF] of June 2017) are Nvidia GPU accelerators and Intel MIC coprocessors. The former are composed of a large number of slim cores while the latter integrate a relatively smaller number of streamlined largish cores relying on SIMD processing. Today, coprocessors allow to achieve peak performance of the order of several TeraFlops. Nevertheless, it is often difficult for the programmers to extract a large portion of the theroretically available performance. Indeed, the specific features of these coprocessors raise several issues including the optimization of data transfer between the processor and its coprocessor, vectorization, data placement optimization, etc. More details on the these hardware features and related challenging issues are given in the next sections.

In this paper, the focus is put on the Branch-and-Bound (B&B) algorithm [START_REF] Lawler | Branch-and-bound methods: A survey[END_REF]. Recently, the parallelization of B&B has been revisited for multi-core (clusters of) processors [START_REF] Barreto | Parallel branch and bound algorithm-a comparison between serial, openmp and mpi implementations[END_REF] and GPU [START_REF] Melab | A gpu-accelerated branch-and-bound algorithm for the flow-shop scheduling problem[END_REF][START_REF] Lalami | Contribution à la résolution de problèms d'optimisation combinatoire: méthodes séquentielles et parallèles[END_REF][START_REF] Chakroun | Reducing thread divergence in a gpuaccelerated branch-and-bound algorithm[END_REF] and their combination [START_REF] Chakroun | Combining multi-core and gpu computing for solving combinatorial optimization problems[END_REF], [START_REF] Vu | Heterogeneity and locality-aware work stealing for large scale branch-and-bound irregular algorithms[END_REF]. In this paper, we investigate the parallelization of Branch-and-Bound (B&B) algorithms for coprocessors (GPU and MIC 1 ). B&B algorithms are well-known methods for solving to optimality NP-hard optimization problems 2 . They are based on an implicit enumeration of the solutions composing the search space associated with the problem to be tackled. The search space is explored by dynamically building a tree whose root node designates the original problem. Each internal or intermediate node represents a subproblem obtained by the decomposition of the subproblem associated with its parent. The leaf nodes designate potential solutions or subproblems that cannot be decomposed. The construction of the B&B tree and its exploration are performed using four operators: branching, bounding, selection and elimination. The algorithm proceeds in several iterations to progressively improve the best solution found so far. The generated and not yet examined subproblems are kept in a pool, that is initialized with the original problem. At each iteration, the selection operator is used to select a subproblem from this pool, according to some strategy (depth-first, best-first,...). The branching operator performs its decomposition into other subproblems. The bounding operator calculates a lower bound of each generated subproblem. Each subproblem having a lower bound higher than the upper bound is discarded by the elimination operator. It means that this subproblem will not be decomposed.

The bounding operator is the most time consuming part of a B&B algorithm. Indeed, in [START_REF] Chakroun | Reducing thread divergence in a gpuaccelerated branch-and-bound algorithm[END_REF], it is shown that this operator represents on average between 98 % and 99 % of B&B applied to the Flow-Shop problem. Such result demonstrates that the bounding operator needs massively parallel computing. The GPU-accelerated bounding has been investigated in [START_REF] Melab | Graphics processing unit-accelerated bounding for branch-and-bound applied to a permutation problem using data access optimization[END_REF]. The reported results show that the CPU-GPU data transfer is costly. Therefore, it is recommended to perform also the branching operator on GPU in order to generate locally the subproblems and evaluate their lower bounds. However, this raises some issues related to the highly irregular (in shape and size) nature of B&B tree: thread divergence, thread mapping, etc. Regarding the parallelization of B&B on MIC coprocessors, it has been addressed in our work proposed in [START_REF] Melab | Parallel branch-and-bound using private ivmbased work stealing on xeon phi MIC coprocessor[END_REF] using the native mode. In this mode, the coprocessor is standalone and executes the whole B&B algorithm. However, to the best of our knowledge, the parallelization of B&B on MIC coprocessors has not been addressed using the offload (GPU-like) mode.

The objective of this paper is to investigate the parallel bounding model combined with the parallel tree exploration model of B&B algorithms to allow highly efficient solving of large instances of the Flow-Shop problem on GPU accelerators and Intel MIC Xeon Phi coprocessors. In Section 2, we first present the general design of the coprocessor-accelerated B&B. In Section 3 (respectively Section 4), we describe the implementation of the GPU-accelerated (respectively Phiaccelerated) approach. In Section 5, we report some experimental results comparing the two coprocessor-based many-core implementations. Finally, some conclusions are drawn in Section 6.

Coprocessor-accelerated B&B: the general design

In this section, the general coprocessor-accelerated parallel approach is first presented. Then, the Flow-Shop permutation scheduling problem is presented as it is used as a use case to validate the approach.

The parallel model

As mentioned above, the coprocessors are many-core devices dedicated to massively parallel computing. Therefore, to take maximum advantage of the computational power provided by these coprocessors these latter should be fed by a large number of computations. In our proposed parallel coprocessor-based approach, as illustrated in Figure 1, several B&B trees are explored in order to generate multiple pools maximizing the use of the coprocessor cores. During the exploration of each pool, except the selection and pruning operators which are performed on the processor the branching and bounding operators are executed on the coprocessor. As shown in Figure 1, on the processor side, at each iteration of the exploration process a set of tree nodes (whose size is a user-parameter) is selected. The selected set of nodes is offloaded to the coprocessor to be processed. On the coprocessor side, as illustrated in Algorithm 10, each parent node is processed by a thread, which starts by getting its identifier thdId (line 2). Using that identifier and a mapping strategy, the thread determines the parent node to branch generating a subproblem child (line 3). If the generated subproblem is a leaf (solution) this latter is evaluated (line 5), otherwise its associated lower bound (LB) is computed (line 7) and inserted in a global pool of nodes with their associated bounds (line 9). This pool is returned back to the CPU host. Every child having a lower bound greater than the cost of the best solution found so far is pruned on CPU. All the non-pruned children are inserted into the pools. The process is iterated until the exploration is completed and the optimal solution is found.

GPU-based implementation of B&B

In this section, we first present the parallelization model on GPU. To do that we recall the hardware view of GPU, its parallel programming model and its associated algorithmic challenging issues. Then, we show how these issues are dealt with in the implementation of the GPU-accelerated B&B.

The Flow-Shop scheduling problem (FSP)

The Flow-Shop scheduling problem is a permutation problem which consists in scheduling n jobs on m machines [START_REF] Garey | The complexity of flow-shop and job-shop scheduling[END_REF]. For instance, in Figure 2, 4 jobs designated by different colors are scheduled on 3 machines.

The scheduling must be done with respect to two constraints: each machine cannot be simultaneously assigned to more than one job (to more than one color), and the execution order of the jobs (the colors) is the same on all the machines. The objective is to minimize the makespan, i.e. the termination date of the last job on the last machine. In this example, the solution consists in scheduling first the red job, then the blue job and orange one, and finally the green job. The solution can be coded as a permutation (3, 4, 2, 1) and its cost is the termination date of the green job on the machine M3. In this work, the lower bound function proposed by Lageweg et al. [START_REF] Lageweg | A General bounding scheme for the permutation flowshop problem[END_REF] is used in our bounding operator (Algorithm 2). This lower bound is mainly based on Johnson's theorem [START_REF] Johnson | Optimal two and three-stage production schedules with setup times included[END_REF] which provides a procedure for finding an optimal solution for Flow-Shop scheduling problem with 2 machines. This bound is known for its good results and has a computational complexity of O(m 2 nlog(n)), where n is the number of jobs and m the number of machines. For large values of the parameters m and n, the problem is time-intensive. More details on the lower bound and its computational complexity can be found in [START_REF] Melab | Graphics processing unit-accelerated bounding for branch-and-bound applied to a permutation problem using data access optimization[END_REF].

Parallelization on GPU

For a long time, GPU computing has been used to speed up image and video processing. Since 2006, with the introduction by Nvidia of its Cuda software toolkit the use of GPUs has been extended to numerous other application domains including combinatorial optimization. The popularity of Cuda is due to its simplicity as it is an extension of the C language with data parallel features. The principle is easy: the programmer writes a code for one thread (kernel) and can instantiate it on a large number of threads to allow massive parallel computing on GPU. In addition, Cuda is portable between successive generations allowing transparent scalability of Cuda applications.

Before the Cuda parallel model is presented, let us recall the hardware architecture of a GPU device. As shown in Figure 3, a GPU is a coprocessor, coupled to a CPU through a PCI Express bus. In the Cuda vocabulary, the processor is called "host" and the GPU is called "device". The GPU is composed by a set of streaming multi-processors (processors) including each a pool of 32-bit or 64-bit SIMD processors (processing cores). For instance, a Kepler GPU device contains 15 processors of 192 Cuda cores for a total of 2880 Cuda cores [START_REF] Nvidia | NVIDIA's Next Generation CUDA Compute Architecture: Kepler[END_REF]. A GPU is also composed of several memories including global and local off-chip memories, and a shared memory, registers and a cache memory. These memories have different characteristics in terms of size and access latency. For instance, the global memory is big and has a long latency while registers are small and fast memories. From software programming point of view, as illustrated in Figure 4 a GPU Cuda-based parallel program is composed of two parts: a "host" part and a "device" part. The host part is a serial or weakly parallel code because the number of CPU cores is small compared to the number of GPU cores. The device part is massively parallel because a GPU contains from hundreds to thousands of processing cores. During the execution of a parallel program the host part offloads streams of threads to the GPU device to be executed according to a two-level parallelism: at the higher level the processors (or SMX) execute the thread kernel according to the Sin-gle Program Multiple Data (or SPMD) model. At the lower level (intra-SMX), the threads are executed according to the Single Instruction Multiple Data (or SIMD) or Single Instruction Multiple Thread (or SIMT) model. Indeed, inside each processor the instruction flow composing a thread kernel is executed according to the SIMD model. The threads offloaded from CPU host to GPU device are organized by the programmer in a hierarchical way into grids of blocks of threads. Grids are arrays 1D or 2D of blocks and blocks are arrays 1D, 2D or 3D of threads. The thread organization corresponds to the organization of application data which are often vectors, matrices or volumes. As shown in Figure 5, the blocks are assigned to the SMXs by the Cuda runtime. Inside each SMX each block is split into warps i.e. pools of 32 threads. Warps are scheduling units i.e. the threads are executed by pools of 32. This allows to overlap the memory access latency by computation. Context switching is very fast as each thread has its own registers.

To sum up, from algorithmic and software programming point of view at least three issues should be addressed: (1) the optimization of the data transfer between CPU and GPU; (2) the optimization of the data placement on the hierarchy of memories of the GPU having different sizes and latencies; and (3) thread or branch divergence management especially for irregular applications.

Parallelization of B&B for GPU

The implementation on GPU of the coprocessor-based B&B according to the general design presented in Section 2 requires to address the challenging issues quoted above. First, to deal with the CPU-GPU data transfer optimization the branching operator, which generates tree nodes or subproblems, is moved to GPU. The execution of the branching operator on the GPU device allows one to avoid the transfer of the branched parent nodes from CPU to GPU which is costly. However, this raises other issues related to thread granularity and mapping. Indeed, if a parent node is processed entirely (branching and bounding) by a single thread there will be a load imbalance leading to thread/branch divergence. In fact, the parent nodes may have different numbers of children as they are located at different levels in the B&B tree. To deal with this problem the processing of each thread is limited to a single node, meaning that each thread generates and evaluates only one child of the parent node.

Second, to tackle the problem related to data placement optimization on GPU we have followed the recommendation proposed in [START_REF] Melab | Graphics processing unit-accelerated bounding for branch-and-bound applied to a permutation problem using data access optimization[END_REF]. Indeed, as illustrated in Algorithm 2, the implementation of the lower bound algorithm includes 6 data structures: the matrix PT M of the processing times of the jobs, the matrix of lags LM, the Johnson's matrix JM, the matrix RM of the earliest starting times of jobs, the matrix QM of their lowest latency times and the matrix MM containing the couples of machines. The algorithm needs as input a subproblem defined as a permutation with some jobs already scheduled at its beginning and/or its end.

The semantics of these data structures is not the focus of this paper. For more details on these ones and on the lower bound please refer to [START_REF] Melab | Graphics processing unit-accelerated bounding for branch-and-bound applied to a permutation problem using data access optimization[END_REF]. The focus is rather put here on the optimization of the placement of these data structures on the different memories of the GPU device. Due to the limited size of the shared memory, the matrices do not fit in all together, especially for large problem instances. Based on the complexity analysis of these data structures and an experimental study conducted in [START_REF] Melab | Graphics processing unit-accelerated bounding for branch-and-bound applied to a permutation problem using data access optimization[END_REF], it is suggested to put in the shared memory the Johnson's and/or processing time matrices (JM and PT M). The other data structures are mapped either to the global memory or to the constant memory. Such data placement allows one to achieve accelerations of more than ×100 compared to a single-CPU core for (k = 0 → N(N-1)

2

) do 6: tmp0, tmp1, ma0, ma1 ← InitFun(k, nbFixed, MM, RM) 7:

for return LB 18: end function serial execution of B&B. In our implementation, MM, JM, PT M and LM are put on the constant memory. A part of PT M is then moved to the shared memory. The other matrices are stored on the global memory.

(j = 0 → n) do ∼ 70%
Third, the parallelization of irregular applications such as B&B applied to permutation problems due to the thread or branch divergence issue [START_REF] Chakroun | Reducing thread divergence in a gpuaccelerated branch-and-bound algorithm[END_REF]. In our implementation, the irregularity is due to two factors: as the tree nodes have different levels they require different amounts of work (number of children) ; on the other hand, the lower bound function, as it can be seen in Algorithm 2, includes several conditional instructions and loops. To deal with the first factor, each thread handles only one child as mentioned above. Therefore, all the threads perform the same amount of work. To tackle the second problem, we have reused the refactoring approach proposed in [START_REF] Chakroun | Reducing thread divergence in a gpuaccelerated branch-and-bound algorithm[END_REF] even if the achieved performance improvement is not significant as the size of the factorized branches is small.

Finally, the mechanisms used on GPU are not experimented individually here because their efficiency has been demonstrated in [START_REF] Chakroun | Reducing thread divergence in a gpuaccelerated branch-and-bound algorithm[END_REF][START_REF] Melab | Graphics processing unit-accelerated bounding for branch-and-bound applied to a permutation problem using data access optimization[END_REF][START_REF] Chakroun | Parallel heterogeneous branch and bound algorithms for multi-core and multigpu environments[END_REF]]. These three citations can be used for further details on the mechanisms. However, the performance of the whole GPU-accelerated B&B is evaluated and compared to the performance of the Xeon Phi-based B&B in Section 5.

MIC-based implementation of B&B

In this section, we first present the parallelization model on MIC architecture. To do that we recall the hardware view of Intel Xeon Phi, its parallel programming model and its associated algorithmic challenging issues. Then, we show how these issues are dealt with in the implementation of the Intel Xeon Phi-accelerated B&B.

Parallelization on Intel Xeon Phi

The market of accelerators has been dominated by Nvidia during several years. Since recently, they are faced to the competition of Intel with its Many Integrated Cores (MIC)-based Xeon Phi. This latter is a coprocessor coupled to the processor through a PCI Express bus [START_REF] Chrysos | Intel Xeon Phi X100 Family Coprocessor -the Architecture[END_REF]. As illustrated in Figure 6, a typical platform consists of one to two Intel Xeon processor(s) (CPUs) and one to eight (two in this figure) Intel Xeon Phi coprocessors per host. Several of such platforms may be interconnected to form a cluster or supercomputer. From the hardware point of view, the Xeon Phi board has one Knights Corner (KNC) processor, the first production chip based on the MIC architecture, and 8 GB of GDDR5 RAM. As illustrated in Figure 7, KNC integrates up to 60 CPU-cores interconnected by a high-speed bi-directional ring, and runs at over 1 GHz. It connects to its private external memory with a peak bandwidth of over 320 Gbps. The cores are based on the Intel Pentium architecture. Each core has 32 KB of L1 data and instruction cache, 512 KB of L2 data cache, and a 512-bit vector Floating Point Unit (FPU). This latter performs fused-multiply-add (FMA) operations. Therefore, the peak performance is about 32 (resp. 16) GFlops in single (resp. double) precision.

Consequently, the KNC delivers a peak performance of about 2 (resp. 1) TFlops in single (resp. double) precision. From a programming standpoint, the key is to treat the Intel Xeon Phi coprocessor as an x86-based SMP-on-a-chip with over 50 cores, with multiple threads per core and 512-bit SIMD instructions. From programming language point of view, Intel Xeon Phi is more accessible than Nvidia GPU because it can be programmed using standard programming environments such as OpenMP, MPI, Cilk Plus and Posix Threads. However, to achieve higher performance one should consider two fundamental features: scaling through locality and Simultaneous Multi-Threading (SMT) and vectorization. On the other hand, as an Intel Xeon Phi coprocessor runs an operating system (Linux) and has its own IP address, there are two ways to involve it in a parallel program:

• A processor-centric "offload" mode where the program is viewed, like for GPUs, as running on processors and offloading work to coprocessors. Another issue which arises with the offload mode is, like for GPU, the optimization of the data transfer between the processor and the coprocessor. • A "native" mode where the program runs natively on only coprocessors or on coprocessors and processors together. In the latter case, the two devices may communicate with each other by various methods. At least two other issues arise in this case: the optimization of the data/task partitioning and the communication between the processor and the coprocessor. The challenge for data/task partitioning is the load balancing between the "big cores" of the processor and the "little" cores of the coprocessor. Regarding the communication, the challenge is to overlap the communications by the computation, to manage data locality, etc.

Offload-based parallelization of B&B for Intel Xeon Phi

In this paper, we investigate the offload mode for the parallelization of B&B algorithms on Intel Xeon Phi coprocessors. In the proposed approach, the processorcoprocessor data transfer optimization approach is the same as for GPU. Indeed, the branching operator is performed on Phi, this allows one to reduce the cost of data transfer between the two devices. On the other hand, all the data structures which are not modified between offloading operations are offloaded once.

One of the major mechanisms allowing performance improvement on Intel Xeon Phi is vectorization. Different levels are provided ranging from compiler-based automatic easy-to-use vectorization to manual and programmer control vectorization. In our work, as quoted previously, the most consuming part of the B&B algorithm is computation of the lower bound function (Algorithm 2). In this paper, we propose a vectorization method of the lower bound function focusing on its most computeintensive portion and the main data-dependencies. This portion of code is the inner for-loop (line 7-13) which consumes about 70% of the bounding time. The body of this inner loop is executed J 2 ×(J-1) 2 times. Regarding data dependencies, the statement in line 11, including a dependency of current tmp1 on tmp1 from previous iteration prevents vectorization (icc does not auto-vectorize it). In addition, except for line 15, the iterations of the outer loop are independent (private variables: tmp0, tmp1, ind0, ind1, current). However, only the inner loop of a loop nest may be vectorized 3 .

In order to allow more vectorization than provided by the compiler, the order of the nested loops must be inverted as illustrated in the vectorized lower bound function (Algorithm 3). For auto-vectorization by the compiler it is preferable to write small separate loops, rather than merging into a single loop. The outer loop is thus split into 3 separate serial loops and a max-reduce operation (line 21) in order to isolate the k-dependent instructions from the inner-loop. The cost to pay for this is to declare the scalars tmp0, tmp1, ma0, ma1 as arrays (resp. T mp0, T mp1, Ma0, Ma1) of size J 2 ×(J-1)

2

. The same strategy on GPU severely breaks down performance due to the memory problem (these intermediate variables are no longer stored into registers). Even with the highest optimization level activated (-O3) the Intel compiler (icc) still needs the hint "#pragma ivdep to vectorize the inner loop (line 10) successfully. The two other for-loop are auto-vectorized.

Experimentation

In this section, we present an experimental study of the proposed many-core approaches using GPU and Intel Xeon Phi and compare them. We first present the for

(k = 0 → n(n-1) 2 ) do 6: Tmp0[k], Tmp1[k], Ma0[k], Ma1[k] ← InitFun(k, nbFixed, MM, RM) 7:
end for 8:

for (j = 0 → J) do permute loop-order 9:

#pragma ivdep 10:

for

(k = 0 → n(n-1) 2 ) do inner loop vectorizable 11: job←JM[j][k] transpose JM 12: if (SM[job]==0) then 13: Tmp0[k] += PTM[Ma0[k]][job] 14: Tmp1[k] ← max(Tmp1[k], Tmp0[k] + LM[k][job]) + PTM[Ma1[k]][job] 15:
end if 16:

end for 17:

end for 18:

for

(k = 0 → n(n-1) 2 ) do 19: Tmp1[k]←EndFun(Tmp0[k], Tmp1[k], k, nbFixed, QM) 20:
end for 21:

LB←max-reduce(Tmp1[]) 22:

return LB 23: end function hardware and software testbeds and some parameter setting used for our experiments. Then we report and discuss some experimental results.

Hardware and software testbed and parameter setting

As shown in Table 1, all the experiments are run on a computer from Grid'5000 [START_REF]Grid'5000: A large scale and highly reconfigurable experimental grid testbed[END_REF] which is a single socket Sandy Bridge machine (Intel Xeon E5-2650, 64 GB RAM), equipped with one MIC coprocessor Intel Xeon Phi 5110P (60 physical cores at 1.053 GHz, 320 GB/s of memory bandwidth). The operating system is a Centos 6 (OS officially supported by Intel for the Xeon Phi) and the compiler is Intel icc for Intel devices. The GPU device is an NVIDIA Tesla K40 based on the Kepler GK110B architecture. It includes 2880 Cuda cores at 0.745 GHz, 288 GB/s of memory bandwidth. For the GPU-accelerated implementation, the NVIDIA CUDA Toolkit release 5.0.35 is used together with the gcc version 4.4.7. For all experiments the compilation level 3 (-O3) is used. On the other hand, the UNIX time command is used to measure the elapsed execution time for each Flow-Shop instance. The GFLOPS (SP) row is obtained using the following computations: • MIC: 32( f lops cycle ) × 60(cores) × 1.053( GHz core ) = 2 022 GFlops. Finally, the last row (Thermal Design Power, TDP) indicates that the hardware configurations of the two compared coprocessors (Tesla K40 and Xeon Phi 5110P) are equivalent in terms of energy consumption. As quoted in Section 2, the coprocessors are many-core devices dedicated to massively parallel computing. Therefore, they need to be fed by a large number of computations. To do that, many B&B trees are explored in order to generate multiple pools maximizing the use of the coprocessor cores. Before the experiments are performed, the number M of pools to be created is calibrated through a series of experiments on the problem instance Ta028. The experimental results are reported in Figure 8, showing the execution time as an average over 10 independent runs. Error bars show the corresponding standard deviation. Based on the figure, the number of pools is fixed to 800 (resp. 2200) for the GPU-accelerated (resp. Xeon Phi-based) approach.

Instance Ta028 (8.1x10 6

Experimental results

The performance analysis of the contributions of this paper we have compared the two many-core approaches: the GPU-accelerated approach and the MIC-accelerated approach. In order to get a fair comparison between them the offload mode is considered for the MIC-based approach. In our experiments, simultaneous multi-threading (SMT with 2 threads per core) is considered for the MIC-based approach. The two approaches have been experimented using the 10 instances (Ta021 -Ta030) of the Taillard's problem 20 jobs on 20 machines. The best solution found so far is initialized to the optimal solution to guarantee that the amount of work (explored nodes) is the same for each of the experimented approaches. Therefore, the objective of the resolution is to prove the optimality of the initial solution. Obviously, to provide the optimal solution the same algorithm is used. One has just to initialize the best solution found so far to infinity for a minimization problem. The obtained experimental results are reported in The first two columns of the table contain respectively the numbers of the 10 solved problem instances and their associated search space sizes in millions of nodes. The following columns designate the exploration time in seconds obtained using respectively the GPU-accelerated approach, and the vectorized and non-vectorized MIC-based approaches using 236 threads.

In [START_REF] Saini | Early multinode performance evaluation of a knights corner (KNC) based NASA supercomputer[END_REF], the authors report some experimental results demonstrating the impact on data transfer overhead of the Coprocessor Offload Infrastructure (COI) daemon in the offload mode. The COI daemon runs the services required to support data transfer for offload on a dedicated core. The reported results show that it is beneficial to avoid using this core for user code, i.e., one should use only 59 cores. Following this recommendation, we have used 236 threads (corresponding to 59 cores) in our experiments for the offload-based MIC approach. Indeed, according to our experimental results, using all 60 cores (240 threads) for computations incurs performance penalties up to 80%.

From the last row of the table, two major observations can be made. First, the GPU-accelerated approach clearly outperforms the MIC-based approach even in its vectorized version. Second, vectorization allows one to speed up the MIC-based approach with a factor of two.

Conclusion

According to the recent Top500 ranking (July 2015), it is confirmed that hybrid HPC platforms including coprocessors is the trend towards the exascale era. On the other hand, it appears that the market of hybrid HPC is dominated by Nvidia followed by Intel with its Xeon Phi. In this paper, we have revisited the parallelization of B&B algorithms for many-core coprocessors, in particular Nvidia GPU and Intel Xeon Phi. From the design point of view, we have combined two hierarchical parallel models: the parallel tree exploration model and the parallel bounding. The bounding operator is performed on the coprocessor because, on the one hand, it is the most time-consuming part of the B&B algorithm. On the other hand, it is massively data parallel and thus well-suited for coprocessors. In addition, the branching operator, which generates tree nodes during the exploration process, is also performed on the coprocessor to minimize the cost of their offloading from the processor to the coprocessor.

Such coprocessor-based design of B&B algorithms gives rise to other issues: thread mapping, thread/branch divergence and data placement optimization on GPU, and vectorization on Intel Xeon Phi. From the GPU side, we have reused some recommendations proposed in [START_REF] Chakroun | Parallel heterogeneous branch and bound algorithms for multi-core and multigpu environments[END_REF]. For the MIC-based approach, we have proposed a vectorization method for the lower bound function. The many-core implementations have been experimented on the 10 instances of the 20 jobs-on-20 machines problem using equivalent hardware configurations in terms of energy consumption. The reported results show that the GPU-accelerated approach outperforms the MIC offload-based one even in its vectorized version. Moreover, vectorization improves the efficiency of the MIC offload-based approach with a factor of two.

In the near future, we plan to extend this work to a cluster-level parallelization combining multi-core GPU-accelerated and MIC-based computing. From an application point of view, the objective is to solve to optimality challenging difficult and unsolved Flow-Shop instances as we did for one 50 × 20 problem instance with grid computing [START_REF] Mezmaz | A grid-enabled branch and bound algorithm for solving challenging combinatorial optimization problems[END_REF]. Finally, we plan to investigate other lower bound functions to deal with other combinatorial optimization problems.

Fig. 1 :

 1 Fig. 1: Coprocessor-based B&B

Fig. 2 :

 2 Fig. 2: Illustration of the Flow-Shop problem.

Fig. 3 :

 3 Fig. 3: Hardware view: GPU = coprocessor of CPU.

Fig. 4 :

 4 Fig. 4: Software view: Parallel program = weakly parallel/serial host code + massively parallel device code.

Fig. 5 :

 5 Fig. 5: Software view: Parallel program = grid(s) of block(s) of threads executed as warps of 32 threads.

Algorithm 2

 2 Computation of lower bound (un-vectorized) input: subproblem = {permutation, nbFixed (#jobs fixed)} constant data (MM, JM, PTM, LM) output: lower bound (LB) of subproblem 1: n := #jobs 2: function COMPUTEBOUND 3: RM, QM, SM ← InitTabs(permutation, nbFixed)

Fig. 6 :

 6 Fig. 6: Hardware view: Intel Xeon Phi = coprocessor of CPU.

Fig. 7 :

 7 Fig. 7: Hardware view: Knights Corner core.

Algorithm 3

 3 Computation of lower bound (vectorized) input: subproblem = {permutation, nbFixed (#jobs fixed)}, constant data (MM, JM, PTM, LM) output: lower bound (LB) of subproblem 1: n := #jobs 2: function COMPUTE LB VECTORIZED 3: RM, QM, SM ← InitTabs(permutation, nbFixed)

•

  K40: 2( f lops cycle ) × 15(SM) × 192( cores SM ) × 0.745( GHz core ) = 4 291 GFlops.

Fig. 8 :

 8 Fig. 8: Elapsed time vs. Number of pools.

Algorithm 1

 1 Kernel of the computation of the lower bounds on the coprocessor.

	input: parent-subproblems, mapping
	output: lower bounds of children subproblems
	1: kernel EVALUATEONCOPROCESSOR
	2:	thdId←getThreadId()
	3:	child-subproblem←generateChild(thdId, mapping, parent-subproblems)
	4:	if isLeaf(child-subproblem) then
	5:	LB←evaluateSolution(child-subproblem)
	6:	else
	7:	LB←computeBound(child-subproblem)
	8:	end if
	9:	poolOfBounds[thdId]←LB
	10: end kernel

Table 1 :

 1 Hardware execution platform

		NVIDIA Tesla K40 Intel Xeon-Phi 5110P
	#Physical cores	15	60
	#Logical cores	2 880	240
	clock(Ghz)	0.745	1.053
	GFLOPS(SP)	4 291	2 022
	GFLOPS(DP)	1 430	1 011
	SIMD	N/A	512-bit
	L2 Cache(MB)	1.5	30
	Mem BW(GB/s)	288	320
	TDP (Watt)	235	225
	Launch date	Oct'13	Oct'12

Table 2 .

 2

		CPU-seq GPU	Xeon-Phi
	Inst. nodes (×10 6 )			Vect No-Vect
	21 41.4	21 752	174 571	1 220
	22 22.1	10 425	90 287	617
	23 140.8	70 386	585 1 856 4 116
	24 40.1	17 635	153 457	1 021
	25 41.4	22 446	177 599	1 203
	26 71.4	30 706	273 829	1 787
	27 57.1	25 073	211 597	1 489
	28 8.1	4 043	34 109	239
	29 6.8	3 262	28	87	193
	30 1.6	795	7	24	50
	Avg 43.1	20 653	173 542	1 193

Table 2 :

 2 Exploration time (in seconds) for solving Flow-Shop instances Ta021-Ta030 initialized at optimal solution

GPU and MIC stand for respectively Graphics Processing Unit and Many Integrated Cores.

An optimization problem consists of minimizing or maximizing a cost function. Without loss of generality, in this paper the minimization case and the permutation Flow-Shop scheduling problem are considered.

http://d3f8ykwhia686p.cloudfront.net/1live/intel/CompilerAutovectorizationGuide.eps

Acknowledgement

Experiments presented in this paper were carried out using the Intel Xeon Phi of Digitalis (http://digitalis.imag.fr) and Grid'5000 testbed (https://www.grid5000.fr). Therefore, we would like to thank the Digitalis staff and especially Pierre Neyron from LIG/CNRS Lab for making the MIC processor fully functional and available. Grid'5000 is supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations.


Céline Merlet

Mathieu Salanne

Benjamin Rotenberg

Paul Madden

Imidazolium Ionic Liquid Interfaces with Vapor and Graphite: Interfacial Tension and Capacitance from Coarse-Grained Molecular Simulations

A recently developed coarse-grained model (J. Phys. Chem. B, 2010, 114, 12629-12631), previously validated against experimental data for a number of bulk properties, is used in molecular dynamics simulations of two different interfaces involving the ionic liquid [BMI][PF 6 ]. First, simulations of the liquid-vapor interface demonstrate that the model is able to predict the surface tension of the fluid (for which we obtain a value of 39.4 mN.m -1 at 400 K). Second, simulations were performed at constant potential differences applied between two graphite electrodes. From simulations with different applied potentials, the differential capacitances of the positive and negative electrodes can be calculated. It appears that both capacitances (C + = 3.9 µF.cm -2 for the positive electrode and C -= 4.8 µF.cm -2 for the negative electrode) agree very well with simulations results obtained with an all-atom model. The coarse-grained model also reproduces accurately the two-dimensional structure observed at the graphite-ionic liquid interface, namely a defective hexagonal lattice with a lattice spacing of approximately 10 Å.

I. INTRODUCTION

Ionic liquids are currently under intense examination as potential electrolytes for future capacitors [1,2,3] as they present a number of highly appropriate properties such as low vapor pressure, low combustibility, high thermal stability, good ionic conductivity and wide electrochemical windows. Moreover, their properties are in principle tunable thanks to the wide variety of cations and anions that can be mixed. Nevertheless, synthesizing, purifying and characterizing new ionic liquids remains expensive and systematic procedures to evaluate new materials are out of reach. Furthermore, a number of experimental observations are not understood, mainly because they are the result of microscopic phenomena which cannot be probed experimentally. Thus, simulation appears as an essential tool to explain what is observed in experiments and to design new materials with improved performances.

One particular experimental result that still suffers from a lack of explanation is the anomalously high value of the capacitance at nanoporous carbon electrodes [2,4,5,6,7]. Understanding this increase could lead to the design of new materials for future capacitors with improved performances. Furthermore, there is an issue about the shape of the capacitance curve when the external potential is varied. Some experimental [8,9] and theoretical [10,11] results suggest camel-shaped and bell-shaped variations for the capacitance but others tend to show that the capacitance is constant [12]. The capacitance curve can depend on several processes with different time constants as seen at the interface between Au(111) and 1-ethyl-3-methylimidazolium tris(pentafluoroethyl)-trifluorophosphate [13]. Thus, a representation of the processes happening at the atomic scale in these systems is needed.

Nowadays, ionic liquids are commonly studied by molecular dynamics and Monte Carlo simulations, and most of the structural, thermodynamic and dynamical bulk properties of these systems are well reproduced by recent force fields [14,15,16]. Interfaces involving ionic liquids are also well represented in the literature and molecular dynamics have already proven efficient in giving insights on liquid-solid interfaces and the capacitive behaviors of systems with simple electrodes [17,18,19]. Most of the liquid-electrode interfaces are studied with the electrode potential generated by holding the electrode at a constant charge [20] which is assigned equal values on equivalent sites. This may well be a reasonable approximation for planar electrodes but will become problematic for more complex geometries for which it is necessary to obtain charges which are consistent with the applied potential and the shape and dielectric properties of the electrode and the charge distribution of the ions across the interface [21].

To explain the capacitance behavior for nanoporous carbons, it is necessary to find a model sufficiently accurate to reproduce structural, dynamical and capacitive properties but simple enough to allow for calculations on complex systems with nanoporous electrodes. A recently developed four-site model for the ionic liquid [BMI][PF 6 ] [22] is examined in this work as a candidate for use in molecular dynamics simulations in complex systems. To deal with the problem of developing a consistent description of the charge on the roughly-shaped electrodes, we examine the use of a method which allows the electrodes to be treated as polarizable and maintained at a constant electrical potential, with the charges being determined self-consistently [23].

The remainder of this article is organized as follows. In section 2 we introduce the model chosen to describe the ionic liquid and the carbon electrodes. In section 3 we show that the model is able to predict liquid-vapor interfacial properties accurately. In sections 4 and 5 we demonstrate that the model is also efficient for the prediction of capacitive and structural properties of the interface between planar graphite electrodes and the ionic liquid.

II. MODEL AND SIMULATION METHODS

The ionic liquid is represented by a coarse-grained model developed by Roy and Maroncelli which accurately predicts a number of bulk properties [22] such as molar volume, isothermal compressibility, viscosity and diffusion. This model is actually an improvement of a previous one published by the same authors which showed poorer dynamical properties [24]. The most important difference between the two models is the use of reduced charges, i.e. all the charges of the first coarse-grained model were scaled by a same factor. This approach was also used successfully in other systems [14,25,26,27] and it is thought that reducing charges compensates for the missing polarization term in an average way in these force fields.

Going beyond the topic of ionic liquids, the Molecular Dynamics in Electronic Continuum model provides a theoretical framework for systematic accounting of the effects of electronic polarization through the use of rescaled charges [28]. A recent publication describes a methodology to calculate the scaling factor from DFT calculations of the electrostatic potential in liquid phase configurations [29]. Van

σ ij r ij ) 12 -( σ ij r ij ) 6 ] + q i q j
4πε 0 r ij where r ij is the distance between sites, ε 0 is the permittivity of free space and crossed parameters are calculated by Lorentz-Berthelot mixing rules. The parameters for the carbon atoms of the graphite electrodes are σ C = 3.37 Å and ε C = 0.23 kJ.mol -1 [30].

described using pairwise additive Lennard-Jones interactions. The parameters of the coarsegrained model developed by Roy and Maroncelli [22] to represent [BMI][PF 6 ] are recalled in table I and a schematic representation of the ionic liquid is given in figure 1. For the graphite electrodes, potential parameters for carbon atoms were set equal to σ C = 3.37 Å and ε C = 0.23 kJ.mol -1 following Cole and Klein [30]. All the cross parameters were obtained using the Lorentz-Berthelot mixing rules.

The electrodes are modelled following a method developed by Reed et al. [23], following a proposal by Siepmann and Sprik [31], which consists in determining the charge on each electrode atom at each MD step by requiring that the potential on this atom is constant and equal to a specified value. This condition therefore means that the electrode polarizes in the same way as a metal, where the potential inside the metal is a constant. The results obtained with this construction will be contrasted with results obtained by assigning constant charges to the electrode atoms. In the constant potential method, the electrode atom charges, represented by gaussian distributions centered on the carbon atoms, are obtained by minimizing the expression:

U = i q i (t) Ψ i ({q j (t)}) 2 + q i (t)κ √ 2π -Ψ 0 , (1)
where q i is the charge on electrode atom i, κ is the width of the gaussian distributions (κ = 0.5055 Å in this work), Ψ i ({q j }) is the potential at position i due to all the ions in the electrolyte and the other electrode atoms, the second term originates from the interaction of the gaussian distribution with itself and Ψ 0 is the externally applied potential. Instead of minimizing this expression using a conjugate gradient method only, we use a two-step procedure where the first step is a prediction of the charges by a polynomial expansion, using the previous steps [32,33], and the second step is a conjugate gradient minimization.

This allows a speed-up by a factor of 2 compared to conjugate-gradient minimization only.

This representation of the electrodes has been successfully applied to the interface between the molten salt LiCl and aluminium electrodes [34,35].

All the simulations were conducted with a time step of 2 fs and by applying a Nosé-Hoover thermostat with a weak relaxation time constant of 10 ps to maintain a temperature of 400 K.

All the simulation cells are orthorhombic with cell lengths such that L x ∼ L y < L z The electrostatic interactions are treated according to the periodic boundary conditions adopted (3D for liquid-vapor simulations, 2D for liquid-graphite simulations). The 2D Ewald summation is done as described in reference [23] and following recommendations of reference [36].

III. LIQUID-VAPOR INTERFACE: SURFACE TENSION

The coarse-grained model developed by Roy and Maroncelli has proven to be reliable for a number of bulk properties but has not yet been tested for interfacial properties. The simplest interface is the one between the ionic liquid and vacuum. The surface tension is a good test of the accuracy of the model for liquid-vapor interfaces as this property has been well studied by experiments [START_REF] Ghatee | Fluid Phase Equilibr[END_REF][START_REF] Freire | [END_REF] and simulations [14,39]. From simulations of an interface between [BMI][PF 6 ] and vacuum, it is possible to calculate the surface tension of the ionic liquid using the pressure tensor components [40,[START_REF] Frenkel | Understanding Molecular Simulation[END_REF]:

γ = Lz 2 < Π zz > - < Π xx + Π yy > 2 .
(2

)
The system consists in a slab (parallel to the xy plane) of ionic liquid sandwiched between two vacuum regions simulated with 3D periodic boundary conditions for 1 ns. The simulations were conducted with an orthorhombic cell with lengths such that L x ∼ L y = L z . The system contains 320 ion pairs for the first simulations when varying L z and 1280 ion pairs for the last simulation when L x and L y are changed.

In a first step, we analyze the effect of 3D periodic boundary conditions on the surface tension. In particular, there is a need to verify that the periodic slabs of ionic liquid do not interact with each other. In practice, when calculating the surface tension with different lengths in the z direction, it appears that there is a minimum length from which interactions can be considered negligible. Simulations were done with three different lengths L z , i.e. increasing the vacuum region size and keeping L x and L y constant. The resulting surface tensions are 14.6 mN.m -1 , 26.0 mN.m -1 and 26.0 mN.m -1 for L z equal to 300 Å, 400 Å and 500 Å respectively. It can thus be concluded that 500 Å is a sufficient length to avoid interactions between the slabs in the z direction.

In a second step, we analyze the finite-size effect in our system. Calculations of interfacial properties of fluids are subject to potential errors due to finite-size effects. The capillary wave theory suggests that the surface tension converges to the macroscopic value following [START_REF] Gelfand | [END_REF]:

γ(L x ) = γ ∞ + C L 2 x , (3)
where C is a constant and γ ∞ is the macroscopic observable (when L x = L y ). This prediction, first established for van der Waals fluids, is also applicable to simple molten salts as demonstrated for KI [43]. Recent studies show that for some liquids, the surface tension behaves somewhat differently [44]:

γ(L x ) = A 0 + A 1 exp(-α c L x ) L x cos(A 2 L x + A 3 ), (4)
where α c is a decay length and A i are constants. In both cases, when the interface area is larger than a threshold value, finite-size effects can be neglected. According to González-Melchor et al. [44], the minimum length one has to consider is approximately equal to 10σ ± with σ ± being the diameter of cations and anions. For the model used in this work, the diameter of the ions is around 5 Å which means that the lengths, L x = 32.2 Å and L y = 34.4 Å which were used in the previous simulations are too small to avoid finitesize dependency. The surface tension was thus evaluated for a system with a four times larger interface area (72.0 Å × 77.0 Å) and with a 500 Å length in the z direction. We then obtained a value of 39.4 mN.m -1 which is in very good agreement with experimental data (40.8 mN.m -1 at 393 K [START_REF] Ghatee | Fluid Phase Equilibr[END_REF], 41.07 mN.m -1 at 343.15 K [START_REF] Freire | [END_REF]). The coarse-grained model performs even better than previously described all-atom models which overestimate (47 mN.m -1 [14], 74 mN.m -1 [45]), or underestimate (37.3 mN.m -1 [39]) the surface tension of [BMI][PF 6 ] at room temperature.

IV. DIFFERENTIAL CAPACITANCE

The [BMI][PF 6 ]-graphite interface has now been well characterized in various studies.

It therefore provides us a good test of the ability of the coarse-grained model to predict capacitive and interfacial properties before tackling more complex interfaces, such as the ones involving nanoporous electrodes. Simulations were conducted with BMI + and PF - We perform simulations with nine potential differences ranging from 0.0 V to 4.0 V. When a new potential difference is set between the electrodes, the system is allowed to equilibrate for 200 ps before collecting data during 1 ns. For all the simulations, the potential difference is applied so that the potential applied to the positive electrode Ψ + is equal to the opposite of the negative electrode Ψ -, i.e.:

Ψ + = -Ψ -= ∆Ψ 0 2 . (5)
The highest considered potential difference was chosen to match with the experimental

electrochemical stability window of [BMI][PF 6 ] [1].
The differential capacitance of the two liquid-graphite interfaces can then be determined from the mean surface charges on the electrodes and the electrostatic potential profile across the cell which is given by Poisson's equation:

Ψ(z) = Ψ q (z) = Ψ q (z 0 ) - 1 ε 0 z z 0 dz z -∞ dz ρ q (z ), (6)
where z 0 is a reference point inside the left-hand electrode with Ψ q (z 0 ) = Ψ + and ρ q (z) is the charge density across the cell, which includes the contributions from the electrode atoms as well as from the charges on the sites of the ionic liquid molecules. The profiles obtained (see figure 3) are consistent with the ones obtained from other simulations of this type with molten salts [34]. The potential of the bulk Ψ bulk , is given by the plateau in the bulk liquid region, and the potential drop across each interface ∆Ψ ± = Ψ ± -Ψ bulk can be extracted from such profiles and used to characterize the interfaces. The differential capacitances of the positive and negative electrodes are obtained by differentiation of the surface charge with respect to the potential drop across the interface (see figure 4):

Ψ (V) ΔΨ 0 = 0.0 V ΔΨ 0 = 2.0 V ΔΨ 0 = 4.0 V
C diff = ∂σ ∂∆Ψ . (7)
It can be noted that the variation of the surface charge as a function of the potential drop across the interface shows two régimes which are clearly linear, indicating that the differential capacitance is constant there, as currently observed in experiments with supercapacitors [4,6]. Another quantity that can be determined from this plot is the potential of zero charge (PZC) which is zero for this system, in agreement with previous simulations 19,46]. Non-trivial values of the PZC were obtained in molten salt simulations of this type where the cations and anions differed in size [35]. This value, which is an absolute potential, cannot be compared directly with experimentally measured PZC as the experimental quantity is always given with respect to a reference electrode.

-3 -2 -1 0 1 2 3 !" +/-(V)
on [BMI][PF 6 ] [
The differential capacitances obtained to the positive (C + ) and negative (C -) sides of the PZC are summarized in table II. The difference between these values is linked to the different shapes of the ions and can also be noted in a plot of the ionic densities at the positive and negative electrodes (see figure 5). The differential capacitances calculated in

Source C + (µF.cm -2 ) C -(µF.cm -2 )
This work 3.9 4.8 this work are in excellent agreement with the ones obtained by Kislenko et al. [17] with an all-atom model and a constant charge of ± 8.2 µC.cm -2 but differ from the results of Feng et al. [19] which correspond to a constant charge of ± 11.2 µC.cm -2 . In particular, C -is greater than C + for this work but the results of reference [19] show the contrary. We underline here the fact that, for the present system, the highest potential difference, ∆Ψ 0 = 4.0 V, induces a mean surface charge of ± 9.4 µC.cm -2 . Because this potential difference corresponds to the electrochemical stability window of [BMI][PF 6 ] [1], it appears that a surface charge of ± 8.2 µC.cm -2 is realistic, while ± 11.2 µC.cm -2 might be too large to evaluate the differential capacitance of the system. To test this assumption, we performed a simulation in which a potential difference of 6.0 V was applied between the electrodes. This potential difference induces a surface charge of 11.9 µC.cm -2 on the graphite electrodes.

Points for this simulation are also reported in figure 4 demonstrating that they are outside the linear régimes. The high surface charge applied in the simulations of Feng et al. is thus the explanation of the contradictory results.

The good agreement between our results obtained with polarizable electrodes and previous results derived from constant charge simulations suggest that, for this system, constant charge simulations are sufficiently accurate to reproduce the capacitive behavior of the ionic liquid-planar graphite interface. To directly test, we performed a simulation with a constant charge of ± 4.5µC.cm -2 on each electrode, the charge being distributed uniformly on all the atoms of the graphite plane closest to the interface. From the very similar value of the differential capacitance obtained in these calculations (figure 4), we confirm that, for planar electrodes, the constant charge simulations compare well with constant potential simulations (see figure 4). We recall here that the present work is directed at validating techniques for a coarse-grained model for [BMI][PF 6 ] for future use in simulations with porous electrodes, in the latter the constant potential treatment, which allows the electrode charges to be determined self-consistently, is crucial [21]. We can thus conclude, in view of the present results, that the constant potential method has not introduced any new aspects to the results obtained at a planar electrode, despite the fact that it implicitly treats the graphite electrode as metallic.

V. STRUCTURAL PROPERTIES AT THE LIQUID-GRAPHITE INTERFACE

We showed in the previous section that the four-site model of Roy and Maroncelli [22] is able to reproduce accurately the differential capacitance of the [BMI][PF 6 ]-graphite system.

Moreover, the ionic densities predicted by this model (see figure 5) have the same characteristics as the ones obtained by all-atom simulations [17,19]. Specifically, high density peaks are observed near the planar graphite electrodes and layers can be defined, centered at approximately 5 Å, 10 Å and 15 Å from the surface at ∆Ψ 0 = 0.0 V. When the potential difference increases, the layering is modified, and an alternation of anionic and cationic layers appears, which corresponds to an overscreening effect [10].

Kislenko et al. report on a two-dimensional structure present at the interface between the uncharged graphite electrodes and the ionic liquid [17]. This type of ordering in the first layer is confirmed experimentally for [BMI][PF 6 ] on mica [47] and was also demonstrated at the interface between LiCl and aluminium electrodes [34,35]. A two-dimensional ordering in the first layer should be reflected by a sharp-peak in the in-plane partial anion-anion and cation-cation structure factors.

The PF - 6 -PF - 6 and BMI + -BMI + (for BMI + , the central site is considered, which corresponds to the imidazolium ring) partial structure factors calculated for the ions in the first adsorbed layer and corresponding snapshots are given in figure 6. When there is no potential difference, it is clear from the partial structure factors that there is a two-dimensional ordering at the interface, not present in the bulk. The first layer is organized in the same defective hexagonal lattice as for the all-atom model [17], with a lattice spacing of approximately 10 Å.

According to the snapshots and partial structure factors, the first layer becomes disordered When the potential difference is not equal to 0.0 V, an asymmetry can be seen for the positive and negative electrodes.

when a potential difference is applied between the electrodes.

VI. CONCLUSION

A coarse-grained model developed by Roy and Maroncelli [22] to represent the ionic liquid

[BMI][PF 6 ] has been tested for the prediction of interfacial properties. It appears that the model is able to reproduce accurately the liquid-vapor interface as demonstrated by the good agreement between the calculated surface tension (39.4 mN.m -1 at 400 K) and the experimental one (40.8 mN.m -1 at 393 K [START_REF] Ghatee | Fluid Phase Equilibr[END_REF], 41.07 mN.m -1 at 343.15 K [START_REF] Freire | [END_REF]). The model will provide information to go beyond mean-field approaches [48] and recent simulations [49] which propose theories to explain the anomalous capacitance increase in these materials.

6

  ions confined between planar graphite electrodes using the following conditions. A typical simulation cell contains a slice of liquid [BMI][PF 6 ] (320 ion pairs) between two graphite electrodes (see figure 2). The graphite electrodes are made of three planes of 416 atoms each with a constant distance between atoms equal to d CC = 1.43 Å. The simulation cell is orthorhombic with dimensions 32.2 Å × 34.4 Å × 123.2 Å.

Figure 2 :

 2 Figure 2: Snapshot of the simulation system: ionic liquid confined between graphite walls.

Figure 3 :

 3 Figure 3: Poisson potential profiles across the simulation for three of the nine potential differences applied. The potential in the bulk region is constant and noted Ψ bulk .

Figure 4 :

 4 Figure4: Surface charge as a function of the potential drop across the interface ∆Ψ = Ψ ± -Ψ bulk . Two linear régimes are observed and differential capacitances for the positive and negative electrode are given by the corresponding slopes. Linear regressions are extrapolated to show that the simulated points for ∆Ψ 0 = 6.0 V do not lie on these lines. Points obtained with a constant charge of ± 4.5 µC.cm -2 compare well with constant potential results.

Figure 5 :

 5 Figure 5: Ionic densities of anions and mass centers of cations for ∆Ψ 0 = 0.0 V and ∆Ψ 0 = 4.0 V.

Figure 6 :

 6 Figure 6: Upper panel: PF - 6 -PF - 6 and BMI + -BMI + partial structure factors for the ions in the first layer. For BMI + , the site considered in the calculation is the central site corresponding to the imidazolium ring. Lower panel: Corresponding representative snapshots, anions are in green, cations are in red and the graphite electrode is represented in blue. Black lines are added in the free potential snapshot to ease the visualization of the hexagonal lattice.

der Waals attraction and repulsion are !"# $% &' ( )% *+ % *,% *-% .%

  Figure 1: Schematic representation of the coarse-grained model of [BMI][PF 6 ].

	Site x	y	z	M	σ i	ε i	q i
	( Å) ( Å) ( Å) (g.mol -1 ) ( Å) (kJ.mol -1 ) (e)
	C1 0 -0.527 1.365 67.07 4.38	2.56	0.4374
	C2 0 1.641 2.987 15.04 3.41	0.36	0.1578
	C3 0 0.187 -2.389 57.12 5.04	1.83	0.1848
	A 0	0	0	144.96 5.06	4.71	-0.7800
	Table I: Force-field parameters of the coarse-grained model of [BMI][PF 6 ] [22]. Site-site interaction
	energies are given by the sum of a Lennard-Jones potential and coulombic interactions u ij (r ij ) =
	4ε ij [(

Table II :

 II Comparison of differential capacitances obtained by different simulations. In this work, the molecular dynamics are conducted with a coarse-grained model at constant potential difference with polarizable electrodes whereas in the other publications cited, the simulations use an all-atom model and constant charge electrodes.

	Kislenko et al. [17]	3.7	4.6
	Feng et al. [19]	4.40	4.09

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR) under grant ANR-2010-BLAN-0933-02 ("Modeling the Ion Adsorption in Carbon Micropores").