Anh The

Thierry Pham

Martin Jéron
email: martin.quinson@irisa.fr

Quinson

Unfolding-based Dynamic Partial Order Reduction of Asynchronous Distributed Programs

Keywords: partial order, unfolding, distributed program, asynchronous

Unfolding-based Dynamic Partial Order Reduction (UDPOR) is a recent technique mixing Dynamic Partial Order Reduction (DPOR) with concepts of concurrency such as unfoldings to efficiently mitigate state space explosion in model-checking of concurrent programs. It is optimal in the sense that each Mazurkiewicz trace, i.e. a class of interleavings equivalent by commuting independent actions, is explored exactly once. This paper shows that UDPOR can be extended to verify asynchronous distributed applications, where processes both communicate by messages and synchronize on shared resources. To do so, a general model of asynchronous distributed programs is formalized in TLA+. This allows to define an independence relation, a main ingredient of the unfolding semantics. Then, the adaptation of UDPOR, involving the construction of an unfolding, is made efficient by a precise analysis of dependencies. A prototype implementation gives promising experimental results.

Introduction

Developing distributed applications that run on parallel computers and communicate by message passing is hard due to their size, heterogeneity, asynchronicity and dynamicity. Besides performance, their correction is crucial but very challenging due to the complex interactions of parallel components.

Model-checking (see e.g. [START_REF] Baier | Principles of model checking[END_REF]) is a set of techniques allowing to verify automatically and effectively some properties on models of such systems. The principle is usually to explore all possible behaviors (states and transitions) of the system model. However, state spaces increase exponentially with the number of concurrent processes. Unfoldings and Partial order reduction (POR) are two candidate alternative techniques born in the 90's to mitigate this state space explosion and scale to large applications.

Unfoldings (see e.g. [START_REF] Esparza | Unfoldings -A Partial-Order Approach to Model Checking[END_REF]) is a concept of concurrency theory providing a representation of the behaviors of a model in the form of an event structure aggregating causal dependencies or concurrency between events (occurrence of actions), and conflicts that indicate choices in the evolution of the program. This representation may be exponentially more compact than an interleaving semantics, while still allowing to verify some properties, such as safety.

POR comprises a set of exploration techniques (see e.g. [START_REF] Godefroid | Partial-Order Methods for the Verification of Concurrent Systems -An Approach to the State-Explosion Problem[END_REF]), sharing the idea that, to detect deadlocks (and, by extension, for checking safety properties) it is sufficient to cover each Mazurkiewicz trace, i.e. a class of interleavings equivalent by commutation of consecutive independent actions. This state space reduction is obtained by choosing at each state, based on the independence of actions, only a subset of actions to explore (ample, stubborn or persistent sets, depending on the method), or to avoid (sleep set). Dynamic partial order reduction (DPOR) [START_REF] Flanagan | Dynamic partial-order reduction for model checking software[END_REF] was later introduced to combat state space explosion for stateless model-checking of software. In this context, while POR relies on a statically defined and imprecise independence relation, DPOR may be much more efficient by dynamically collecting it at run-time. Nevertheless, redundant explorations, named sleep-set blocked (SSB) [START_REF] Abdulla | Optimal dynamic partial order reduction[END_REF], may still exist that would lead to an already visited interleaving, and detected by using sleep sets.

In the last few years, two research directions were investigated to improve DPOR. The first one tries to refine the independence relation: the more precise, the less Mazurkiewicz traces exist, thus the more efficient could be DPOR. For example [START_REF] Albert | Constrained Dynamic Partial Order Reduction[END_REF] proposes to consider conditional independence relations where commutations are specified by constraints, while in [START_REF] Aronis | Optimal dynamic partial order reduction with observers[END_REF] independence is built lazily, conditionally to future actions called observers. The other direction proposes alternatives to persistent sets, in order to minimize the number of explored interleavings. Optimality is obtained when exactly one interleaving per Mazurkiewicz trace is built. In [START_REF] Abdulla | Optimal dynamic partial order reduction[END_REF] authors propose source sets that outperform DPOR, but optimality requires expensive computations of wake-up trees to avoid SSB explorations. In [START_REF] Rodríguez | Unfolding-based partial order reduction[END_REF] the authors propose unfolding-based DPOR (UDPOR), an optimal DPOR method combining the strengths of PORs and unfoldings with the notion of alternatives. The approach is generalized [START_REF] Nguyen | Quasioptimal Partial Order Reduction[END_REF] with a notion of kpartial alternative allowing to balance between optimal DPOR and sometimes more efficient sub-optimal DPOR. Some approaches already try to use DPOR techniques for the verification of asynchronous distributed applications, such as MPI programs (Message Passing Interface). In the absence of model, determining global states of the system and checking equality [START_REF] Pham | Verifying MPI applications with simgridmc[END_REF] are already challenging. In [START_REF] Palmer | Semantics driven dynamic partialorder reduction of MPI-based parallel programs[END_REF], an approach is taken that is tight to MPI. A significant subset of MPI primitives is considered, formally specified in order to define the dependency relation, and then to use the DPOR technique of [START_REF] Flanagan | Dynamic partial-order reduction for model checking software[END_REF]. In [START_REF] Sharma | A sound reduction of persistentsets for deadlock detection in MPI applications[END_REF], the efficiency is improved by focusing on particular deadlocks, but at the price of incompleteness.

We propose first steps to adapt UDPOR for asynchronous distributed applications. In [START_REF] Rosa | A simple model of communication APIsapplication to Dynamic Partial Order Reduction[END_REF] authors proposed an abstract model of distributed applications with a small set of primitives, sufficient to express most communication actions.

We revise and extend this model with synchronization primitives and formally specify it in TLA+ [START_REF] Lamport | Specifying Systems, The TLA+ Language and Tools for Hardware and Software Engineers[END_REF]. A clear advantage of this model is its abstraction: it remains concise, but its generality allows e.g. the encoding of MPI primitives. Already defining a correct independence relation from this formal model is difficult, due to the variety and complex semantics of actions. In addition, making UDPOR and in particular the computation of unfoldings and extensions efficient cannot directly rely on solutions of [START_REF] Nguyen | Quasioptimal Partial Order Reduction[END_REF], which are tuned for concurrent programs with only mutexes, thus clever algorithms need to be designed. For now we prototyped our solutions in a simplified context, but we target the Sim-Grid tool which allows to run HPC code (in particular MPI) in a simulation environment [START_REF] Degomme | Simulating MPI applications: The SMPI approach[END_REF]. The paper is organized as follows. Section 2 recalls notions of interleaving and concurrency semantics, and how a transition system is unfolded into an event structure with respect to an independence relation. In Section 3 the programming model is presented together with a sketch of the independence relation. Section 4 presents the UDPOR algorithm, its adaptation to our programming model, and how to make it efficient. Finally we present a prototype implementation and its experimental evaluation.

Interleaving and Unfolding Semantics

The behaviors of a distributed program can be described in an interleaving semantics by a labelled transition system, or in a true concurrency semantics by an event structure. An LTS equipped with an independence relation can be unfolded into an event structure [START_REF] Rodríguez | Unfolding-based partial order reduction[END_REF]. This is a main step for UDPOR.

Definition 1 (Labelled transition system). A labelled transition system (LTS) is a tuple T = S , s 0 , Σ, → where S is the set of states, s 0 ∈ S the initial state, Σ is the alphabet of actions, and →⊆ S × Σ × S is the transition relation.

We note s a -→ s when (s, a, s ) ∈ → and extend the notation to execution sequences:

s a1•a2•••an ------→ s if ∃s 0 = s, s 1 , . . . s n = s with s i-1 a -→ i s i for i ∈ [1, n].
For a state s in S , we denote by enabled (s) = {a ∈ Σ : ∃s ∈ S , s a -→ s } the set of actions enabled at s.

Independence is a key notion in both POR techniques and unfoldings, linked to the possibility to commute actions: Definition 2 (Commutation and independence). Two actions a 1 , a 2 of an LTS T = S , s 0 , Σ, → commute in a state s if they satisfy the two conditions:

executing one action does not enable nor disable the other one:

a 1 ∈ enabled (s) ∧ s a1 -→ s =⇒ (a 2 ∈ enabled (s) ⇔ a 2 ∈ enabled (s )) (1)
their execution order does not change the overall result:

a 1 , a 2 ∈ enabled (s) =⇒ (s a1•a2 ---→ s ∧ s a2•a1 ---→ s =⇒ s = s ) (2)
A relation I ⊆ Σ × Σ is a valid independence relation if it under-approximates commutation, i.e. ∀a 1 , a 2 , I (a 1 , a 2 ) implies that a 1 and a 2 commute in all states. Conversely a 1 and a 2 are dependent and we note D(a 1 , a 2 ) when ¬(I (a 1 , a 2 )).

A Mazurkiewicz trace is an equivalence class of executions (or interleavings) of an LTS T obtained by commuting adjacent independent actions. By the second item of Definition 2, all these interleavings reach a unique state. The principle of all DPOR approaches is precisely to reduce the state space exploration while covering at least one execution per Mazurkiewicz trace. If a deadlock exists, a Mazurkiewicz trace leads to it and will be discovered. More generally, safety properties are preserved.

The UDPOR technique that we consider also uses concurrency notions. A classical model of true concurrency is prime event structures: Definition 3 (Prime event structure). Given an alphabet of actions Σ, a Σ-prime event structure (Σ-PES) is a tuple E = E , <, #, λ where E is a set of events, < is a partial order relation on E , called the causality relation, λ : E → Σ is a function labelling each event e with an action λ(e), # is an irreflexive and symmetric relation called the conflict relation such that, the set of causal predecessors or history of any event e, e = {e ∈ E : e < e} is finite, and conflicts are inherited by causality: ∀e, e , e ∈ E , e#e ∧ e < e =⇒ e#e .

Intuitively, e < e means that e must happen before e , and e#e that those two events cannot belong to the same execution. Two distinct events that are neither causally ordered nor in conflict are said concurrent. The set [e] := e ∪ {e} is called the local configuration of e. An event e can be characterized by a pair < λ(e), H > where λ(e) is its action, and H = e its history.

We note conf (E ) the set of configurations of E, where a configuration is a set of events C ⊆ E that is both causally closed (e ∈ C =⇒ e ⊆ C ) and conflict free (e, e ∈ C =⇒ ¬(e#e )). A configuration C is characterized by its causally maximal events maxEvents(C ) = {e ∈ C : e ∈ C , e < e }, since it is exactly the union of local configurations of these events: C = e ∈ maxEvents(C ) [e]; conversely a conflict free set K of incomparable events for < defines a configuration config(K ) and C = config(maxEvents(C )).

A configuration C , together with the causal and independence relations defines a Mazurkiewicz trace: all interleavings are obtained by causally ordering all events in the configuration C but commuting concurrent ones. The state of a configuration C denoted by state(C ) is the state in T reached by any of these executions, and it is unique as discussed above. We write enab(C ) = enabled (state(C )) ⊆ Σ for the set of actions enabled at state(C ), while actions(C ) denotes the set of actions labelling events in C , i.e. actions(C ) = {λ(e) : e ∈ C }.

The set of extensions of C is ex (C ) = {e ∈ E \ C : e ⊆ C }, i.e. the set of events not in C but whose causal predecessors are all in C . When appending an extension to C , only resulting conflict-free sets of events are indeed configurations. These extensions constitute the set of enabled events en(C ) = {e ∈ ex (C ) : e ∈ C , e#e } while the other ones are conflicting extensions cex (C ) := ex (C ) \ en(C ).

Parametric Unfolding Semantics. Given an LTS T and an independence relation I , one can build a prime event structure E such that each linearization of a maximal (for inclusion) configuration represents an execution in T , and conversely, to each Mazurkiewicz trace in T corresponds a configuration in E [START_REF] Nguyen | Quasioptimal Partial Order Reduction[END_REF].

Definition 4 (Unfolding). The unfolding of an LTS T under an independence relation I is the Σ-PES E = E , <, #, λ, incrementally constructed from the initial Σ-PES ∅, ∅, ∅, ∅ by the following rules until no new event can be created:

-for any configuration C ∈ conf (E ),

Programming Model and Independence Relation

In this section we introduce the abstract model of asynchronous distributed systems that we consider. While abstract, this model is sufficient to represent concrete MPI programs, as it encompasses all building blocks of the SMPI implementation of the standard [START_REF] Degomme | Simulating MPI applications: The SMPI approach[END_REF]. We formalized this model in the specification language TLA+ [START_REF] Lamport | Specifying Systems, The TLA+ Language and Tools for Hardware and Software Engineers[END_REF], to later infer an independence relation. In our model an asynchronous distributed system P consists in a set of n actors Actors = {A 1 , A 2 , ...A n } that perform local actions, communicate asynchronously with each others, and share some resources. We assume that the program is terminating, which implies that all actions are terminating. All local actions are abstracted into a unique one LocalComp. Communication actions are of four types: AsyncSend, AsyncReceive, TestAny , and WaitAny . Actions on shared resources called synchronizations are of four types: AsyncMutexLock, MutexUnlock, MutexTest and MutexWait.

Abstract Model

At the semantics level, P is a tuple P = Actors, Network, Synchronization where Network and Synchronization respectively describe the abstract objects, and the effects on these of the communication and synchronizations actions. The Network subsystem provides facilities for the Actors to asynchronously communicate with each others, while the subsystem Synchronization allows the synchronization of actors on access to shared resources.

Network subsystem. The state of the Network subsystem is defined as a pair Mailboxes, Communications , where Mailboxes is a set of mailboxes storing unpaired communications, while Communications stores only paired ones. Each communication c has a status in {send, receive, done}, ids of source and destination actors, data addresses for those. A mailbox is a rendez-vous point where send and receive communications meet. It is modelled as an unbounded FIFO queue that is either empty, or stores communications with all same send or receive status, waiting for a matching opposite communication. When matching occurs, this paired communication gets a done status and is appended to the set Communications. We now detail the effect in actor A i of the communication actions on Mailboxes and Communications: Synchronization subsystem. The Synchronization subsystem consists in a pair Mutexes, Requests where Mutexes is a set of asynchronous mutexes used to synchronize the actors, and Requests is a vector indexed by actors ids of sets of requested mutexes. Each mutex m j is represented by a FIFO queue of actors ids i who declared their interest on a mutex m j by executing the action AsyncMutexLock(m j ). A mutex m j is free if its queue is empty, busy otherwise. The owner is the actor whose id is the first in the queue. In actor A i , the effect of the synchronization actions on Mutexes and Requests is as follows:

-c = AsyncSend(m,
-AsyncMutexLock(m j ) requests a mutex m j with the effect of appending the actor id i to m j 's queue and adding j to Requests[i ]. A i is waiting until owning m j but, unlike classical mutexes, waiting is not necessarily blocking. -MutexUnlock(m j ) removes its interest to a mutex m j by deleting the actor id i from the m j 's queue and removing j from

Requests[i ]. -MutexTest(M ) returns true if actor A i owns some previously requested mutex m j in M (i is first in FIFO m j ∈ M s.t. j in Requests[i ]).
-MutexWait(M ) blocks until A i owns some mutex m j in M . Note that Mutex-Test (resp. MutexWait) are similar to TestAny (resp. WaitAny ) and could be merged. We keep them separate here for simplicity of explanations.

Beside those actions, a program can have local computations named Local-Comp actions. Such actions do not intervene with shared objects (Mailboxes, Mutexes and Communications), and they can be responsible for I/O tasks.

We specified our model of asynchronous distributed systems in the formal language TLA+ [START_REF] Lamport | Specifying Systems, The TLA+ Language and Tools for Hardware and Software Engineers[END_REF]. Our TLA+ model1 focuses on how actions transform the global state of the system. An instance P of a program is described by a set of actors and their actions (representing their source code). Following the semantics of TLA+, and since programs are terminating, the interleaving semantics of a program P can be described by an acyclic LTS representing all its behaviors. Formally, the LTS of P is a tuple T P = S , s 0 , Σ, → where Σ represent the actions of P ; a state s =< l , g > in S consists of the local state l of all actors (i.e. local variables, Requests) and g the state of all shared objects including Mutexes, Mailboxes and Communications; in the initial state s 0 all actors are in their initial local state, sets and FIFO queues are empty; a transition s a -→ s is defined if, according to the TLA+ model, the action encoded by a is enabled at s and executing a transforms the state from s to s .

Notice that when verifying a real program, we only observe its actions and assume that they respect the proposed TLA+ model and the independence relation discussed below. These assumptions are necessary to suppose that the LTS correctly models the actual program behaviors.

Additional Property of the Model

The model presented in the previous section may appear unusual, because the lock action on mutexes is split into an AsyncMutexLock and a MutexWait while most works in the literature consider atomic locks. Our model does not induce any loss of generality, since synchronous locks can trivially be simulated with asynchronous locks. One reason to introduce this specificity is that this entails the following lemma, that is the key to the efficiency of UDPOR in our model.

Lemma 1 (Persistence). Let u be a prefix of an execution v of a program in our model. If an action a is enabled after u, it is either executed in v or still enabled after v .

Intuitively, persistence says that once enabled, actions are never disabled by any subsequent action, thus remain enabled until executed. Persistence does not hold for classical synchronous locks, as some enabled lock(m) action of an actor may become disabled by the lock(m) of another actor. This persistence property has been early introduced by Karp and Miller [START_REF] Karp | Parallel program schemata[END_REF], and later studied for Petri Nets [START_REF] Landweber | Properties of conflict-free and persistent Petri Nets[END_REF]. It should not be confused with the notion of persistent set used in DPOR 2 . Persistent sets are linked to independence, while persistence is not.

Proof. When a is a LocalComp, AsyncSend, AsyncReceive, TestAny , AsyncMutex-Lock, MutexUnlock, or MutexTest action, a cannot be disabled by any new action. Indeed, these actions are never blocking (e.g. AsyncMutexLock comes down to the addition of an element in a FIFO, which is always enabled) and only depend on the execution of the action right before them by the same actor.

WaitAny and MutexWait may seem more complex. If a is a WaitAny , being enabled after u means that one communication it refers to was paired. Similarly, if a is a MutexWait, being enabled after u means that the corresponding actor is first in the FIFO of a mutex it refers to. In both cases these facts cannot be modified by any subsequent action, so a remains enabled until executed.

Independence Theorems

In order to use DPOR algorithms for our model of distributed programs, and in particular UDPOR that is based on the unfolding semantics, we need to define a valid independence relation for this model. Intuitively, two actions in distinct actors are independent when they do not compete on shared objects, namely Mailboxes, Communications, or Mutexes. This relation is formally expressed in TLA+ as so-called "independence theorems". We use the term "theorem" since the validity of the independence relation with respect to commutation should be proved. We proved them manually and implemented them as rules in the model-checker. These independence theorems are as follows 3 :

1. A LocalComp is independent with any other action of another actor. 2. Any synchronization action is independent of any communication action of a distinct actor. 3. Any pair of communication actions in distinct actors concerning distinct mailboxes are independent. 4. An AsyncSend is independent of an AsyncReceive of another actor.

2 A set of transitions T is called persistent in a state s if all transitions not in T and, either enabled in T or enabled in a state reachable by transitions not in T , are independent with all transitions in T . As a consequence, exploring only transitions in persistent sets is sufficient to detect all deadlocks. 3 Some independence theorems could be enlarged but we give these ones for simplicity.

Any pair of actions in

Adapting UDPOR

This section first recalls the UDPOR algorithm of [START_REF] Rodríguez | Unfolding-based partial order reduction[END_REF] and then explains how it may be adapted to our context, in particular how the computation of extensions, a key operation, can be made efficient in our programming model.

The UDPOR Algorithm

Algorithm 1 presents the UDPOR exploration algorithm of [START_REF] Rodríguez | Unfolding-based partial order reduction[END_REF]. Like other DPOR algorithms, it explores only a part of the LTS of a given terminating distributed program P according to an independence relation I , while ensuring that the explored part is sufficient to detect all deadlocks. The particularity of UDPOR is to use the concurrency semantics explicitly, namely unfoldings, which makes it both complete and optimal: it explores exactly one interleaving per Mazurkiewicz trace, never reaching any sleep-set blocked execution. The algorithm works as follows. Executions are represented by configurations, thus equivalent to their Mazurkiewicz traces. The set U , initially empty, contains all events met so far in the exploration. The procedure Explore has three parameters: a configuration C encoding the current execution; a set D (for disabled) of events to avoid (playing a role similar to a sleep set in [START_REF] Godefroid | Partial-Order Methods for the Verification of Concurrent Systems -An Approach to the State-Explosion Problem[END_REF]), thus preventing revisits of configurations; a set A (for add) of events conflicting with D and used to guide the search to events in conflicting configurations in cex (C ) to explore alternative executions.

First, all extensions of C are computed and added to U (line 4). The search backtracks (line 6) in two cases: when C is maximal (en(C ) = ∅), i.e. a deadlock (or the program end) is reached, or when all events enabled in C should be avoided (en(C ) ⊆ D), which corresponds to a redundant call, thus a sleepset blocked execution. Otherwise, an enabled event e is chosen (line 7-10), in

Explore(C ∪ {e}, D, A \ {e}) 12 if ∃J ∈ Alt(C , D ∪ {e}) then 13 Explore(C , D ∪ {e}, J \ C ) 14 U := U ∩ QC,D
A if this guiding information is non empty (line 10), and a "left" recursive exploration Explore(C ∪ {e}, D, A \ {e}) is called (line 11) from this extended configuration C ∪ {e}, it continues trying to avoid D, but e is removed from A in the guiding information. When this call is completed, all configurations containing C and e have been explored, thus it remains to explore those that contain C but not e. In this aim alternatives are computed (line 12) with the function call Alt(C , D ∪ {e}). Alternatives play a role similar to "backtracking sets" in the original DPOR algorithm, i.e. sets of actions that must be explored from the current state. Formally, an alternative to D = D ∪ {e} after C in U is a subset J of U that, does not intersect D , forms a configuration C ∪ J after C , and such that all events in D conflict with some event in J . If an Alternatives J exists, a right "recursive" exploration is called Explore(C , D ∪ {e}, J \ C ): C is still the configuration to extend, but e is now also to be avoided, thus added to D, while events in J \ C are used as guides. Upon completion (line 14), U is intersected with Q C ,D which includes all events in C and D as well as every event in U conflicting with some events in C ∪ D.

In order to avoid sleep-set blocked executions (SSB) and obtain the optimality of DPOR, the function Alt(C , D ∪{e}) has to solve an NP-complete problem [START_REF] Nguyen | Quasioptimal Partial Order Reduction[END_REF]: find a subset J of U that can be used for backtracking, conflicts with all D ∪ {e} thus necessarily leading to a configuration C ∪J that is not already visited. In this case en(C ) ⊆ D can then be replaced by en(C ) = ∅ in line 5. Note that with a different encoding, Optimal DPOR must solve the same problem [START_REF] Abdulla | Optimal dynamic partial order reduction[END_REF] as explained in [START_REF] Nguyen | Quasioptimal Partial Order Reduction[END_REF]. In [START_REF] Nguyen | Quasioptimal Partial Order Reduction[END_REF], a variant of the algorithm is proposed for the function Alt that computes k -partial alternatives rather than alternatives, i.e. sets of events J conflicting with only k events in D, not necessarily all of them. Depending on k , (e.g. k = ∞ (or k = |D| + 1) for alternatives, k = 1 for source sets of [START_REF] Abdulla | Optimal dynamic partial order reduction[END_REF]) this variant allows to tune between an optimal or a quasi-optimal algorithm that may be more efficient.

Computing Extensions Efficiently

Computing the extensions ex (C ) of a configuration C may be costly in general. It is for example an NP-complete problem for Petri Nets since all subconfigurations must be enumerated. Fortunately this algorithm can be specially tuned for sub-classes. In particular for the programming model of [START_REF] Rodríguez | Unfolding-based partial order reduction[END_REF][START_REF] Nguyen | Quasioptimal Partial Order Reduction[END_REF] it is tuned in an algorithm working in time O(n 2 log(n)), using the fact that events have a maximum of two causal predecessors, thus limiting the subsets to consider.

This section tunes the algorithm to our more complex model, using the fact that the amount of causal predecessors of events is also bounded. Next section shows how to incrementally compute ex (C ) to avoid recomputations. Figure 2 illustrates some aspects of an extension. This section mandates some additional notations. Given a configuration C and an extension with action a, let pre(a) denote the action right before a in the same actor, while preEvt(a, C ) denotes the event in C associated with pre(a) (formally e = preEvt(a, C ) ⇐⇒ e ∈ C , λ(e) = pre(a)). Given a set F of events F ⊆ E , Depend (a, F ) means that a depends on all actions labeling events in F .

The definition of ex (C ) (set of extensions of a configuration C ) {e ∈ E \ C : e ⊆ C } can be rewritten using the definitions of section 2 as follows:

{e = a, H ∈ E \ C : a = λ(e) ∧ H = e ∧ H ∈ 2 C ∩ conf (E ) ∧ a ∈ enab(H )}.
Fortunately, it is not necessary to enumerate all subsets H of C , that are in exponential numbers, to compute this set. According to the unfolding construction in Definition 4, an event e = a, H only exists in ex (C ) if the action a is dependant with the actions of all maximal events of H . This gives:

ex (C ) = {e = a, H ∈ E \ C : a = λ(e) ∧ H = e ∧ H ∈ 2 C ∩ conf (E ) ∧ a ∈ enab(H ) ∧ Depend (a, maxEvents(H ))}. Now ex (C )
can be simplified and decomposed by enumerating Σ, yielding to:

ex (C ) = a ∈ Σ { a, H : H ∈ S a,C } \ C where S a,C = {H ∈ conf (E ) : H ⊆ C ∧ a ∈ enab(H ) ∧ Depend (a, maxEvents(H ))}.
The above formulation of ex (C ) iterates on all actions in Σ. However, interpreting the persistence property (Lemma 1) for configurations entails that for two configurations H and C with H ⊆ C , an action a in enab(H ) is either in actions(C ) or enab(C ).

Therefore, ex (C ) can be rewritten by restricting a to actions(C ) ∪ enab(C ) :

ex (C ) = ( a ∈ actions(C )∪enab(C ) { a, H : H ∈ S a,C }) \ C (3)
Now, instead of enumerating possible configurations H ∈ S a,C , we can enumerate their maximal sets K = maxEvents(H ). Hence,

ex (C ) = ( a ∈ actions(C )∪enab(C ) { a, config(K ) : K ∈ S max a,C }) \ C (4)
with

S max a,C = {K ∈ 2 C : K is maximal ∧a ∈ enab(config(K ))∧Depend (a, K ))} and K is maximal if ( e, e ∈ K : e < e ∨ e#e ).
One can then specialize the computation of ex (C ) according to the type of action a. Due to space limitations, we only detail the computation for AsyncSend actions, the other ones being similar.

Computing extensions for AsyncSend actions. Let C be a configuration, and a an action of type c = AsyncSend(m, ) of an actor A i . We want to compute the set S max a,C of sets K of maximal events from which a depends. According to independence theorems (see 3.3), a only depends on the following actions: pre(a), all AsyncSend(m, ) actions of distinct actors A j which concern the same mailbox m, and all WaitAny (resp. TestAny ) actions that wait (resp. test) a AsyncReceive which concerns the same communication c. Considering this, we now examine the composition of maximal events sets K in S max a,C . First, two events labelled by AsyncSend(m, ) actions cannot co-exist in K , formally e, e ∈ K : λ(e), λ(e ) in AsyncSend(m, ): indeed, if two such events exist in a configuration, they are dependent but cannot conflict, thus are causality related and cannot be both maximal.

Second, if a WaitAny (Com) action concerns communication c, there are two cases: (i) either c is not the first done communication in Com, then WaitAny (Com) and the action a are independent. (ii) or c is the first done communication in Com and WaitAny is enabled only after a. Thus the only possibility for a maximal event to be labelled by a WaitAny is when pre(a) is a WaitAny of the same actor. We can then write: e ∈ K : λ(e) in WaitAny ∧λ(e) = pre(a).

Third, all AsyncReceive events for the mailbox m are causally related in configuration C , and c can only be paired with one of them, say c . Thus a can only depend on actions TestAny (Com ) such that c ∈ Com and c and c form the first done communication in Com , and all those TestAny events are ordered. Thus, there is at most one event e labelled by TestAny in K such that λ(e) = pre(a).

To conclude, K contains at most three events: preEvt(a, C ), some event labelled with an action AsyncSend on the same mailbox, and some TestAny for some matching AsyncReceive communication. There is thus only a cubic number of such sets, which is the worse case among considered action types. Algorithm 2 generates all events in ex (C ) labelled by an AsyncSend action a. ). We then iterate on all AsyncSend events in C , combining them with e 2 to create maximal event sets K (lines 2-8). We only have one AsyncSend event e 3 . Since ¬(e 2 < e 3 ) and ¬(e 3 < e 2 ), we form a first set K = {e 2 , e 3 }, and add e 7 =< AsyncSend, {e 2 , e 3 } > to ex (C ). Next all TestAny events that concern the mailbox m should be considered. Events e 2 and e 5 can be combined to form a new maximal event set K = {e 2 , e 5 }, but since a and λ(e 5 ) are not related to the same communication, D(a, λ(e 5 )) is not satisfied and no event is created. Finally combinations of e 2 with an AsyncSend event and a TestAny event are examined (lines 9-17). We then get K = {e 2 , e 5 , e 3 }, and e 8 is added to ex (C ) since D(a, λ(e 5 )) holds in the configuration config({e 2 , e 5 , e 3 }).

Computing Extensions Incrementally

In the UDPOR exploration algorithm, after extending a configuration C by adding a new event e, one must compute the extensions of C = C ∪ {e}, thus resulting in redundant computations of events. The next theorem improves this by providing an incremental computation of extensions.

Theorem 1. Suppose C = C ∪ {e} where e is the last event added to C by the Algorithm 1. We can compute ex (C ) incrementally as follows: Proof. With the definition of S a,C as above, recall that

ex (C ) = (ex (C ) ∪ a ∈ enab(C ) {< a, H >: H ∈ S a,C }) \ {e} (5)
ex (C ) = ( a ∈ actions(C )∪enab(C ) { a, H : H ∈ S a,C }) \ C (6)
Applying the same equation ( 6) to C we get: (⊆) Suppose there exists some event e n =< a, H > belonging to the left but not the right set. If a = λ(e n ) = λ(e), then H ∈ S a,C ∩S a,C , so e n is in both sets, resulting in a contradiction. If a = λ(e n ) = λ(e), there are two cases: (i) either e / ∈ H then H ∈ S a,C and e n belongs to the right set, a contradiction. (ii) or e ∈ H , then λ(e n ) ∈ actions(C ) \ {λ(e)} = actions(C ), thus there is another event e ∈ C such that λ(e ) = λ(e n ), then e cannot belong to H (one action a cannot appear twice in e n ). Besides, e is the last event explored in C , thus a depends on λ(e) by Definition 4. Then, e conflicts with e, contradicting their membership to the same configuration C . This proves [START_REF] Flanagan | Dynamic partial-order reduction for model checking software[END_REF].

ex (C ) = ( a ∈ actions(C )∪enab(C ) { a, H : H ∈ S a,C }) \ C Now,
Second, since C ⊆ C , according to persistence of the programming model (Lemma 1), (enab(C ) \ {λ(e)}) ⊆ enab(C ). We thus have:

a ∈ enab(C )\{λ(e)} {< a, H >| H ∈ S a,C } ⊆ a ∈ enab(C ) {< a, H >| H ∈ S a,C } (8)
Now, using equations ( 7) and ( 8), ex (C ) can be rewritten as follows:

ex (C ) = ( a ∈ actions(C )∪(enab(C )\λ(e)) { a, H : H ∈ S a,C } ∪ a ∈ enab(C ) { a, H : H ∈ S a,C }) \ (C ∪ {e}) (9)
But since no event in a ∈ enab(C ) { a, H : H ∈ S a,C } is in (C ∪ {e}), equation ( 9) can be rewritten as Equation ( 5) in Theorem 1.

Experiments

We implemented the quasi-optimal version of UDPOR with k -partial alternatives [START_REF] Nguyen | Quasioptimal Partial Order Reduction[END_REF] in a prototype adapted to the distributed programming model of Section 3, i.e. with its independence relation. The computation of k -partial alternatives is essentially inspired by [START_REF] Nguyen | Quasioptimal Partial Order Reduction[END_REF]. Recall the algorithm reaches optimality when k = |D| + 1, while k = 1 corresponds to Source DPOR [START_REF] Abdulla | Optimal dynamic partial order reduction[END_REF]. The prototype is still limited, not connected to the SimGrid environment, thus can only be experimented on simple examples. We first compare optimal UDPOR with an exhaustive stateless search on several benchmarks (see Table 1). The first five benchmarks come from Umpire Tests5 , while DTG and RMQ-receiving belong to [START_REF] Khanna | Dynamic symbolic verification of MPI programs[END_REF] and [START_REF] Rosa | A simple model of communication APIsapplication to Dynamic Partial Order Reduction[END_REF], respectively. The last benchmark is an implementation of a simple Master-Worker pattern. We expressed them in our programming model and explored their state space with our prototype. The experiments were performed on an HP computer, Intel Core i7-6600U 2.60GHz processors, 16GB of RAM, and Ubuntu version 18.04.1. Table 1 presents the number of explored traces and running time for both an exhaustive search and optimal UDPOR. In all benchmarks UDPOR outperforms the exhaustive search. For example, for RMQ-receiving with 4 processes, the exhaustive search explores more than 20000 traces in around 8 seconds, while UDPOR explores only 6 traces in 0.2 second. Besides, UDPOR is optimal, exploring only one trace per Mazurkiewicz trace. For example in RMQ-receiving with 5 processes, with only 4 AsyncSend actions that concern the same mailbox, UDPOR explores exactly 24 (=4!) non-equivalent traces. Similarly, the DTG benchmark has only two dependent AsyncSend actions, thus two non-equivalent traces. Furthermore, deadlocks are also detected in the prototype.

We also tried to vary the value of k . When k is decreased, one gains in efficiency in computing alternatives, but looses optimality by producing more traces. It is then interesting to analyse, whether this can be globally more efficient than optimal UDPOR. Similar to [START_REF] Nguyen | Quasioptimal Partial Order Reduction[END_REF], we observed that in some cases, fixing smaller values of k may improve the efficiency. For example with RMQ-receiving, k = 7 is optimal, but reducing to k = 4 still produces 24 traces (thus is optimal) a bit more quickly (2.3 s), while for k = 3 the number of traces and time diverge. We have to analyse this more precisely on more examples in the future.

Note that with our simple prototype, we do not yet make experiments with concrete programs (e.g. MPI programs), for which running time may somehow diverge. We expect to make it in the next months and then experiment the algorithms in more depth. However, we believe that the results are already significant and that UDPOR is effective for asynchronous distributed programs.

Benchmarks

Conclusion and Future Work

The paper adapts the unfolding-based dynamic partial order reduction (UD-POR) approach [START_REF] Rodríguez | Unfolding-based partial order reduction[END_REF] to the verification of asynchronous distributed programs.

The programming model we consider is generic enough to properly model a large class of asynchronous distributed systems, including e.g. MPI applications, while exhibiting some interesting properties. From a formal specification of this model in TLA+, an independence relation is built, that is used by UDPOR to partly build the unfolding semantics of programs. We show that, thanks to the properties of our model, some usually expensive operations of UDPOR can be made efficient. A prototype of UDPOR has been implemented and experimented on some benchmarks, gaining promising first results.

In the future we aim at extending our model of asynchronous distributed systems, while both preserving good properties, getting a more precise independence relation, and implementing UDPOR in the SimGrid model-checker and verify real MPI applications. Once done, we should experiment UDPOR more deeply, and compare it with state of the art tools on more significant benchmarks, get a more precise analysis about the efficiency of UDPOR compared to simpler DPOR approaches, analyse the impact of quasi-optimality on efficiency.

  any action a ∈ enabled (state(C )), if for any e ∈ maxEvents(C ), ¬I (a, λ(e )), add a new event e = a, C to E ; for any such new event e = a, C , update <, # and λ as follows: λ(e) := a and for every e ∈ E \ {e}, consider three cases: (i) if e ∈ C then e < e, (ii) if e / ∈ C and ¬I (a, λ(e )), then e#e , (iii) otherwise, i.e. if e / ∈ C and I (a, λ(e )), then e and e are concurrent.

Fig. 1 .

 1 Fig. 1. Main elements of the model: Actors, Network and Synchronization

Fig. 2 .

 2 Fig.2. A configuration C , extended by event e, its history H and maximal events K .

Algorithm 2 : 4 K := ∅ 5 - 6 - 11 where 12 K := ∅ 13 -Example 1 .

 24561112131 createAsyncSendEvt(a, C ) 1 create e = < a, config(preEvt(a, C )) >, and ex (C ) := ex (C ) ∪ {e } 2 foreach e ∈ C s.t. λ(e) ∈ {AsyncSend(m, ), TestAny(Com)} 3 where Com contains a matching c = AsyncReceive(m, ) with a do if ¬(e < preEvt(a, C ) then K := K ∪ {e} if ¬(preEvt(a, C ) < e) then K := K ∪ {preEvt(a, C )} 7 if D(a, λ(e)) then 8 create e = < a, config(K ) > and ex (C ) := ex (C ) ∪ {e } 9 foreach es ∈ C s.t. λ(es ) = AsyncSend(m, ) do 10 foreach et ∈ C s.t. λ(et ) = TestAny(Com) Com contains a matching c = AsyncReceive(m, ) with a do if ¬(es < preEvt(a, C )) and ¬(es < et ) then K := K ∪ {es } 14 -if ¬(et < preEvt(a, C ) and ¬(et < es ) then K := K ∪ {et } 15 -if ¬(preEvt(a, C ) < es ) and ¬(preEvt(a, C ) < et ) then K := K ∪ {preEvt(a, C )} 16 if D(a, λ(et )) then 17 create e = < a, config(K ) >, and ex (C ) := ex (C ) ∪ {e } We illustrate the Algorithm 2 by the example of Figure 3. Suppose we want to compute the extensions of C = {e 1 , e 2 , e 3 , e 4 , e 5 } associated with a, the action c2 = AsyncSend(m, ) of Actor 2 . First e 6 =< AsyncSend, {2} > ∈ ex (C ) because preEvt(a, C ) = e 2 (line 1

Fig. 3 .

 3 Fig. 3. The pseudo-code of a distributed program (left) and the configuration C .

  exploring e from C leads to C , which entails that λ(e) belongs to enab(C ) and actions(C ) ∪ λ(e) = actions(C ), thus the range of a in ex (C ) which is actions(C ) ∪ enab(C ) can be rewritten actions(C ) ∪ (enab(C ) \ λ(e)). First, separating action(C ) from the rest in both ex (C ) and ex (C ) we prove: a ∈ actions(C ) {< a, H >: H ∈ S a,C } = a ∈ actions(C ) {< a, H >: H ∈ S a,C } (7) (⊇) This inclusion is obvious since C ⊇ C , and thus S a,C ⊇ S a,C .

  TestAny (Com) tests a set of communications Com of A i . It returns a boolean, true if and only if some communication in Com with done status exists. -WaitAny (Com) waits for a set of communications Com of A i . The action is blocking until at least one communication in Com has a done status.

data) drops an asynchronous send communication c to the mailbox m. If pending receive communications exist in the mailbox, c is paired with the oldest one c to form a communication with done status in Communications, the receive communication is removed from m and the data is copied from the source to the destination. Otherwise, a pending communication with send status is appended to m.

c = AsyncReceive(m, d ) drops an asynchronous receive communication to mailbox m; the way a receive communication is processed is similar to send.

If pending send communications exist, c is paired with the oldest one c to form a communication with done status in Communications, the send communication is removed from m, and the data of the send is copied to d . Otherwise, a pending communication with receive status is appended to m. -

  {TestAny , WaitAny } in distinct actors is independent. 6. Any action in {TestAny (Com), WaitAny (Com)} is independent with any action of another actor in {AsyncSend, AsyncReceive} as soon as they do not both concern the first paired communication in the set Com 4 . 7. Any pair of synchronization actions of distinct actors concerning distinct mutexes are independent. 8. An AsyncMutexLock is independent with a MutexUnlock of another actor. 9. Any pair of actions in {MutexWait, MutexTest} of distinct actors is independent. 10. A MutexUnlock is independent of a MutexWait or MutexTest of another actor, except if the same mutex is involved and one of the two actors owns it. 11. An AsyncMutexLock is independent of any MutexWait and MutexTest of another actor.

Table 1 .

 1 Comparing exhaustive exploration and UDPOR. TO: timeout after 30 minutes; #P: number of processes; Deadlock: deadlock exists; #Traces: number of traces

		#P Deadlock	Exhaustive search #Traces Time (second) #Traces Time (second) UDPOR
	wait-deadlock	2	yes	2	<0.01	1	<0.01
	complex-deadlock	3	yes	36	0.03	1	<0.01
	waitall-deadlock	3	yes	1458	1.2	1	<0.01
	no-error-wait-any src 3	no	21	0.02	1	0.01
	any-src-can-deadlock3 3	yes	999	0.65	2	0.03
	DTG	5	yes	-	TO	2	0.07
	RMQ-receiving	4 5	no no	20064 -	8.15 TO	6 24	0.2 2.52
	Master-worker	3 4	no no	1356444 -	1038 TO	2 6	0.2 2.5

https://github.com/pham-theanh/simixNetworks

Intuitively, WaitAny (Com) needs only one done communication (the first paired (AsyncSend, AsyncReceive)) in Com to become enabled. Similarly, the effect of TestAny (Com) only depends on this first done communication.

http://formalverification.cs.utah.edu/ISP-Tests/

Acknowledgement: we wish to thank the reviewers for their constructive comments to improve the paper.

This work has been supported by INRIA collaborative project IPL HAC-SPECIS (http://hacspecis.gforge.inria.fr/)


Keywords: water, heat capacity, molecular simulation

Molecular simulation study of the heat capacity of metastable water between 100K and 300K

Molecular simulations have been used to study the heat capacity of metastable liquid water at low temperature adsorbed on a smooth surface. These calculations aim at modelling water properties measured by experiments performed on water films adsorbed on Vycor nanoporous silica at low temperature. In particular, the study focuses on the non-monotonous variation of the heat capacity around between 100 and 300 K.

Introduction

In bulk water, local energetics and long-range connectivity come along with a very specific local three-dimensional (3D) tetrahedral organization of the hydrogen bond (HBond) network.

Altogether these properties drive the numerous so-called "water anomalies", [START_REF] Debenedetti | Supercooled and Glassy Water[END_REF][START_REF] Caupin | Escaping the No Man's Land: Recent Experiments on Metastable Liquid Water[END_REF] for example, the apparent power-law divergence of the specific heat (Cp) at low temperature.

In interfacial situations, as the water molecules are forced to organize in only two dimensions, the three dimensional extension of the HBond network is frustrated. This profoundly affects the physics of water and a number of events are experimentally detected at 150, 220 and 250 K. [START_REF] Zanotti | Competing Coexisting Phases in 2d Water[END_REF] Thermal events are detected at these three temperatures in experimental Cp measurement, [START_REF] Tombari | Dynamic and Apparent Specific Heats During Transformation of Water in Partly Filled Nanopores During Slow Cooling to 110 K and Heating[END_REF] but no clear-cut interpretation of the phenomenon has been reached so far.

Computer simulation of water is widely used to provide important information on the structure at the molecular level, as well as dynamical and thermodynamic properties. The quality of the intermolecular potential is of course mandatory for reliable results. Many water models have been developed, the simplest being rigid with fixed partial charges, to more sophisticated ones including flexible bonds and polarizability. [START_REF] Guillot | A Reappraisal of What We Have Learnt During Three Decades of Computer Simulations on Water[END_REF][START_REF] Ouyang | Modelling Water: A Lifetime Enigma[END_REF] One of the most popular models is probably the SPC model. [START_REF] Berendsen | Intermolecular Forces; B. Pullman[END_REF] It has been intensively used in molecular simulations involving large systems because it is a fast model in terms of computational cost. The heat capacity of SPC water has been shown to be in reasonable agreement with other simulation models as well as experiments. [START_REF] Jorgensen | Temperature Dependence of TIP3P, SPC, and TIP4P Water from NPT Monte Carlo Simulations: Seeking Temperatures of Maximum Density[END_REF] The water adsorption in nanoporous silica has been extensively studied by molecular simulations, in particular in the grand canonical ensemble, or the generalization of the Gibbs ensemble. [START_REF] Puibasset | Water Adsorption on Hydrophilic Mesoporous and Plane Silica Substrates: A Grand Canonical Monte Carlo Simulation Study[END_REF][START_REF] Puibasset | A Grand Canonical Monte Carlo Simulation Study of Water Adsorption on Vycor-Like Hydrophilic Mesoporous Silica at Different Temperatures[END_REF][START_REF] Puibasset | Thermodynamic Characterization of Fluids Confined in Heterogeneous Pores by Monte Carlo Simulations in the Grand Canonical and the Isobaric-Isothermal Ensembles[END_REF][START_REF] Puibasset | Phase Coexistence in Heterogeneous Porous Media : A New Extension to Gibbs Ensemble Monte Carlo Simulation Method[END_REF][START_REF] Mischler | Water Adsorption on Amorphous Silica Surfaces: A Car-Parrinello Simulation Study[END_REF][START_REF] Malani | Influence of Hydrophilic Surface Specificity on the Structural Properties of Confined Water[END_REF][START_REF] Argyris | Dynamic Behavior of Interfacial Water at the Silica Surface[END_REF][START_REF] Bonnaud | Molecular Simulation of Water Confined in Nanoporous Silica[END_REF][START_REF] Yamashita | Molecular Simulations of Water Adsorbed on Mesoporous Silica Thin Films[END_REF][START_REF] Yamashita | Grand Canonical Monte Carlo and Molecular Dynamics Simulations of Capillary Condensation and Evaporation of Water in Hydrophilic Mesopores[END_REF] It is well-known that the fluid adsorption is influenced by the chemical disorder of the adsorbing surface. [START_REF] Kuchta | Monte Carlo Simulations of Krypton Adsorption in Nanopores: Influence of Pore-Wall Heterogeneity on the Adsorption Mechanism[END_REF][START_REF] Puibasset | Capillary Condensation in a Geometrically and a Chemically Heterogeneous Pore: A Molecular Simulation Study[END_REF][START_REF] Puibasset | Influence of Surface Chemical Heterogeneities on Adsorption/Desorption Hysteresis and Coexistence Diagram of Metastable States within Cylindrical Pores[END_REF][START_REF] Kuchta | Microscopic Mechanism of Adsorption in Cylindrical Nanopores with Heterogeneous Wall Structure[END_REF][START_REF] Feng | Adsorption of Fluids in a Pore with Chemical Heterogeneities: The Cooperative Effect[END_REF] In particular, the structure and dynamics of an adsorbed water film has been shown to be influenced by the nature of the silica surface. [START_REF] Argyris | Dynamic Behavior of Interfacial Water at the Silica Surface[END_REF][START_REF] Puibasset | A Comparison of Water Adsorption on Ordered and Disordered Silica Substrates[END_REF] This is expected to have an impact on the heat capacity of adsorbed layers of water, to be studied now.

Numerical details

Model and interactions

The system to be studied is a water film adsorbed on a plane surface. The interactions between water molecules are described by the SPC model. [START_REF] Berendsen | Intermolecular Forces; B. Pullman[END_REF] This fast computable model is well suited as it reproduces quite well the thermodynamic properties of water around ambient temperature, like vapor pressure, enthalpy of vaporization and heat capacity. [START_REF] Jorgensen | Comparison of Simple Potential Functions for Simulating Liquid Water[END_REF] In this model, a Lennard-Jones center is located on the oxygen atom, and positive and negative charges are placed on the three atoms. The water-water contribution to energy is calculated directly with minimal image convention, with a cutoff for the LJ contribution situated at 0.95 nm. Long range corrections are not applied for the LJ term. However, long range electrostatic contributions are handled with the Ewald summation. The direct space parameter kappa = 2.288 nm -1 and 7 [START_REF] Zanotti | Competing Coexisting Phases in 2d Water[END_REF] reciprocal vectors are chosen. The total expected precision is around 1%.

The simulation box is cubic for the bulk and the films configurations. The number of molecules for the bulk simulations is 216, with dimensions L = 1.863 nm, while for the film configurations there are 220 molecules on a surface with dimensions L=3.72 nm.

The film configurations correspond to water adsorbed on a hydrophilic surface, typically silica. The interaction between water molecules and silica species contain both van der Waals and electrostatic contributions due to hydroxyls at the silica surface. However, in this preliminary study we will neglect the atomic roughness of the adsorbing surface (smooth wall approximation), as well as the presence of hydroxyl groups. As a result, the external potential felt by the adsorbed water molecules is that produced by a continuum of silica atoms in a half-space interacting with the water oxygen species via the 12-6 LJ potential. The integration results in the 9-3 potential

( ) [ ] (1)
where is the density of the solid interaction sites, and and are the parameters of the 12-6 LJ interaction between water molecules and silica. [START_REF] Puibasset | Water Adsorption on Hydrophilic Mesoporous and Plane Silica Substrates: A Grand Canonical Monte Carlo Simulation Study[END_REF] However, the external potential obtained with this procedure is far too weak to be representative of a hydrophilic surface, because of the ignored surface dipoles associated with hydroxyls which are responsible for a major contribution to the interaction. As a consequence, we have chosen to enhance the external potential so that the minimum equals -15 kJ/mol, representative of a hydrophilic silica surface.

Adsorption isotherm and heat capacity

The adsorption isotherm is calculated as follows. The initial configuration is the empty simulation box. We set the chemical potential at its lowest value, corresponding to an external vapor pressure of one tenth of the saturating pressure P sat (P/P sat = 0.1). The grand canonical Monte Carlo (GCMC) procedure is applied until the equilibrium is reached. After equilibration is reached, a new GCMC is started for a larger value of the chemical potential (corresponding to P/P sat = 0.2), the initial configuration of the new run being the final configuration of the previous one. The same procedure is repeated until the saturating chemical potential is reached. These GCMC calculations are devised to produce and equilibrate the water film used to calculate the heat capacity.

The heat capacity is calculated in the canonical ensemble, i.e. for a given film thickness, or, equivalently, for a fixed number of molecules. Note however that the system contains a liquid/vapor interface, which is free to move to equilibrate pressure and chemical potential between the liquid film and vapor.The heat capacity is obtained from the standard fluctuations formula of the intermolecular enthalpy during the Monte Carlo run. [START_REF] Jorgensen | Convergence of Monte Carlo Simulations of Liquid Water in the NPT Ensemble[END_REF] The quantum corrections are taken into account so that the simulation results can be compared to experiments. [START_REF] Owicki | Monte Carlo Calculations in the Isothermal-Isobaric Ensemble. 1. Liquid Water[END_REF] The simulations are started at the highest temperature, the final configuration of a given run being used as the initial configuration for a lower temperature simulation. This procedure is thought to be the most efficient to equilibrate the system at low temperature. Long Monte Carlo runs are however required to reach the equilibrium, in particular at low temperature where mobility of the molecules is extremely low. We have performed 10 6 Monte Carlo steps per molecule for each temperature for equilibration, followed by 5×10 [START_REF] Ouyang | Modelling Water: A Lifetime Enigma[END_REF] steps per molecule for data acquisition.

Results and discussion

Adsorption isotherm at 300K

The GCMC results for the adsorption isotherm are given in Figure 1. The amount of water adsorbed increases with the chemical potential, and exhibits a step around P/P sat = 0.3. This corresponds to the formation of the first layer of water. A second step is visible around P/P sat = 0.8, which is associated with the formation of a second layer of water molecules on the surface. This evolution of the amount adsorbed is in agreement with previous calculations done by taking into account explicitly the electrostatic contributions due to partial charges of the silica species, in particular the presence of hydroxyls. [START_REF] Puibasset | A Comparison of Water Adsorption on Ordered and Disordered Silica Substrates[END_REF] The extremely simple model used in the present study to describe the surface thus catches the main features of a hydrophilic surface.

The heat capacity of the water film is calculated for two different film thicknesses corresponding approximately to 22 mol/m 2 (monolayer, P/P sat = 0.5) and 44 mol/m 2 (two layers, P/P sat = 0.8). Visual inspection of the molecular configurations (see Figure 2) shows that the first water layer is highly structured, with OH parallel to the adsorbing surface. As expected, the second layer is much less structured.

Heat capacity

The heat capacity dependence with temperature is given in Figure 3, for the two water film thicknesses adsorbed on the plane surface. As can be seen, the two curves exhibit non monotonic variations with the temperature. The general trend is that the heat capacity is independent of the film thickness at low and high temperature (below 150 K and above 300 K), while strong differences appear in the intermediate region. Both curves show an increase of the heat capacity, with a maximum around 225 K. However, the peak for the thin film is more pronounced, and probably narrower than for the thick film. How do these results compare with the bulk and experiments? The bulk simulation data show an almost monotonic increase of the heat capacity with temperature, except for a small excess between 200 and 250 K. Taking into account the rather large error bars of the simulations (10-20%) this excess might be the signature of an underlying peak around 225 K. Note that it is well known that the SPC water model predicts an increase of the heat capacity when reducing the temperature in a much lesser extent than experiments. [START_REF] Jorgensen | Temperature Dependence of TIP3P, SPC, and TIP4P Water from NPT Monte Carlo Simulations: Seeking Temperatures of Maximum Density[END_REF] An interesting feature emerges when comparing the results for the bulk and the two film thicknesses: the narrower the film, the higher the peak.

This suggests a strong correlation between the appearance of this peak and the structure of the first layer of adsorbed water. It is furthermore emphasized that the smooth surface used to model silica enhances the structuration of the first layer. These results are thus expected to be altered for a more realistic atomistic silica surface.

Another feature to be noticed is that at low temperature, the heat capacity of water converges towards the same value, independent of the film thickness, the value being that of the bulk water. This value remains significantly higher than the expected value for ice: the SPC water model is frozen in a metastable state even for the lowest temperature. The simulations are unable to cross the nucleation barrier associated with the crystallization of water.

An interesting point would be to check for the robustness of these results with the model. It is thus planned to perform similar calculations with other water models (TIP4P, TIP4P/2005) and another silica model including an atomistic description and silanols.

Figure 1. Grand Canonical Monte Carlo results of SPC water adsorption isotherm (T = 300 K) on a smooth hydrophilic surface (eq. 1). The horizontal line at 22 mol/m 2 correspond to a single water layer (P/P sat = 0.5) and at 44 mol/m 2 to two water layers (P/P sat = 0.8).

Figure 2 .Figure 3 . 29 temperature

 2329 Figure 2. Molecular simulation configurations of SPC water at 300K adsorbed on the smooth


Roger Mayer

Isabelle Picard

Philippe Lawton

Philippe Grellier

Christine Barrault

Michel Monsigny

Joseph Schrével

Peptide derivatives specific for a Plasmodium falciparum proteinase inhibit the

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Joël Puibasset

Fluid adsorption in linear pores: a molecular simulation study of the influence of heterogeneities on the hysteresis loop and the distribution of metastable states

Fluid adsorption in linear pores: a molecular simulation study of the influence of heterogeneities on the hysteresis loop and the distribution of metastable states

Joël Puibasset

INTRODUCTION

Essentially all surfaces adsorb fluids, due to molecular interactions [START_REF] Cushman | The Physics of Fluids in Hierarchical Porous Media: Angstroms to Miles[END_REF][START_REF] Guéguen | Introduction to the physics of rock[END_REF][START_REF] Gregg | Adsorption, Surface Area and Porosimetry[END_REF][START_REF] Rouquerol | Adsorption by Powders and Porous Solids[END_REF]. These interactions may be weak (van der Waals or electrostatic) or much stronger. For strong interactions, chemical bonds may be broken in the adsorbed molecules and/or established between the adsorbent and the adsorbate: this is chemisorption. On the other hand, for weak interactions, the adsorption process occurs reversibly without breaking any molecule: this is physisorption.

The paper will focus on this situation. Despite weak, the fluid/substrate interactions induce observable modifications of the fluid properties in the vicinity of the surface. Stronger effects may be induced when fluids are confined in porous materials, because the local curvature associated to porosity enhances fluid/substrate interactions. As an example, one observes that porous materials condense gases below their saturating point. The measured adsorption isotherm (giving the amount of fluid adsorbed versus fluid pressure) thus provides a powerful experimental tool for porous materials characterization [START_REF] Gregg | Adsorption, Surface Area and Porosimetry[END_REF][START_REF] Rouquerol | Adsorption by Powders and Porous Solids[END_REF]. In practice, it is not a straightforward task to assess the pore geometry from the isotherms. For instance, the local surface chemistry may enhance or disfavor adsorption: a famous example is given by the hydroxyl groups OH which influence the hydrophilic property of silica [START_REF] Legrand | The surface properties of silicas[END_REF][START_REF] Puibasset | A comparison of water adsorption on ordered and disordered silica substrates[END_REF][START_REF] Pellenq | Molecular simulations of water in hydrophobic microporous solids[END_REF][START_REF] Puibasset | Grand Canonical Monte Carlo simulation study of water adsorption in silicalite at 300K[END_REF][START_REF] Puibasset | Water confined in mesoporous silica glasses: Influence of temperature on adsorption/desorption hysteresis loop and fluid structure[END_REF][START_REF] Bussai | On the diffusion of water in silicalite-1 : MD simulations using ab initio fitted potential and PFG NMR measurements[END_REF][START_REF] Findenegg | Freezing and melting of water confined in silica nanopores[END_REF][START_REF] Takei | Changes in density and surface tension of water in silica pores[END_REF][START_REF] Trzpit | The effect of local defects on water adsorption in silicalite-1 Zeolite: a joint experimental and molecular simulation study[END_REF][START_REF] Zhuravlev | The surface chemistry of amorphous silica. Zhuravlev model[END_REF][START_REF] Morishige | X-ray diffraction studies of freezing and melting of water confined in a mesoporous adsorbent (MCM-41)[END_REF][START_REF] Bellissent-Funel | Structural study of water confined in porous glass by neutron scattering[END_REF][START_REF] Brovchenko | Effect of pore size on the condensation/evaporation transition of confined water in equilibrium with saturated bulk water[END_REF]. Another important case for environmental purpose regards the adsorption of condensable gases like carbon dioxide or methane on porous carbon-based materials. See for example the recent work by Wilcox et al. [START_REF] Liu | Effects of surface heterogeneity on the adsorption of CO2 in microporous carbons[END_REF][START_REF] Liu | Molecular simulation of CO2 adsorption in micro-and mesoporous carbons with surface heterogeneity[END_REF][START_REF] Liu | Molecular simulation studies of CO2 adsorption by carbon model compounds for carbon capture and sequestration applications[END_REF][START_REF] Firouzi | Molecular simulation of carbon dioxide transport and storage in porous carbon-based materials[END_REF][START_REF] Firouzi | Slippage and viscosity predictions in carbon micropores and their influence on CO2 and CH4 transport[END_REF] on the influence of surface chemistry and pore network on adsorption.

An important feature of such adsorption isotherms, which increases the difficulty of their analysis, is their possible irreversibility: desorption do not coincide with adsorption (existence of a hysteresis, to be discussed later). Furthermore, the size of the reservoir used to perform adsorption experiments has to be chosen carefully, i.e. large enough to avoid spurious effects, a condition which is not so frequently fulfilled since "large" should be understood in terms of number of molecules between the system and the reservoir and not in terms of volumes [START_REF] Kierlik | Effect of the reservoir size on gas adsorption in inhomogeneous porous media[END_REF][START_REF] Puibasset | Influence of reservoir size on the adsorption path in an ideal pore[END_REF][START_REF] Nguyen | Monte Carlo simulation of the gas-phase volumetric adsorption system: effects of dosing volume size, incremental dosing amount, pore shape and size, and temperature[END_REF]. Data analyses also have to take into account corrections arising from an accurate reconsideration of physics of adsorption [START_REF] Coasne | Domain theory for capillary condensation hysteresis[END_REF][START_REF] Coasne | Grand canonical Monte Carlo simulation of argon adsorption at the surface of silica nanopores: effect of pore size, pore morphology, and surface roughness[END_REF][START_REF] Coasne | Adsorption, intrusion and freezing in porous silica: the view from the nanoscale[END_REF].

Focusing on heterogeneous surfaces, we show how surface chemistry and geometry influence the adsorption properties of a condensable fluid and hence its liquid/vapor coexistence diagram. In most cases, molecular simulation [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Nicholson | Computer Simulation and the Statistical Mechanics of Adsorption[END_REF], taking into account the fluid/substrate interactions at atomistic scale, has proven to be a powerful tool to understand adsorbed fluid properties in the vicinity of a surface. In particular for studying the effect of surface chemistry or corrugation [START_REF] Bock | Phase behavior of a simple fluid confined between chemically corrugated substrates[END_REF][START_REF] Kuchta | Microscopic mechanism of adsorption in cylindrical nanopores with heterogeneous wall structure[END_REF], local curvature, pore morphology [START_REF] Coasne | Molecular simulation of gas adsorption in realistic models of silica nanopores[END_REF][START_REF] Coasne | Adsorption of simple gases in MCM-41 materials: the role of surface roughness[END_REF], or any other nanometerscale property of the surface [START_REF] Gelb | Phase separation in confined systems[END_REF].

We show that, increasing the surface heterogeneities has drastic effects on coexistence diagram. Not only a shift of the critical point (i.e. quantitative changes), but also the appearance of many metastable states that introduce qualitative changes [START_REF] Evans | Capillary condensation and adsorption in cylindrical and slit -like pores[END_REF][START_REF] Morishige | Adsorption hysteresis and pore critical tempereture in a single cylindrical pore[END_REF][START_REF] Sonwane | Analysis of criticality and isotherm reversibility in regular mesoporous materials[END_REF][START_REF] Puibasset | Pseudocritical or hysteresis temperature versus pore size for simple fluids confined in cylindrical nanopores[END_REF][START_REF] Puibasset | Influence of surface chemical heterogeneities on adsorption/desorption hysteresis and coexistence diagram of metastable states within cylindrical pores[END_REF]. A multiscale approach has been used to explore this issue, showing the importance of the interdependence between the various heterogeneities for the adsorption and desorption mechanisms, the observed hysteresis, and the distribution of the underlying metastable states [START_REF] Kierlik | Adsorption hysteresis and capillary condensation in disordered porous solids: a density functional study[END_REF][START_REF] Kierlik | Capillary condensation in disordered porous materials : hysteresis versus equilibrium behavior[END_REF]. Such approach, including both an atomistic and a large-scale description (i.e. from interatomic forces to large scale domain network) is expected to offer a powerful tool for experimental data analysis.

2.

MOLECULAR MODELS

The main objective of this simulation work is to describe trends in adsorption due to the presence of heterogeneities, ranging from one single heterogeneity to a large number. The model needs to be simple to catch the main features of adsorption, avoiding unnecessary (and possibly misleading) details. This is why the interactions are chosen to be extremely simple, with however the realistic orders of magnitude. Similarly, the heterogeneities will be generated in the simplest possible way, as described in the following sections.

All (fluid-fluid and fluid-wall) interactions are described by the (12-6) Lennard-Jones interatomic potential cut at 2.5 σ (interactions are not calculated beyond that distance) and shifted (the potential reaches zero at cutoff). The adsorbate is chosen to be argon, a simple fluid commonly used for porous materials characterization. The parameters for argon are ε Ar- Ar /k = 119.8 K where k is Boltzmann's constant, and σ Ar-Ar = 0.3405 nm. The fluid-fluid parameters provide natural units to measure the physical quantities introduced in this paper (the unit for distances is σ Ar-Ar ; for energies it is ε Ar-Ar ; for temperatures it is ε Ar-Ar /k; for pressures it is ε Ar-Ar /σ Ar-Ar 3 ; etc). The reduced quantities will be denoted by an asterisk. The parameters for argon/wall interactions are ε Ar-W /k = 100 K and σ Ar-W = 0.333 nm. These parameters are typical for argon-oxygen interactions, corresponding to the main contribution for argon/silica interaction. The adsorbent (pore) is chosen to be a cylinder of 3 nm in diameter. The smooth wall approximation may be used: it consists in averaging the fluid-wall interaction over the atomic roughness (the potential is purely radial in this case). Differences in adsorption isotherms are negligible (small differences at low pressure only), but it allows to save computing time.

We focus on simple fluid adsorption in the frequently encountered silica-based nanoporous materials like MCM-41 or oxidized porous silica. These materials exhibit parallel and noninterconnected tubular pores with a high aspect ratio (few nanometers in diameter and few microns of length) [START_REF] Beck | A new family of mesoporous molecular sieves prepared with liquid cristal templates[END_REF][45][START_REF] Uhlir | Electrolytic shaping of germanium and silicon[END_REF]. Such materials are of greatest fundamental interest since they exhibit unexpected hysteresis shapes [START_REF] Coasne | Adsorption in noninterconnected pores open at one or at both ends: a reconsideration of the origin of the hysteresis phenomenon[END_REF][START_REF] Grosman | Nature of capillary condensation and evaporation processes in ordered porous materials[END_REF][START_REF] Grosman | Capillary condensation in porous materials. Hysteresis and interaction mechanism without pore blocking/percolation process[END_REF]. Due to fabrication procedure, one expects heterogeneities at nanometer scale along the pore axis [START_REF] Maddox | Characterization of MCM-41 using molecular simulation: heterogeneity effects[END_REF][START_REF] Edler | Small-angle neutron scattering studies on the mesoporous molecular sieve MCM-41[END_REF][START_REF] Berenguer-Murcia | About the exclusive mesoporous character of MCM-41[END_REF][START_REF] Fenelonov | Comparative textural study of highly ordered silicate and aluminosilicate mesoporous mesophase materials having different pore sizes[END_REF][START_REF] Sonwane | A model for the structure of MCM-41 incorporating surface roughness[END_REF]. These heterogeneities result from variations of local surface chemistry (the density of surface hydroxyls groups may vary for instance) or from variations in geometric characteristics (constrictions in pores). These heterogeneities essentially result in a modulation of the affinity of the adsorbed fluid (argon in our case) with the wall. We have in particular shown that pore size variations actually predominantly result in variations of fluid/wall interactions, with moderate "purely geometric" effects [START_REF] Puibasset | Capillary condensation in a geometrically and a chemically heterogeneous pore: a molecular simulation study[END_REF][START_REF] Puibasset | Phase coexistence in heterogeneous porous media : a new extension to Gibbs Ensemble monte carlo simulation method[END_REF][START_REF] Puibasset | Grand potential, Helmholtz free energy, and entropy calculation in heterogeneous cylindrical pores by the grand canonical Monte Carlo simulation method[END_REF][START_REF] Puibasset | Thermodynamic characterization of fluids confined in heterogeneous pores by Monte Carlo simulations in the grand canonical and the isobaric-isothermal ensembles[END_REF]. As a simplification, the corresponding heterogeneities are characterized by a simple parameter which correlates directly to the local intensity of the fluid/wall interaction. More precisely, this parameter is an extra factor applied to the fluid/wall interaction. Values less than 1 (down to 0.7 in our model) produce domains with lower affinity for the fluid. Values larger than 1 (up to 1.3 in our model) produce domains with higher affinity for the fluid. This parameter is varied along the pore axis to produce heterogeneities. Several cases will be considered: no heterogeneity, one single heterogeneity, two heterogeneities with different amplitudes, and a large number of heterogeneities with random amplitude taken from a given (Gaussian) distribution. The first case corresponds to a perfectly cylindrical pore, chosen to be of diameter 3 nm and length L =4.086 nm = 12 σ Ar-Ar (see Fig. 1-a). For the second case, the heterogeneity is produced by modulating the fluid-wall potential of the perfectly cylindrical pore by the extra factor

( ) { } L z π 2 cos 25 . 0 1 +
, where L is the pore length and z is the coordinate along the pore axis. This corresponds to a sinusoidal modulation along the pore axis of 25 percent maximum. The pore is thus a chain (due to periodic boundary conditions) of more or less attractive regions (see Fig. 1-b). In the third case, the modulation factor is

( ) { } L z z a π 4 sin ) ( 1+
where a = 0.3 for z<0 and a = 0.2 for z>0, and L = 8.172 nm = 24 σ Ar-Ar is the pore length. This produces four regions per period in the pore with high, medium high, medium low and low fluid/wall affinities (see Fig. 1-c). In the last case, many modulations are introduced in the pore. For simplicity a step function has been used: the modulation factor is constant within each domain (of length 7.36 nm = 21 σ Ar-Ar ),

and takes uncorrelated values between two adjacent domains. Its value is supposed to follow a Gaussian distribution (see Fig 1-d).

The simulations are performed using the grand canonical Monte Carlo (GCMC) algorithm [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Nicholson | Computer Simulation and the Statistical Mechanics of Adsorption[END_REF]. The chemical potential µ of argon is imposed and related to the gas pressure in the reservoir using the ideal gas relation since it is a good approximation at the temperature at which the simulations are performed (for instance This local minimum is the analog in heterogeneous cylindrical pores of the bridgelike phase observed in corrugated slit pores [START_REF] Bock | Phase behavior of a simple fluid confined between chemically corrugated substrates[END_REF]. See schematic representations in Fig. 2.

In the pore with two sinusoidal heterogeneities, five local minima are observed, which correspond respectively to: a gaslike fluid adsorbed in the pore; one liquidlike meniscus in the most attractive domain of the pore; two liquidlike menisci in the two attractive regions of the pore; one single large meniscus bridging the two attractive domains; a liquidlike fluid completely filling the pore (Fig. 2). This system shows very interesting features. As an example, for a reduced temperature of 0.8, starting from the empty system and gradually increasing the chemical potential results in continuous filling of the system up to µ = -10.84

where the system jumps into the second local minimum. Increasing further the chemical potential makes the system jump successively into the third, fourth and finally the fifth local minimum. Upon desorption, the system does not visit the local minima in the reverse order.

After the first step, which consists in the formation of a bubble in the less attractive domain (i.e. the system jumps from the fifth into the fourth local minimum), further decrease in chemical potential makes the system jump directly into the first local minimum (gaslike branch) without visiting the third and second minima. This is an important point which shows that it is possible to miss metastable states when performing simple adsorption or desorption cycles [START_REF] Sethna | Hysteresis and hierarchies : dynamics of disorder -driven first order phase transformations[END_REF].

This simple model shows that the hysteresis loop is related to a complex free energy landscape, and may present a dissymmetry in adsorption and desorption. It is also shown that the main adsorption and desorption isotherms are made of several branches with small discontinuities. Furthermore, a simple picture of adsorption/desorption subloops (scanning curves) arises, as irreversible paths between the various local minima.

PORE WITH A LARGE NUMBER OF HETEROGENEITIES: MULTISCALE APPROACH

Let us now consider pores with a large number of domains. Such systems are too large to be treated as a whole with usual molecular simulation tools. Let us introduce the multi-scale molecular simulation approach, which enables to study pores with a large number of heterogeneities. As introduced by Everett [60-62] a porous material can be seen as a collection of domains corresponding to elemental subdivision of the porous network (elemental cavities). Most of the time, the domains are considered to be independent to simplify the calculations (see Fig. 3). The global thermodynamic properties are thus obtained as an average over the domain distribution. Molecular simulation is able to treat separately each domain of nanometric scale, and thus provides a powerful analysis of fluid adsorption after averaging over the appropriate distribution of domains. This approach leads to an extremely accurate description of the low pressure region of the adsorption isotherm, where the hypothesis of independence between domains is well justified. But the hysteresis region is not well described by this approach because the disorder associated to the complex porous network plays an important role [START_REF] Kierlik | Capillary condensation in disordered porous materials : hysteresis versus equilibrium behavior[END_REF][START_REF] Mason | The effect of pore space connectivity on the hysteresis of capillary condensation in adsorption desorption isotherms[END_REF][START_REF] Mason | A model of adsorption -desorption hysteresis in which hysteresis is primarily developed by the interconnections in a network of pores[END_REF][START_REF] Swift | Phase transitions in a model porous medium[END_REF][START_REF] Kierlik | Equilibrium and out-of-equilibrium (hysteretic) behavior of fluids in disordered porous materials : theoretical predictions[END_REF][START_REF] Detcheverry | Local mean-field study of capillary condensation in silica aerogels[END_REF][START_REF] Detcheverry | Gas adsorption and desorption in silica aerogels: a theoretical study of scattering properties[END_REF][START_REF] Naumov | Understanding capillary condensation and hysteresis in porous silicon: Network effects within independent pores[END_REF].

As a matter of fact, in the hysteresis region, the amount of fluid adsorbed in any domain is not a single valued function of the chemical potential. It also depends on the fluid state in the neighboring domains (this is closely connected to the existence of multiple metastable states).

It is possible to approximately rationalize this dependence from a classical thermodynamic analysis of fluid adsorption in pores. However, significant improvement can be achieved if one uses molecular simulation to extract them, with appropriate explicit boundary conditions to take into account the presence of neighboring domains [START_REF] Puibasset | Adsorption/desorption hysteresis of simple fluids confined in realistic heterogeneous silica mesopores of micrometric length: a new analysis exploiting a multiscale Monte Carlo approach[END_REF][START_REF] Puibasset | Monte Carlo multiscale simulation study of argon adsorption/desorption hsyteresis in mesoporous heterogeneous tubular pores like MCM-41 or oxidized Porous silicon[END_REF]. This is possible because interactions and correlations are short range: as a consequence, only the (two) nearest neighboring domains need to be considered. Actually, only the portions of the neighboring domains within interaction range need to be considered for the calculations (the remaining part beyond the interaction range has no direct influence on the central domain). Furthermore, the main effect of the neighboring domains does not originate in their detailed physico-chemical properties: the relevant information is essentially the absence or presence of fluid in the neighboring domains, or, more quantitatively, the corresponding density. Therefore, these neighboring domains may be taken as homogeneous cylindrical pores with the same fluid/wall interactions as the central domain to be studied. The simulation box then contains the central domain (length: 7.36 nm = 21.6 σ Ar-Ar ) and two adjacent portions of domains (taken to be 3 σ Ar-Ar long) which play the role of boundaries; extra portions are added, where the fluid state (gaslike or liquidlike) is imposed during the calculations (see Fig. 4).

domain to be considered substrate adsorbed argon transition regions frozen region frozen region Figure 4: Schematic representation of the simulation box with explicit boundary conditions: the domain to be considered is in the center (green zone); the explicit boundaries are made of frozen regions where the fluid state (gaslike or liquidlike) is imposed and transition regions where the fluid relaxes between its imposed state in the frozen region and its actual state in the domain. Blue circles: adsorbed argon; red circles: substrate species. Note that it is not yet known if the neighboring domains are actually filled or empty: we just perform pre-calculations considering all possible hypotheses. In other words, these precalculations aim at determining how the considered domain behaves as a function of the neighboring domains. Essentially, the central domain is bistable (gaslike or liquidlike), and the imposed boundaries forces the system to switch to one state, depending on the chosen boundaries. It's only when all pre-calculations are performed that a real adsorption in the whole pore can be performed. The initial condition being the empty pore, we know that all domains have empty neighbors. The chemical potential is then slowly increased and the precalculated behavior of each domain is used to determine the filling status of each domain, taking into account the status of the neighboring domains. The total amount of fluid adsorbed versus µ is then deduced (adsorption isotherm). A similar procedure is used for desorption.

PORE WITH A LARGE NUMBER OF HETEROGENEITIES: SIMULATION

RESULTS

As for previous systems with few heterogeneities, the main adsorption and desorption isotherms have been calculated by slowly increasing the chemical potential from the initial empty pore up to saturation, and then decreasing the chemical potential until the pore empties completely. Scanning isotherms are obtained in a similar way. However, as seen previously, not all metastable states can be reached with this procedure [START_REF] Sethna | Hysteresis and hierarchies : dynamics of disorder -driven first order phase transformations[END_REF]. To be sure to determine all metastable states for our system, we proceed as follows. We first generate all possible states, characterized by a set of gaslike and liquidlike status for all domains. The total number of possible states is necessarily less than 2 N (where N is the number of domains) since all domains have at most two possible states. We then eliminate all states which are not compatible with the pre-calculated isotherms. From a physical point of view, these configurations correspond to situations where the fluid state is unstable and has to evolve towards a metastable state (local minimum of the free energy). When a new metastable state is found, its stability limits are easily determined by slowly decreasing and/or increasing the chemical potential. The corresponding branches can be drawn in the usual graph giving the amount adsorbed versus chemical potential.

Figure 5 gives the results for a pore containing 15 domains. In this particular system, a total of N = 135 metastable states where found. The lowest and highest branches correspond to the usual gaslike and liquidlike states, i.e. fluid adsorbed only at surface or completely filling the pore respectively. The other (metastable) states correspond to partial filling of the pore (some domains may be filled while others are still essentially empty. They appear as portions of lines corresponding to their metastability range. These lines can be seen as the projection in the µ-ρ plane of the local minima in the free energy landscape. As expected, the slopes are always positive due to increasing film thickness with µ. These lines are almost straight and essentially parallel, which is a remarkable feature. As a consequence, despite eventually great differences in fluid configuration between two adjacent metastable states, the response of the system to a chemical potential variation is essentially the same in both cases. It is also observed that the main hysteresis loop corresponds to the external envelope of the metastable states. This feature is observed in other systems and seems to be quite general [START_REF] Detcheverry | Metastable states and T = 0 hysteresis in the random-field Ising model on random graphs[END_REF][START_REF] Rosinberg | Numerical approach to metastable states in the zero-temperature random-field Ising model[END_REF]. These results are typical for a pore with a dozen domains, but details in the number and positions of the metastable states are dependent on the particular distribution of the heterogeneities in the pore. This effect is expected to be minimized for large pores containing a large number of heterogeneities, closer to experimental situations. As a consequence, an appropriate finite size analysis has been done to extrapolate to large systems. To do this analysis, and for each system size considered, an appropriate average over disorder (distribution of heterogeneities) has also to be performed to minimize spurious effects due to the choice of particular configurations for the heterogeneities. Previous works [START_REF] Detcheverry | Metastable states and T = 0 hysteresis in the random-field Ising model on random graphs[END_REF][START_REF] Rosinberg | Numerical approach to metastable states in the zero-temperature random-field Ising model[END_REF][START_REF] Cherrier | The number of metastable states in the generalized random orthogonal model[END_REF][START_REF] Illa | Influence of the driving mechanism on the response of systems with athermal dynamics: the example of the random-field Ising model[END_REF][START_REF] Puibasset | Counting metastable states within the adsorption/desorption hysteresis loop: a molecular simulation study of confinement in heterogeneous pores[END_REF] have

shown that the number of metastable states is expected to increase exponentially with system size (number of domains), and that the appropriate self-averaging quantity is the corresponding normalized logarithm. We thus introduce the quenched complexity as

N ln 1 lim N N ∞ → = Σ
where N is the number of metastable states for each system, N is the number of domains (measure of the system size), and the bar denotes the average over disorder. In order to get insight in the distribution of the metastable states within the hysteresis loop, the local quenched complexity is introduced: Figure 6 shows the results for N = 25. This system size allows a good compromise between large N and a large number of disorder realizations (5×10 4 ) to increase statistical accuracy. As can be seen, the complexity presents a typical shape with maximal values in the hysteresis loop. An interesting feature is that the location where ( ) 0 ; = Σ µ ρ closely follows the gaslike and liquidlike adsorption/desorption branches, within statistical accuracy and discretization: this result is in agreement with what has been observed in other systems [START_REF] Detcheverry | Metastable states and T = 0 hysteresis in the random-field Ising model on random graphs[END_REF][START_REF] Rosinberg | Numerical approach to metastable states in the zero-temperature random-field Ising model[END_REF][START_REF] Illa | Influence of the driving mechanism on the response of systems with athermal dynamics: the example of the random-field Ising model[END_REF].

6.

CONCLUSION

A molecular simulation approach has been used to model simple fluid adsorption in heterogeneous tubular pores mimicking mesoporous materials like MCM-41 or porous silicon (with typical argon/silica interaction parameters). A hysteresis loop is observed in adsorption / desorption cycles, which is closely connected to the appearance of many metastable states.

The density of these metastable states is studied in the µ-ρ plane. Although this quantity is not directly accessible to experiments, it is of primary theoretical importance to understand the location of the hysteresis and the shape of the scanning curves as determined experimentally.

As a conclusion, the main adsorption/desorption hysteresis loop appears as the envelope of the numerous metastable states appearing in the system. This work focuses essentially on one-dimensional systems. From this point of view, the results are complementary of other fundamental studies focusing for instance on the effect of surface and finite size on the occurrence of true phase transitions. [START_REF] Winkler | Capillary condensation in cylindrical pores: Monte Carlo study of the interplay of surface and finite size effects[END_REF] This fundamental work can also, on the other hand, be considered for more applicative studies. As has been suggested, it could be applied to strongly adsorbing systems like carbon dioxide or methane in carbonaceous materials for industrial applications (storage, gas shale, etc.).

Figure 1 :

 1 Figure1: Presentation of the four pore models. a) pore 0: fluid/wall potential in reduced units for the perfectly cylindrical pore (length 12σ Ar-Ar ); r* is the reduced radial coordinate and z* is the reduced coordinate along the pore axis; b) pore I: fluid/wall potential in reduced units for the pore with one single sinusoidal heterogeneity deduced from pore 0 by sinusoidal modulation along z* (length 12σ Ar-Ar ); c) pore II: fluid/wall potential in reduced units for the pore with two sinusoidal heterogeneities of different strength, giving rise to four domains with different fluid/wall affinities (length: 24σ Ar-Ar ); d) pore III (many heterogeneities): modulation factor along the domains (140 in this example, corresponding to a 1 micrometer long pore); the values are taken from the Gaussian distribution shown below.

Figure 2 :

 2 Figure 2: Adsorption/desorption isotherms for argon confined in the three systems with few heterogeneities obtained by Grand Canonical Monte Carlo simulation (symbols). The reduced temperatures are given in the inset. Lines are guide to the eye. The vertical lines materialize the system jumps when the metastability limits are exceeded. The large circles contain schematic representation of the fluid distribution in the pores (black: high density liquidlike fluid; white: low density gaslike fluid).

  to 70 K). The adsorption/desorption isotherms are given as a function of the reduced chemical potential µ * = µ / Ar Ar-ε . Long simulations are performed in order to reduce the uncertainties down to 7% (10 8 Monte Carlo trials per argon where required).

Figure 2

 2 Figure2shows the adsorption/desorption isotherms for fluid adsorbed in a cylindrical pore

Figure 3 :

 3 Figure 3: Schematic representation of a heterogeneous pore (in the center; four different domains labeled A to D can be seen). The system may either be treated as if the domains where independent (left, all domains are treated independently with periodic boundary conditions, PBC) or as a whole (right, the interdependence between domains is taken into account by means of explicit boundary conditions, see text).

Figure 5 :

 5 Figure 5: Monte Carlo simulation result of fluid adsorption in a heterogeneous pore comprising 15 domains, taking into account the interdependence between domains. Black continuous lines: main adsorption/desorption isotherms. Red portions of line: all (135) possible metastable states for this system.

Figure 6 :

 6 Figure 6: Local quenched complexity (see Eq 1) obtained by Monte Carlo simulation result of fluid adsorption in a heterogeneous pore comprising 15 domains, taking into account the interdependence between domains. The roughness of the periphery is due to the discretization used for numerical calculations. Solid lines: main adsorption/desorption isotherms.

ACKNOWLEDGMENTS

Dr. E. Kierlik and Dr. G. Tarjus are gratefully acknowledged for helpful discussions. The simulations were performed thanks to grant from the Institut de Développement des Ressources en Informatique Scientifique IDRIS-CNRS, Orsay, France.