Sylvain Conchon
email: sylvain.conchon@lri.fr

Giorgio Delzanno
email: giorgio.delzanno@unige.it

Angelo Ferrando
email: angelo.ferrando@unige.it

Parameterized Verification of Topology-sensitive Distributed Protocols goes Declarative

, an SMT-based infinite-state model checker, can be applied as a verification engine for GLog, a logic-based specification language for topology-sensitive distributed protocols with asynchronous communication. Existential coverability queries in GLog can be translated into verification judgements in Cubicle by encoding relational updates rules as unbounded array transitions. We apply the resulting framework to automatically verify a distributed version of the Dining Philosopher mutual exclusion protocol formulated for an arbitrary number of nodes and communication buffers.

Introduction

Automated verification of distributed systems is a difficult task for standard model checkers [START_REF] Bloem | Decidability of Parameterized Verification[END_REF][START_REF] Bloem | Decidability in parameterized verification[END_REF]. Protocols designed to operate in distributed systems are often defined for an arbitrary number of nodes, arbitrary connection topology, and asynchronous communication. protocol rules typically depend on the current network configuration (e.g, presence of a communication link, state of all connections, etc). Several formal languages have been proposed to specify this class of systems, e.g., communicating state machines, automata, process algebraic languages, (graph) rewriting, etc. In this setting safety properties can be nicely formulated by lifting decision problems based on reachability and coverability, in which the initial configuration is typically fixed a priory, to formulations that are existentially quantified over an infinite set of initial configurations. Existentially quantified coverability problems have been considered in [START_REF] Delzanno | Parameterized verification of ad hoc networks[END_REF][START_REF] Delzanno | On the power of cliques in the parameterized verification of ad hoc networks[END_REF][START_REF] Delzanno | Verification of ad hoc networks with node and communication failures[END_REF][START_REF] Bertrand | Distributed local strategies in broadcast networks[END_REF][START_REF] Bertrand | On the decidability status of reachability and coverability in graph transformation systems[END_REF] in order to reason on parameterized formulation of distributed protocols with broadcast communication. The coverability decision problem [START_REF] Abdulla | Parameterized verification[END_REF] is typically used to formulate reachability of bad configurations independently from the number of components of a system. Therefore, a constructive way to solve an existentially quantified coverability problem for a formal specification of a distributed algorithm provides a characterization of initial configurations from which it is possible to reach a bad configuration (e.g. an anomaly in the protocol). Existentially quantified coverability problems turn out to be undecidable even for systems with a static communication topology and very basic interaction primitives like atomic broadcast communication [START_REF] Abdulla | Parameterized verification[END_REF][START_REF] Delzanno | A unified view of parameterized verification of abstract models of broadcast communication[END_REF][START_REF] Delzanno | Parameterized verification of ad hoc networks[END_REF][START_REF] Delzanno | On the power of cliques in the parameterized verification of ad hoc networks[END_REF][START_REF] Bloem | Decidability of Parameterized Verification[END_REF]. As mentioned before, communicating state machine (automata), adopted in [START_REF] Delzanno | Parameterized verification of ad hoc networks[END_REF][START_REF] Delzanno | On the power of cliques in the parameterized verification of ad hoc networks[END_REF][START_REF] Delzanno | Verification of ad hoc networks with node and communication failures[END_REF], and graph rewriting, adopted in [START_REF] Bertrand | On the decidability status of reachability and coverability in graph transformation systems[END_REF], are two examples of formal description languages for this kind of systems. communicating state machines can be considered a standard design and protocol specification language adopted in several verification tools like Uppaal and Spin. Graph rewriting systems are well-suited for representing topology-sensitive rules as shown by the examples of the Groove tool suite. In several case-studies, protocol rules require complex guards that require the inspection of the state of nodes, links, paths, vicinity etc. and tables to store information collected during the execution of protocol phases. The combination of these features seems to require more general specification formalisms. To this aim, in [START_REF] Delzanno | A logic-based approach to verify distributed protocols[END_REF] we proposed to adopt a logic-based declarative language, named GLog, a fragment of both DCDCs [START_REF] Delzanno | A logic-based approach to verify distributed protocols[END_REF] and MSR [22]. GLog can be viewed as a logic-based presentation of graph update rules with global conditions expressed using quantified first order formulas. GLog is based on a quantified predicate logic in a finite relational signature with no function symbols. Configurations are represented here as sets of ground atomic formulas (instances of unary and binary predicates). Update rules consist of a guard and two sets of first order predicates that define resp. deletion and addition of state components. Differently from specification languages based on extension of Petri nets like transfer and broadcast protocols, guards are checked atomically but update transitions have only local effect. In other words, we forbid simultaneously update of the state of all nodes in a graph. Update rules can be applied to update a global configuration node by node and to operate on the vicinity of a node by restricting updates to given predicates. Termination of an update subprotocol can then be checked via a global condition. Similar specification patterns have been applied to model non-atomic consistency protocol and mutual exclusion protocols with non-atomic global conditions. GLog has been applied to manually analyze distributed protocols in [START_REF] Delzanno | A logic-based approach to verify distributed protocols[END_REF]. In the present paper we show that Cubicle [START_REF] Conchon | Cubicle: A parallel smtbased model checker for parameterized systems -tool paper[END_REF][START_REF] Mebsout | Inférence d'invariants pour le model checking de systèmes paramétrés[END_REF], an SMT-based infinite-state model checker based on previous work by Ghilardi et al. [START_REF] Alberti | A framework for the verification of parameterized infinite-state systems[END_REF], can be applied as automated verification engine for existentially quantified coverability queries in GLog. In Cubicle parameterized systems can be specified as unbounded arrays in which individual components can be referred to via an array index. The Cubicle verification engine performs a symbolic backward reachability analysis using an SMT solver for computing intermediate steps (preimage computation, entailment and termination test) and applies overapproximates predecessors using upward closed sets as in monotone abstractions [START_REF] Abdulla | Monotonic abstraction: on efficient verification of parameterized systems[END_REF]. A peculiar feature of Cubicle w.r.t. MCMT [START_REF] Alberti | A framework for the verification of parameterized infinite-state systems[END_REF] is that the tool can handle unbounded matrices. This is particularly relevant when modeling topology-sensitive protocols as done in GLog using binary relations defined over component identifiers. Furthermore, existentially quantified coverability decision problems in GLog can directly be mapped into Cubicle. More specifically, the encoding transforms GLog update rules into array-based update formulas in Cubicle. Classes of initial configurations are specified by using partial specifications of initial configurations in Cubicle verification judgements. Infinite sets of bad configurations can be expressed using unsafe configurations in Cubicle verification judgements.

As a case-study, we consider the distributed version of the Dining Philosopher mutual exclusion protocol (DDP) recently studied in [START_REF] Namjoshi | Analysis of dynamic process networks[END_REF][START_REF] Delzanno | Parameterized verification of graph transformation systems with whole neighbourhood operations[END_REF][START_REF] Delzanno | A logic-based approach to verify distributed protocols[END_REF]. The protocol deals with an arbitrary, finite number of nodes and buffers that act as singleplace communication channels, and arbitrary link topology between nodes and buffers. Ownership of buffers is specified using asynchronous rules. Global conditions over linked buffers are used as enabling conditions for acquiring access to resources shared among neighbors. The GLog formal specification of DDP is mapped to a Cubicle verification problem in a natural way. In our preliminary experiments, Cubicle verified the correctness of the protocols in negligible time. Furthermore, as expected it reports potential error traces when introducing network reconfiguration rules (e.g. dynamic link creation and deletion) unrelated to the state of the corresponding involved nodes. The application of declarative specification languages and SMT-based engine seems a very promising research line for dealing for a larger class of distributed algorithms.

Contents

The paper is organized as follows: we first present the GLog declarative language and, in that context, the existential coverabiliity problem; we then introduce Cubicle and exhibit a general encoding of existential coverability into Cubicle; we then discuss the experimental evaluation on the DDP case-study, and, finally, discuss other examples and future directions.

GLog

GLog [START_REF] Delzanno | A logic-based approach to verify distributed protocols[END_REF] formulas are based on a simple relational calculus that can be used to express updates of sets of ground atoms. A set of ground atoms can be interpreted as the current state or configuration of the system we are modeling. Update rules contain a formula working as a condition and deletion and addition sets that specify ground atoms to be deleted and added to the current state. More formally, let P be a finite set of names of (unary and binary) predicate names, N a denumerable set of node identifiers equipped with a total order <, V be a denumerable set of variables. Predicates in P are used to model current configurations. In addition to predicates in P , we interpret the binary relation lt as the total order < in our model. Our logic has no function symbols but can be instantiated with elements from N . An atomic formula is either a formula p(x), lt(x, y) or p(x, y), where p ∈ P , x, y ∈ V ∪ N A ground atom is a either a p(n), lt(n, m), or p(n, m), where n, m ∈ N . A literal is either an atomic formula or the negation ¬A of an atomic formula A. A formula is a first order formula built on literals, namely, any literal is a formula, conjunctions, disjunctions, universally and existentially quantified formulas are still formulas. Multiple occurrences of the same variable implicitly model equality constraints. The set of free variables of a formula F , namely F V (F ), is the minimal set satisfying

F V (p(x, y)) = {x, y}, F V (A ∨ B) = F V (A) ∪ F V (B), F V (A ∧ B) = F V (A) ∩ F V (B), F V (¬A) = F V (A), F V (∀v.A) = F V (A) \ {v}, and F V (∃v.A) = F V (A) \ {v}. Given S = {F 1 , . . . , F n }, we define F V (S) = F V (F 1 ) ∪ . . . ∪ F V (F n ).
Quantified formulas we will be used as application conditions of rules.

Configurations, Interpretations and Update Rules As mentioned before a set of ground atoms will be used to model a configuration. Formally, a configuration is a finite set ∆ of ground atomic formulas with predicates in P . A configuration implicitly defines a graph in which directed edges are represented by atomic formulas whose predicate name acts as edge label. Configurations can also be viewed as models in which to evaluate a conditions. An interpretation is a mapping σ from V to N . We use here a fixed interpretation of variables. The interpretation domain however consists of a denumerable set of node identifiers. For a formula F we use F σ as an abbreviation for σ(F ), where σ is the natural extension of σ to terms. For a set S = {A 1 , . . . , A n }, we use Sσ to denote the set {A 1 σ, . . . , A n σ}. Update rules consists of conditions defined by quantified formulas with no function symbols, a deletion and an addition set. The deletion (resp. addition) set defines the set of ground atoms that have to be cancelled from (resp. added to) the current configuration. A rule has the following form C, D, A , where C is a quantified formula, D and A are two sets of atomic formulas with variables in V and predicates in P , and such that

F V (A) ∪ F V (D) ⊆ F V (C). A protocol P is defined as a set of rules.
Operational Semantics To fix an operational semantics for our language we need a support for the interpretation of relations and variables. We use ∆ |= A to define the satisfiability relation of a quantified formula A s.t. F V (A) = ∅. Let A[n/X] denote the formula obtained by replacing each free occurrence of X with n. The relation is defined by induction as follows.

∆ |= p(n), if p(n) ∈ ∆, ∆ |= lt(n, m), if n < m, ∆ |= p(n, m) for p ∈ P , if p(n, m) ∈ ∆, ∆ |= A ∧ B, if ∆ |= A and ∆ |= B, ∆ |= ¬A, if ∆ |= A, ∆ |= ∀X.A, if ∆ |= A[n/X] for each n ∈ N , and ∆ |= ∃X.A, if ∆ |= A[n/X]
for some n ∈ N . Given a configuration ∆, we say that the quantified formula A is satisfied in ∆, if there exists an interpretation σ s.t. Aσ is satisfiable. In order to apply a rule C, D, A to ∆, there must be an interpretation σ that satisfies the quantified formula C. The same interpretation σ is then applied to the atomic formulas in D and A. The resulting sets of atoms, say D and A respectively, are deleted from and added to ∆, respectively.

The operational semantics of a protocol P is given by a transition system T P = C, → , where C is the set of possible configurations, i.e., finite subsets of ground atoms with predicates in P , and →⊆ C ×C is a relation defined as follows. For ∆, ∆ ∈ C and a rule C, D, A ∈ P, ∆ → ∆ if there exists σ s.t. ∆ |= Cσ and

∆ = (∆ \ Dσ) ∪ Aσ. A computation is a sequence of configurations ∆ 0 ∆ 1 . . . s.t. ∆ i → ∆ i+1 for i ≥ 0.
We use → * to denote the reflexive and transitive closure of →. In a single step of the operational semantics a rule is evaluated in the current configuration by taking a sort of closed-world assumption, i.e., ground atomic formulas that do not occur in a configuration are evaluated to false. Furthermore, ground atomic formulas that are not deleted are transferred from the current to the successor configuration. The latter property can be viewed then as a sort of frame axiom. It is important to notice that, in general, a configuration ∆ has several possible successors. Indeed, depending of the chosen interpretation of free variables the same rule can be applied to different subsets of ground atoms contained in the same configuration. Furthermore, the choice of the rules to be applied at a given step is non-deterministic.

As an example, we consider possible application of GLog to the specification of distributed protocols. The key ingredient of the specification language is the combination of complex conditions and update rules to reason on graphs in which predicates can be viewed as labels of links between agents and communication buffers. We have shown that we can also add labels to individual agents and buffers, e.g., to represent their current state. Update rules can be used to dynamically reconfigure the graph, i.e., change labels, topology and add or delete agents. The separation between agents and buffers is convenient to model asynchronous communication. For instance, let us consider a protocol in which two agents need to establish a connection via a shared buffer.

-An agent n 1 of type A connects to a buffer e 1 in idle state (the buffer is free) and sets the state of the buffer to ready. -An agent n 2 of type B connects to e 1 in state ready and changes the state to ack. -Agent n 1 sends message m by changing the state of e 1 to msg m .

-Agent n 2 receives message m and updates the state of the channel to ack for further communications.

The protocol can be specified as follows. We use unary predicates to associate states to edges. send messages are non-deterministically generated. An initial configuration has the form idle(b 1 ), . . . , idle(b k ), where b i < b j for i = j, i, j : 1, . . . , k. For the sake of simplicity, we do not model the state of agents but only their capabilities (req, rec, send).

R C D A 1 idle(B) ∧ ¬req(A, B) {idle(B)} {ready(B), req(A, B)} 2 ready(B) ∧ ¬rec(A, B) {ready(B)} {ack(B), rec(A, B)} 3 true {} {send(A, B, M )} 4 ack(B) ∧ send(A, B) {ack(B), send(A, B)} {msg(B, M )} 5 msg(B, M ) ∧ rec(A, B) {msg(B, M ), rec(A, B)} {idle(B)}
In rule 1 a buffer B is locked by a non-deterministically generated request req(A, B) from sender agent A (a variable). In rule 2 a buffer B is locked by a non-deterministically generated request rec(A, B) from receiver agent A (a variable). Rule 3 nondeterministically generates a send action from agent A. Rule 4 synchronizes a send action from agent A with a buffer locked by the same agent. The (non deterministically generated) message M is stored in the buffer. Rule 5 synchronizes and consumes a message in the buffer with the receiver agent, releasing the buffer.

The model provides other form of interactions. For instance, we can model ordered buffers by forming lists of messages attached to a given edge as in the representation of the tape of the Turing machine.

We can also model synchronous communication as in the following example

C link(A, B) ∧ s 1 (A) ∧ link(E, B) ∧ s 2 (E) D {s 1 (A), link(A, B), link(E, B), s 2 (E)} A {link(A, B), s 1 (A), link(E, B), s 2 (E)}
Here s(A) and s (A) denote agent A resp. in state s and s , s 1 (E) and s 1 (E) denote agent E resp. in state s 1 and s 1 , and link(A, B) and link(Y, B) denote links to a common buffer B.

Existential Coverability

We consider here decision problems that generalize the standard notion of reachability between configurations. The key point is to reason about an infinite set of initial configurations in order to prove properties for protocol instances with an arbitrary number of nodes. For a set S of configurations, we first define the P ost and P re operators as follows P ost(S) = {∆ | ∃∆ ∈ S, ∆ → ∆ } and P re(S) = {∆ | ∃∆ ∈ S, ∆ → ∆}. We use P ost * (S) (resp. P re * (S)) to denote the reflexive-transitive closure of P ost (resp. P re).

We now introduce the ∃-coverability problem as follows.

Definition 1 ( ∃-coverability ). Given a protocol P, a set of target configurations T and a possibly infinite set of initial configurations I, ∃-coverability is satisfied for P, I and T , written ∃Reach(P, I, T ), if there exists ∆ ∈ T and a configuration ∆ 1 s.t. ∆ 1 ∈ P ost * (I) and ∆ ⊆ ∆ 1 .

By expanding the definition of P ost * , ∃Reach(P, I, T ) holds if there exists a configuration ∆ 0 ∈ I s.t. ∆ 0 → * ∆ 1 and ∆ ⊆ ∆ 1 for some ∆ ∈ T . The target T can be interpreted as a pattern to match or avoid in computations starting from initial configurations. If the set I consists of configurations consisting of an arbitrary, finite number of components than ∃-coverability formally describes a parameterized verification decision problem for specifications given in GLog.

The ∃-coverability problem turns out to be undecidable [START_REF] Delzanno | A logic-based approach to verify distributed protocols[END_REF].

From GLog to Cubicle

Cubicle is a model checker that can be applied to verify safety properties of arraybased systems, a syntactically restricted class of parametrized transition systems with states represented as arrays indexed by an arbitrary number of processes [START_REF] Conchon | Cubicle: A parallel smtbased model checker for parameterized systems -tool paper[END_REF][START_REF] Mebsout | Inférence d'invariants pour le model checking de systèmes paramétrés[END_REF] The transition req(i) is enabled if there exists index i (a process) such that Want[i]=false. Its effect is to set Want to true for index i and leave the array unchanged in all other positions. A system execution is defined by an infinite loop that at each iteration: (1) non-deterministically chooses a transition instance whose guard is true in the current state; and (2) updates state variables according to the action of the fired transition instance. Infinite sets of unsafe states (bad configurations) are defined by using unsaf e constraints. For instance, the judgement

unsafe (x y) { Crit[x] = True && Crit[y] = True }
specifies the infinite set of arrays Crit (with any size) in which there there exist two cells with value T rue.

The Cubicle verification engine is based on symbolic backward exploration. Cubicle operates over sets of existentially quantified formulas called cubes. Formulas containing universally quantified formulas (generated during the computation of predecessors) are over-approximated by existentially quantified formulas. The class of formulas manipulated by the backward reachability loop of Cubicle in not closed by pre-image in presence of universally quantified guards. To handle such formulas, Cubicle implements a safe but over-approximate pre-image computation. Given a cube ∃ ī.Φ and a guard G of the form ∀ j.Ψ ( j), the pre-image replaces G by the conjunction σ∈Σ( j, ī) Ψ ( j)σ of instances over the permutation of Σ( j, ī). In other words, in order to handle universally quantified guards, Cubicle applies monotone abstraction [START_REF] Abdulla | Monotonic abstraction: on efficient verification of parameterized systems[END_REF] and over-approximates predecessors via upward-closed sets of configurations. The search procedure maintains a set V and a priority queue Q resp. of visited and unvisited cubes. Initially, let V be empty and let Q contain the cubes representing bad states. At each iteration, the procedure selects the highest-priority cube Φ from Q and checks for intersection with the formula denoting the initial configurations (satisfiability of conjunction of Φ and formulas in the initial conditions). If the test fails, it terminates reporting a possibile error trace. If the test passes, the procedure proceeds to the subsumption check, i.e., implication between formulas. If subsumption fails, then add Φ to V , compute all cubes in pred t (for every t), add them to Q, and move on to the next iteration. If the subsumption check succeeds, then drop Φ from consideration and move on. The algorithm terminates when a safety check fails or Q becomes empty. When an unsafe cube is found, Cubicle actually produces a counterexample trace. Safety checks, being ground satisfiability queries, are easy for SMT solvers. The challenge is in the subsumption check because of their size and the existential implies existential logical form. Cubicle applies the heuristics described in [START_REF] Conchon | Cubicle: A parallel smtbased model checker for parameterized systems -tool paper[END_REF] to handle subsumption. The BRAB algorithm introduce in [START_REF] Conchon | Invariants for finite instances and beyond[END_REF] automatically computes over-approximations of backward reachable states that are checked to be unreachable in a finite instance of the system (using Murφ). The resulting approximations (candidate invariants) are model checked together with the original safety properties. Completeness of the approach is ensured by a mechanism for backtracking on spurious traces introduced by too coarse approximations.

Encoding GLog in Cubicle

In this section we present an encoding of GLog into an array-based specification language. The encoding is quite natural. The interpretation domain of variables is that of process indexes. For each unary predicate p ∈ P , we introduce a corresponding Boolean array variable array p.

array p[proc] : bool

For each binary predicate q, we introduce a two-dimensional Boolean array q array q[proc,proc] : bool Encoding of guards is straightforward. Free variables occurring in GLog update rules become parameters of transition definitions. A literal q(x, y) [resp. q(x)] is mapped to the formula q(x, y) = true [resp. q(x) = true]. A literal ¬q(x, y) [resp. ¬q(x)] is mapped to the formula q(x, y) = f alse [resp. q(x) = f alse]. Compound/quantified require conditions are mapped to compound/quantified formulas over literals.

Some care has to be taken in the encoding of GLog update rules. Transitions in Cubicle operate simultaneously on every cell of an array to provide support for global operations like reset and transfer. This kind of operations are not provided in GLog since the focus is on asynchronous behavior, i.e., we assume that global operations are split into several asynchronous operations equipped with guards that can be used to check for the current state of the protocol phase under consideration.

To encode a deletion rule, operating on the atomic formula A(x, y), we use auxiliary variables u, t and a case analysis on indexes: for the case x = u, y = t we add the action A[x, y] := f alse, and A[x, y] := A[x, y] in all other cases. To encode an addition rule, operating on the atomic formula A(x, y), we use again auxiliary variables u, t and a case analysis on indexes: for x = u, y = t we add the action A[x, y] := true, and A[x, y] := A[x, y] in all other cases.

To encode ∃-coverability, we also need to specify initial and unsafe configurations. Unsafe configurations can be described as in Cubicle using an existentially quantified formula over array cells. To select classes of initial states, we can use init declarations in which we specify only partial conditions on array cells.

Case Study: Distributed Dining Philosophers

We consider here a distributed version of the dining philosopher mutual exclusion problem presented in [START_REF] Namjoshi | Uncovering symmetries in irregular process networks[END_REF]. Agents are distributed on an arbitrary graph and communicate asynchronously via point-to-point channels. Channels are viewed as buffers with state. Distributed Dining Philosophers (DDP) is defined as follows. The goal is to ensure that agents can access a resource shared in common with their neighbors in mutual exclusion. The protocol from the perspective a single agent consists of the following steps:

-Initially, all agents are in idle state.

-When an agent A wants to get a resource, A has to acquire the control of each buffer shared with his/her neighbors. -To acquire a channel, A marks the channel with its identifier. If the channel is already marked, A has to wait. -A acquires the resources when all channels shared with neighbors are marked with his/her identifier. -To release a resource, A first resets each buffer. When all buffers are reset, A moves back to idle state.

In a statically defined topology, agent A gets access to a resource when all neighbors are either idle or are waiting for acquiring some channel. Communication between two neighbors is asynchronous. Indeed, they interact by reading and writing on the shared channel. The protocol should guarantee that two agents that share the same channel cannot acquire and use a resource simultaneosly. The protocol should be robust under dynamic reconfigurations of the network.

Formal Specification of DDP

In this section we present a formal specification of the DDP protocol. Network configurations are expressed as GLog configurations. The dynamics in a protocol interaction is expressed via a finite set of update rules. We use a predicate link to represent connections from an agent to a possibly shared buffer. We model buffers with states using unary predicates. Asynchronous communication is modeled as in the previous example, i.e., agents interact only via a common buffer. Communication between two agents is not atomic. Instead of modeling identifiers and buffers with data, we introduce a special relation own that is used to model ownership of a given buffer to which a agent is linked. Ownership is normed in the same way as the labeling of buffers in the original protocol, i.e., an agent can acquire ownership only if the buffer is not owned by other agents.

Ownership can be realeased when in idle state. We also model non-deterministic creation (in idle state) and deletion of links. We model this behavior using the following predicates and rules (rules have the form (C i , D i , A i ) for i : 1, . . . , 6): R C D A getE link(X, E) ∧ ∀Z.¬own(Z, E) ∅ {own(X, E)} relE {idle(X), own(X, E)} {own(X, E)} ∅ acquire idle(X) ∧ ∀E.(link(X, E) ⊃ own(X, E)) {idle(X)} {busy(X)} release {busy(X)} {busy(X)} {idle(X)}

An initial state configuration has the following form idle(n 1 ), . . . , idle(n k ), where n i = n j for i = j, i, j : 1, . . . , k and k ≥ 1.

Encoding in Cubicle

Following the encoding rules specified for GLog, we can now obtain a Cubicle specification for DDP. To simplify a little bit the specification, we introduce an enumeration type state over node states.

type state = Idle | Busy

This way, we can use a single array State with cell type state instead of two Boolean arrays. We also need a link array and a own array to specify link and ownership relations between nodes and buffers.

array State[proc] : state array Link[proc,proc] : bool array Own[proc,proc] : bool

The initial configuration consists of all possible topologies in which nodes are in idle state. We also enforce the ownership relation to be false for each pair of node and buffer. Two node are in mutex state while pointing to at the same buffer. The transitions are obtained via the encoding of the GLog specification into Cubicle input language shown in Appendix A. When applying Cubicle to the above described problem, the tool proves the model correct in few seconds without need to apply multicore optimizations via the Functory library. More specifically, Cubicle visits 19 nodes with at most 3 process indexes, 529 fixpoint checks, and 176 calls to the Alter-Ego SMT solver. Since Cubicle operates over unbounded arrays, the above result provides a formal correctness proof of the considered model for any number of nodes and links and any topology. The proof certificate can be obtained by taking the set of assertions (formulas) collected during the fixpoint computation. By applying the BRAB algorithm with parameter 3, the number of visited nodes reduces to 10 with 312 fixpoint tests, and 88 calls to the Alter-Ego solver. Furthermore, the BRAB algorithm infers the invariant ¬(Own[#1, #3] = T rue ∧ Own[#2, #3] = T rue).

Dynamic Reconfiguration

To model dynamic reconfigurations, we can non-deterministically add and remove link predicates between pairs of agents and buffers. We first consider the non-deterministic rules link and unlink defined below.

R C D A Link ¬link(X, E) ∅ {link(X, E)} U nlink link(X, E) {link(X, E)} ∅
When the model extended with the above rules is checked with Cubicle (see appendix B), the tool reports the error trace acquire(#1) → link(#1, #3) → get(#2, #3) → acquire(#2) → unsaf e [START_REF] Abdulla | Parameterized verification[END_REF]. This trace is a real error trace. Indeed, process p 1 can acquire ownership when there are no links to buffer b 3 . Since the link rule has no condition on p 1 , a link can then be added from p 1 to b 3 . However, Process p 2 can now become owner of b 3 and then move to state Busy. Two processes are linked to the same buffer b 3 while in state Busy. We can modify the model and restrict addition of new link connected to node X only when X is in state Idle as follows.

Conclusions

We have studied a possible application of SMT-based infinite-state model checker to the verification of topology-sensitive distributed protocols, i.e., protocols defined over network graphs and in which rules have guards and effects that depend on communicaton links. Starting from a logic-based presentation of distributed protocols based on the GLog relational update language, we have shown how to encode existential coverability queries in GLog as Cubicle verification judgements. As a case-study, we have shown that the declarative approach supported by GLog+Cubicle provides a very effective way to verify protocols operating on graphs. For instance, in previous work DDP required complex verification methodologies like assume-guarantee reasoning or ad hoc algorithms for graph rewriting systems. In the present paper DDP is verified using a very simple declarative specification and a general purpose model checker. Cubicle verifies correctness in negligible execution time (without need of multicore optimizations via the Functory library).

The proposed methodology can be applied to other types of distributed protocols. More specifically, we are currently studying how to deal with routing protocols for arbitrary topologies and hierarchical protocols for reference counting (e.g. garbage collections disciplines etc). Another interesting direction is related to the possible application of Cubicle for verification of protocol specifications in parametereized multi-agent systems [START_REF] Ancona | Parametric runtime verification of multiagent systems[END_REF].

var

  Turn : proc array Want[proc] : bool array Crit[proc] : bool init (z) { Want[z] = False && Crit[z] = False } unsafe (x y) { Crit[x] = True && Crit[y] = True } The system state is defined by a set of global variables and arrays. The initial states are defined by a universal conjunction of literals characterizing the values for some variables and array entries. init (z) { Want[z] = False && Crit[z] = False } A state of our example consists of a process identifier Turn and two boolean arrays Want and Crit; a state is initial iff every cell of both arrays are set to false Transitions are given in the usual guard/action form and may be parameterized by (one or more) process identifiers. Guards are expressed via required expressions. They are quantified formulas. Quantification is defined only over variables of type proc. As an example, consider the following rule. transition req (i) requires { Want[i] = False } { Want[j] := case | i = j : True | _: Want[j] }

  init (n m) { State[n] = Idle && Own[n,m] = False } This way we do not put any constraints on link topology. The bad configuration are defined by graphs of the following form. unsafe (n m e) { State[n] = Busy && State[m] = Busy && Link[n,e] = True && Link[m,e] = True }

  X), ¬link(X, E) ∅ {link(X, E)} U nlink link(X, E) {link(X, E)} ∅In this model we assume that nodes have some form of control over connections with buffers (i.e. a new link is detected by a node in state Idle). When the model extended with the link and unlink rules is checked with Cubicle (see appendix B), the tool verifies correctness by visiting 28 nodes with at most 3 process indexes, invoking 261 times the SMT solver, performs 1449 fixpoint tests, and deletes 6 redundant nodes. The above result provides a formal correctness proof of the considered model for any number of nodes and links and any topology. Using dfs search, the tool verifies the property but the number of visited nodes is 39 with 414 calls to the SMT solver and 1996 fixpoint checks. Using the BRAB algorithm with parameter 3, the number of visited nodes reduces to 14 with 735 fixpoint tests, and 153 calls to the Alter-Ego solver. As for static topologies, the BRAB algorithm infers the invariant ¬(Own[#1, #3] = T rue ∧ Own[#2, #3] = T rue). A summary of the results with the considered models and heuristics are shown in Fig 1.

  denoted by the built-in type proc. Standard types int, real, and bool are also built in. Additionally, the user can specify abstract types and enumerations with simple declarations like type data and type msg = Empty | Req | Ack. As an example consider the following declaration.

	. Cache coherence protocols and mutual exclusion algorithms are typical
	examples of such systems. Cubicle model-checks by a symbolic backward reach-
	ability analysis on infinite sets of states represented by specific simple formulas,
	called cubes. Cubicle is written in OCaml. The SMT solver is a tightly integrated,
	lightweight and enhanced version of Alt-Ergo [20]; and its parallel implementa-
	tion relies on the Functory library [21]. Cubicle input language is a typed version
	of Murphi similar to the one of Uclid. A system is described in Cubicle by: (1) a
	set of type, variable, and array declarations; (2) a formula for the initial states;
	and (3) a set of transitions. It is parametrized by a set of process identifiers,

  Experimental results: V=visited nodes, F=fixpoint tests, S=solver calls, M=max process number, D=deleted node, I=number of invariants (brab), C=property checked (Yes/No).

	Model	dfs brab(3) V	F S M D I C
	DPP " "	√	√	19 529 176 3 0 0 Yes 19 529 176 3 0 0 Yes 10 312 88 3 0 1 Yes
	DPP+Link+Unlink " "	√	√	38 379 528 3 10 0 No 63 2660 961 3 12 0 No 33 365 617 3 9 1 No
	DPP+iLink+Unlink " "	√	√	28 1449 261 3 6 0 Yes 39 1996 414 3 6 0 Yes 14 735 153 3 2 1 Yes
	Fig. 1.


Lisa B Weiss

Vincent Dahirel
email: vincent.dahirel@sorbonne-universite.fr

Virginie Marry

Marie Jardat
email: marie.jardat@sorbonne-universite.fr

Computation of the Hydrodynamic Radius of Charged Nanoparticles from Non-equilibrium Molecular Dynamics

HAL is a

Introduction

Charged nanoparticles are used in numerous technologies. Modeling the properties of dispersions of nanoparticles plays an important role to design applications. [START_REF] Stark | Industrial applications of nanoparticles[END_REF] The numerical simulation of their properties usually requires the use of coarse-grained models, especially when nanoparticles are dispersed in aqueous solutions with added salt. Effective descriptions of the system averaged over water molecules and ions are justified, since the time scale of the nanoparticles transport is frequently much larger than typical solvent time scales.

However, the definition of the parameters of coarse-grained models is often a complicated task. Indeed, in order to perform simulations of the dynamical properties of nanoparticles in solution, one needs to determine their effective size and the behavior of water molecules in the vicinity of their surface. Both properties control the hydrodynamic radius of the nanoparticle. In most mesoscopic simulations techniques, like Brownian Dynamics (BD), [START_REF] Ermak | Brownian dynamics with hydrodynamic interactions[END_REF][START_REF] Dahirel | Two-scale Brownian dynamics of suspensions of charged nanoparticles including electrostatic and hydrodynamic interactions[END_REF] Multi-Particle Collision Dynamics (MPCD) [START_REF] Padding | Hydrodynamic interactions and Brownian forces in colloidal suspensions: Coarse-graining over time and length scales[END_REF][START_REF] Gompper | Multi-particle collision dynamics: A particle-based mesoscale simulation approach to the hydrodynamics of complex fluids[END_REF] or Dissipative Particle Dynamics (DPD), [START_REF] Español | Perspective: Dissipative particle dynamics[END_REF] the boundary condition for the velocity of the solvent at the nanoparticles surface can be modeled as a slip or a stick boundary condition. For colloidal particles, the definition of the hydrodynamic radius and of the boundary condition are relatively easy, because the structural radius, deduced for instance from microscopy or other structural investigations often matches the one deduced from the measurement of the diffusion coefficient using the Stokes-Einstein law (see for instance [START_REF] Kops-Werkhoven | Concentration dependence of the self diffusion coefficient of hard, spherical particles measured with photon correlation spectroscopy[END_REF] ). For nanoparticles, whose size is not large compared to the solvent size, the whole definition of the hydrodynamic radius is not obvious anymore.

Does the solvent perfectly stick on the surface, or partially slip [START_REF] Schmidt | Hydrodynamic boundary conditions, the Stokes-Einstein law, and long-time tails in the Brownian limit[END_REF][START_REF] Schmidt | Brownian motion of a rough sphere and the Stokes-Einstein law[END_REF][START_REF] Morrone | Interplay between hydrodynamics and the free energy surface in the assembly of nanoscale hydrophobes[END_REF][START_REF] Ishii | Molecular insights into the boundary conditions in the Stokes-Einstein relation[END_REF] ? If it sticks, is there an immobile layer of molecules at the surface of the nanoparticle that modifies the effective hydrodynamic size of the particle ? Does the charge influence the friction on the nanoparticle ? To answer these questions, modeling the system at the atomic scale is mandatory.

For simple ions, the interpretation of individual transport coefficients in term of an effective radius raises many questions. [START_REF] Pau | Application of Stokes law to ions in aqueous solutions[END_REF] It is known that the effective hydrodynamic radius depends on all interactions with the solvent molecules, and in particular that electrostatic forces give rise to a dielectric contribution to the friction. [START_REF] Hubbard | Dielectric dispersion and dielectric friction in electrolyte solutions[END_REF][START_REF] Wolynes | Dynamics of electrolyte solutions[END_REF] In the present study, our goal is to compute the flow of water molecules around a charged nanoparticle from non-equilibrium (NEQ) atomistic simulations, in order to derive quantities like the hydrodynamic radius that can be used afterwards in coarse-grained simulations, and that can be used to challenge typical electrolyte or colloidal concepts and theories for small nanoparticles. Such non-equilibrium molecular dynamics simulations have to be carried out and analysed with great care.

A frequently employed strategy in molecular dynamics to create a solvent flow is to apply an external force to the system, either on all solvent atoms, like in Poiseuille flow, or only on charged particles in the solvent as in the case of electroosmosis. [START_REF] Yoshida | Molecular dynamics simulation of electrokinetic flow of an aqueous electrolyte solution in nanochannels[END_REF][START_REF] Bhadauria | Multiscale modeling of electroosmotic flow: Effects of discrete ion, enhanced viscosity, and surface friction[END_REF] In both cases, some energy is added to the system, and is converted into viscous energy. It results in heating of the system. Thus, a thermostat must be used in order to extract the added energy and to ensure a constant temperature during the simulation. Simulations of water flows between parallel walls, or parallel layers of inorganic materials modeled at the atomic level are well described in the literature. [START_REF] Joly | Liquid friction on charged surfaces: From hydrodynamic slippage to electrokinetics[END_REF][START_REF] Botan | Hydrodynamics in clay nanopores[END_REF][START_REF] Botan | How electrostatics influences hydrodynamic boundary conditions: Poiseuille and electro-osmostic flows in clay nanopores[END_REF][START_REF] Bonthuis | Beyond the continuum: How molecular solvent structure affects electrostatics and hydrodynamics at solid/electrolyte interfaces[END_REF] Several studies have focused on the different thermostat strategies for the simulation of fluids confined in a channel. [START_REF] Bernardi | Lyapunov spectra and conjugate-pairing rule for confined atomic fluid[END_REF][START_REF] Yong | Thermostats and thermostat strategies for molecular dynamics simulations of nanofluidics[END_REF][START_REF] De Luca | A new and effective method for thermostatting confined fluids[END_REF] In such cases, one can take advantage either of the high symmetry of the solvent flow or of the presence of the numerous atoms of the walls to thermostat the system. For example, an efficient strategy consists in coupling only the wall atoms with the thermostat. [START_REF] De Luca | A new and effective method for thermostatting confined fluids[END_REF] However, when only one fixed obstacle is put in the center of a simulation box filled with solvent molecules, the solvent velocity field varies in all directions and the use of usual thermostat strategies is tricky. We focus here on two different procedures to thermostat the simulation box: The Profile Unbiased Thermostat [START_REF] Evans | Nonlinear-response theory for steady planar Couette flow[END_REF] (PUT) and the Partial Thermostat (PT). The principle of the PUT is the following. The simulation box is divided in cells, and velocities of solvent molecules relative to the average solvent velocity in this cell are rescaled so that the target temperature is obtained, via a Nosé-Hoover algorithm. [START_REF] Martyna | Nosé-Hoover chains: The canonical ensemble via continuous dynamics[END_REF] When the symmetry of the velocity field is low, using PUT demands a particular study of the choice of the size of the cells. If they are too large, the flow field varies too much inside the cells. If they are too small, the number of degrees of freedom to be thermostatted together is too small. The other procedure, called in what follows the Partial Thermostat (PT) is simpler in principle: Only a fraction of the simulation box is thermostatted, and only part of the degrees of freedom of water molecules are thermostatted. This is similar to the aforementioned thermostatting of the atoms of a wall. [START_REF] De Luca | A new and effective method for thermostatting confined fluids[END_REF] The challenge in this case is to define a region where the geometry of the velocity field is rather uniform, and to check that this region is large enough to thermostat the rest of the simulation box.

In this paper, we present a robust way to study the flow around a fixed, approximately spherical nanoparticle from non-equilibrium molecular dynamics. We thoroughly studied the efficiency of the PUT and of the PT in this case, and defined the conditions where they lead to a proper control of the temperature of the system. From these NEQ simulations, we computed the friction coefficient exerted by the solvent on the nanoparticle. Then, we deduced the hydrodynamic radius of the nanoparticle by using the approximate analytical result of Hasimoto, [START_REF] Hasimoto | On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres[END_REF][START_REF] Yeh | System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions[END_REF] which gives the flow field across an infinite periodic array of spheres. Our simulation strategy is applicable to any nanoparticle that can be approximated by a sphere. Nevertheless, as water-interface interactions are highly dependent on subtle surface details, we focused here on a specific nanoparticle. Indeed, even charged generic nanoparticles can have spurious hydrophobic properties. [START_REF] Dzubiella | Competition of hydrophobic and coulombic interactions between nanosized solutes[END_REF] We focused on polyoxometalate ions (POM), which are inorganic charged nanoparticles used as a standard in electrokinetic measurements, especially in electroacoustic experiments. [START_REF] Hodne | Verification of the electroacoustic calibration standard: comparison of the dynamic mobility of silicododecamolybdate and silicododecatungstate Acids and Salts[END_REF] More precisely, we took the example of the phosphotungstate anion [PW 12 O 40 ] 3-. This family of systems is already well studied in the literature, but a description of their dynamic characteristics at the atomic scale is still lacking. Once the simulation procedure has been validated, we could make a preliminary study of the influence of the charge of the nanoparticle on its dynamic properties. We computed the friction coefficient and the hydrodynamic radius of a neutralized POM, and of nanoparticles of charge -3e -4e, and -5e while keeping the surface charges unchanged.

We also studied the influence of the presence of an added salt on the hydrodynamic radius of these nanoparticles.

The paper is organized as follows. Section II describes the force fields and the procedures used to impose a flow and to compute the friction coefficient. In Section III we compare the efficiency of the profile unbiased thermostat (PUT) and that of the partial thermostat procedure. The calculation of the friction coefficient of the inorganic ion is discussed in Section IV. [START_REF] Maestre | Electronic and magnetic properties of R-Keggin anions : A DFT[END_REF] Some of these authors developed force fields for Keggin anions containing phosphorus, silicon and aluminium as central atom. [START_REF] López | Polyoxometalates in media[END_REF][START_REF] Leroy | Keggin polyoxoanions in aqueous solution: ion pairing and its effect on dynamic properties by molecular dynamics simulations[END_REF] The proposed force fields were successfully engaged to analyse ion pair formation between these POMs and potassium, sodium, and lithium as counterions. For a box size with 1000 water molecules and one POM compared to 8000 water molecules and 8 POMs no size dependence of the radial distribution functions was noticed. [START_REF] Leroy | Keggin polyoxoanions in aqueous solution: ion pairing and its effect on dynamic properties by molecular dynamics simulations[END_REF] These force fields were also successfully employed to study the aggregation of multiples POMs as function of counterions, charge and solvent. [START_REF] Chaumont | Ion aggregation in concentrated aqueous and methanol solutions of polyoxometallates Keggin anions: the effect of counterions investigated by molecular dynamics simulations[END_REF][START_REF] Chaumont | Interactions between Keggin anions in water: The higher their charge, the higher their condensation? A simulation study[END_REF] Therefore, we have taken the force field parameters from Ref. [START_REF] Leroy | Keggin polyoxoanions in aqueous solution: ion pairing and its effect on dynamic properties by molecular dynamics simulations[END_REF] The parameters of the force field are given in Table 1. As the coordinates of the DFT optimized POM structure used in Ref. [START_REF] Leroy | Keggin polyoxoanions in aqueous solution: ion pairing and its effect on dynamic properties by molecular dynamics simulations[END_REF] are not available, the crystal coordinates of H 3 PW 12 O 40 • 6H 2 O were taken from reference. [START_REF] Noe-Spirlet | Neutron diffraction studies of three hydrates of phosphotungstic acid[END_REF] We checked that the intramolecular radial distribution functions are in very good agreement with those reported in Ref. [START_REF] Leroy | Keggin polyoxoanions in aqueous solution: ion pairing and its effect on dynamic properties by molecular dynamics simulations[END_REF] The extended simple point charge model of water (SPC/E) was employed. [START_REF] Berendsen | The missing term in effective pair potentials[END_REF] The force field parameters for potassium and chloride ions were taken from references. [START_REF] Lee | Molecular dynamics simulation of ion mobility. 2. Alkali metal and halide ions using the SPC/E Model for water at 25C[END_REF][START_REF] Koneshan | Solvent structure, dynamics, and ion mobility in aqueous solutions at 25C[END_REF] Lorenz-Berthelot rules were used to calculate the interatomic interactions for different atom types. All simulations were performed using the molecular dynamics simulation package LAMMPS. 39,40 The SHAKE algorithm, which is implemented in LAMMPS as described in reference 41 was used to restrict bond lengths and angles of the water molecules. Moreover, the POM was kept fixed in the center of the box during the equilibrium and non-equilibrium simulations.

Methods

Modeling of the system

The forces acting on all the atoms of the POM were thus set to zero.

Long range coulomb interactions were taken into account via a particle-particle particle mesh solver. A cut-off radius of 15.0 Å was set for the Lennard-Jones interactions and the real space Coulomb interactions. For each initial configuration of the system, a combination of equilibration runs in the N V T and N pT ensembles were performed to obtain bulk water at the expected density, for approximately 1 ns, with a time step δt = 1.0 fs. The target values were T = 300 K and p = 1.0125 MPa. During equilibration, pressure and temperature control were reached by a Nosé Hoover barostat and thermostat. These are implemented in LAMMPS as described in reference. [START_REF] Martyna | Constant pressure molecular dynamics algorithms[END_REF][START_REF] Martyna | Nosé-Hoover chains: the canonical ensemble via continuous dynamics[END_REF] No thermostat was used for ions. We verified that ions were effectively thermostatted by the surrounding solvent.

To induce a flow of water molecules in the simulation box, two strategies are described in the literature: Either adding a force to the fluid particles (within a slice or for all fluid molecules), or imposing a constant velocity to part of the atoms of the fluid. We have applied a constant force to all atoms in the x-direction, which is equivalent to a pressure driven flow.

We have varied the force until the flow corresponds to a low Reynolds number R e 1, defined as:

R e = v H 2 O×d P OM ν (1)
with a diameter of the POM d P OM ≈ 12 Å and the kinematic viscosity denoted by ν. We have then ensured that the flow is in the linear response regime. The characteristic velocity of the fluid v H 2 O can be defined as v

H 2 O = | Fx | ξ S
, where | F x | is the total external force acting on solvent molecules along the x-axis, and is therefore an input of the simulation, and ξ S is the friction felt by a perfectly spherical obstacle of radius R with stick boundary conditions computed by the Stokes law ξ S = 6πηR. If we take the viscosity of SPC/E water (η water = 0.68 × 10 -3 Pa.s at 300 K, see reference [START_REF] Tazi | Diffusion coefficient and shear viscosity of rigid water models[END_REF] ), and a force equal to f = 0.1 × 10 -13 N per water molecule (N 1700 water molecule in the simulation box), we obtain a Reynolds number R e 0.003.

During the non-equilibrium simulations, the central nanoparticle is fixed in the center of the simulation box and is not allowed to rotate. We have checked that the velocity profiles of water molecules are the same whatever the orientation of the nanoparticle within the statistical error.

Use of a thermostat

Our first goal is to establish a strategy to thermostat a flow around an approximately spherical obstacle in molecular dynamics simulations. Hence, we simulate a system out of equilibrium, but still in the linear response regime. As already mentioned, such non-equilibrium molecular dynamics requires special thermostatting strategies. In this study, NEQ simulations were performed employing two different thermostatting strategies, namely the profile unbiased thermostat (PUT) [START_REF] Evans | Nonlinear-response theory for steady planar Couette flow[END_REF] and a partial thermostat (PT).

Profile Unbiased Thermostat

The PUT has been developed by Evans et al. [START_REF] Evans | Nonlinear-response theory for steady planar Couette flow[END_REF] to thermostat a turbulent flow of Lennard-Jones particles. In this procedure, the simulation box is divided into a fixed number of bins.

Subsequently, the center of mass velocity of every bin is calculated. The velocities of the solvent particles minus the center of mass velocity are computed, and these velocities are employed for temperature definition. [START_REF] Evans | Nonlinear-response theory for steady planar Couette flow[END_REF] More precisely, for every bin of position r and at each time step t, the local kinetic energy E kin (r, t) is calculated by summing the square of all velocities v i minus the average streaming velocity u over all atoms N in the bin, which allows one to compute a local and instantaneous temperature T (r, t):

E kin (r, t) = (D • n(r, t) -d) k B T (r, t) 2 (2) = N i=1 1 2 m i [v i (r, t) -u(r, t)] 2 δ(r i (t) -r) (3)
Here the number of degrees of freedom per particle D is multiplied by the instantaneous number of molecules n(r, t) in the bin. The dimensionality d of the system is subtracted to account for the conservation of the d components of the momentum of the simulation box; k B is Boltzmann constant. The number of degrees of freedom per particle is D = 6

for the rigid SPC/E water model in combination with the SHAKE algorithm. The delta function δ(r) ensures that all atoms are within the specific bin. The system temperature is then calculated, and the Nosé-Hoover thermostat with a relaxation time τ =1000 fs is used.

We have checked that our results were not influenced by the value of the relaxation time:

a value twice as large leads to the same results within the statistical error. In the PUT procedure, the number of bins in the simulation box is a crucial parameter: The resolution of the grid should be high enough to have a small velocity gradient within every bin, but small enough to ensure that thermal fluctuations are statistically significant inside each bin.

This common idea is a classical result for random variables, as quantified for instance by Chebyshev's inequality for the random variable [v i (r, t) -u(r, t)] 2 . To illustrate this, one can simplify the above equation writing for a fluid at rest and particles of equal masses m i :

T (r, t) = A N N i=1 v i (r, t) 2 δ(r i (t) -r)
where A does not depend on N . Chebyshev's inequality then reads:

P (|T (r, t) -T | > δv 2 ) < 1 N ( 4
)
where P is a probability, T is the temperature of the system (i.e the temperature computed for a large number of degrees of freedom), and δv 2 the standard deviation of the individual square velocities v i (r, t). T (r, t) converges towards a good measure of the exact temperature

T for large values of N (large bins), for a fluid at rest.

Partial Thermostat

Another way to regulate the temperature of the simulation box during NEQ simulations consists in using a partial thermostat. Typically, the simulation box is divided in two parts:

The first one is thermostatted by using a NH algorithm and the second one is not thermostatted. It means that molecules located in the second part of the simulation box would evolve in a (N V E) ensemble if no heat transfer were allowed with the first part of the box.

This strategy has been used for simulating a flow in a porous media, for which the atoms of the solid walls are thermostatted, but the atoms of the fluid are not. [START_REF] De Luca | A new and effective method for thermostatting confined fluids[END_REF] In our case, we chose to thermostat the fluid in a region where the averaged fluid velocity is very weak in the x and y directions. Therefore, the thermostatted region should be as far as possible from the spherical obstacle, and sufficiently large to contain a significant number of fluid molecules.

More precisely, the simulation box has been separated in two parts: (1) The molecules that are within a cylinder of radius R cyl , the axis of which is parallel to the x-axis passing through the central atom of the POM and the length of which is equal to that of simulation box, and

(2) the water molecules that are outside the cylinder. The molecules that are outside the cylinder are thermostatted by a Nosé Hoover (NH) thermostat restricted to the two components of the velocity perpendicular to the flow (v y and v z ). For such a thermostat (T yz ) the number of degrees of freedom is 4 for rigid water molecules, so that the total number of degrees of freedom to be accounted for by the thermostat is 4 n molec,out -2. This expression takes into account the fact that the overall transport in y and z direction should be zero due to the geometry of the flow profile in the thermostatted region as soon as R cyl is high enough.

Computation of the friction coefficient of the solute

The friction coefficient ξ of a solute can be defined as the ratio of the force exerted by the fluid on the solute to the characteristic velocity of the fluid. This force is indeed proportional to the fluid velocity in the linear regime, under the conditions of a low Reynolds number.

The friction on a sphere is unambiguously defined for a sphere immersed in a fluid at infinite dilution, when the fluid is moving along the x-direction far from the object with a velocity v x,∞ :

ξ = | F x | v x,∞ (5)
In eq. 5, F x is the total force acting on the sphere, averaged over the molecular dynamics trajectory. In a low Reynolds number flow, for stick boundary conditions, the classical Stokes calculation yields ξ = 6πηR, where R is by definition the hydrodynamic radius of the spherical obstacle. [START_REF] Guyon | Hydrodynamique physique[END_REF] In our case, there is one solute fixed in the center of the simulation box with periodic boundary conditions, as it is commonly done in molecular dynamics, so that the solute concentration is finite. None of the fluid molecules inside the simulation box can be considered to be far from the solute. However, the long-range nature of the hydrodynamic interactions influences dynamical quantities such as the Stokes friction, even when structural equilibrium properties are close to those of an infinite fluid around a unique spherical object. [START_REF] Dünweg | Molecular dynamics simulation of a polymer chain in solution[END_REF] We are actually dealing with a flow of solvent in a periodic cubic array of spherical objects. Since the fluid is incompressible, the total mass flow across the (y, z) planes perpendicular to the flow is conserved for all values of x. In the (y, z) planes that do not cross the spherical obstacle, the average fluid velocity v av does not depend on the value of x. The friction coefficient ξ can then be unambiguously defined as the ratio between the averaged fluid velocity v av and the total applied force | F x | :

ξ = | F x | v av (6)
In this case, the friction coefficient ξ in the Stokes regime does not scale as the hydrodynamic radius. Corrections related to the size of the particle and to the distance between the periodic images (here the box length L) must be taken into account. This hydrodynamic problem has been studied in various works, and the impact of periodic images on friction has been first estimated by Hasimoto. [START_REF] Hasimoto | On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres[END_REF] In the case of a periodic cubic array of spheres of same radius R, with R being very small compared to the distances between spheres R L, one can neglect quadratic corrections and gets the following expression of the force acting on a sphere:

F = 6πηR 1 + 2.8373R/L v av (7)
In the regime where R is not very small compared to L, Hasimoto derives a more precise analytical expression of the force (see Section 5 of Hasimoto's paper [START_REF] Hasimoto | On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres[END_REF] ). In what follows, we have used this refined expression to deduce the hydrodynamic radius of the charged nanoparticle from our non-equilibrium simulations. Indeed, as we perform explicit solvent calculations with many solvent molecules, we are constrained in a regime where the size of the simulation box is not much higher than the size R of the nanoparticle. Then, the linear correction in R/L box is not accurate, and one needs to consider higher order terms.

To compute the friction coefficient from the simulations, we have to evaluate the total force acting on the nanoparticle. When the flow is stationary, the force acting on the nanoparticle is exactly equal to the total external force applied to solvent molecules. Indeed, when the fluid does not accelerate, the total sum of all forces on the fluid should be zero. The forces acting between water molecules, between ions, and between ions and water molecules cancel out due to Newton's third principle. The remaining contributions are thus the forces acting on the solvent due to the presence of the central nanoparticle, and the external forces.

These two forces are then opposite. The force felt by the nanoparticle should therefore be exactly equal to the external force. This result was numerically confirmed in our simulations. As we proceed to show in Sec. , the chosen thermostat strategy has no influence on the balance of forces.. In practice, we have first computed the friction coefficient as the ratio between the total force acting on the nanoparticle and the average fluid velocity v av .

We have then deduced the hydrodynamic radius of the nanoparticle by using Hasimoto's expression given in Section 5 of reference. [START_REF] Hasimoto | On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres[END_REF] In the following, we define a characteristic friction ξ o as the friction felt by a perfectly spherical obstacle of radius 0.6 nm in a periodic cubic array of period 3.7 nm, which is the value of the simulation box length in our calculations. With a viscosity equal to 0.68 × 10 -3

Pa.s, which is the viscosity at 300 K of pure SPC/E water, [START_REF] Tazi | Diffusion coefficient and shear viscosity of rigid water models[END_REF] we obtain the value ξ o = 1.381 10 -11 kg.s -1 using Hasimoto's expression given in Section 5 of reference. [START_REF] Hasimoto | On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres[END_REF]

Results

Efficiency of the thermostat to compute the flow around a nanoparticle Profile Unbiased Thermostat

In this part, we test the ability of the Profile Unbiased Thermostat to regulate temperature so that the flow of solvent in non-equilibrium molecular dynamics is correct. The flow depends on the temperature for various reasons. The most obvious one is that the fluid viscosity highly depends on temperature. The effect of the thermostat on the fluid viscosity is thus indirectly evaluated through the measure of the average velocity v av in the (y, z) planes perpendicular to the flow. Indeed, for a given external force, this average velocity is inversely proportional to the friction coefficient ξ, which is itself exactly proportional to the viscosity. Therefore, in a Stokes flow, the viscosity is inversely proportional to the average velocity v av .

As already mentioned, the correct use of the PUT requires to choose the number of bins used to compute the local temperature carefully. In this section, we report the results we obtained for a flow of water around a POM nanoparticle at small Reynolds numbers, using the PUT procedure with several resolutions of the grid of bins. In this case, as the POM is described at the atomic scale, the velocity field of water is highly inhomogeneous in the vicinity of the POM. In principle, we should use small bins in the PUT to account for spatial variations of the velocity field. But, there should be enough molecules per bin so that thermal fluctuations inside each bin converge close to the average fluctuations in a bulk fluid at the correct temperature. We therefore need to evaluate the effect of the bin size on the flow for reasonable values of the bin size.

As far as the flow inhomogeneity is concerned, choosing bins with a size of half the box length L box , i. e. taking a grid of 2 3 bins, does not improve the problem, because of obvious symmetry reasons. We thus chose bins of size equal to L box /4 as the largest possible bins. In our case, there are about 1700 water molecules in a cubic box of about 37 Å length, so that a grid of 10 3 bins corresponds to about 2 water molecules per bin in average. For 8 [START_REF] Dahirel | Two-scale Brownian dynamics of suspensions of charged nanoparticles including electrostatic and hydrodynamic interactions[END_REF] bins there are about 3 molecules per bins, 8 for 6 3 bins, and 27 molecules for 4 3 bins. We have chosen as greatest resolution, i. e. as the smallest bins, a case where the average number of water molecules is greater than 2, else the calculation makes no sense at all. Finally, we have thus restricted our study to the cases with 8 3 , 6 [START_REF] Dahirel | Two-scale Brownian dynamics of suspensions of charged nanoparticles including electrostatic and hydrodynamic interactions[END_REF] , and 4 3 bins. In all described simulations a force of f = 0.1 × 10 -13 N was added to all water molecules. We give in Fig. 2 the velocity profiles of water molecules as a function of the x coordinate, for different resolutions of the PUT grid. The results weakly depend on the grid resolution.

The standard uncertainty of the average velocity has been computed from a block analysis of the fluctuations of the velocity of the fluid over several temporal averages of 4 ns. It is always smaller than 0.35 m.s -1 so that the error estimate on v/v o is equal to 0.03. When the grid size changes from 4 3 to 6 3 bins, the differences in the averaged flow velocity v av are smaller than these fluctuations. There is no exact analytical result with which we could directly compare the values of the velocity v av . Nevertheless, if the POM nanoparticle were a perfectly spherical particle of radius equal to 6.0 Å with stick boundary conditions, one would get v av = v o , where v o is computed for the Stokes flow within a cubic array of spherical obstacles, with spheres of radius 6.0 Å close to the radius of POM. The results are close to this perfectly spherical case. The velocities obtained with thermostat grids of 8 3 bins are nevertheless smaller than the velocities obtained with larger bins. This difference is slightly larger than the uncertainty of the velocity. It can be explained by the small size of the bins.

As indicated before, with such a small number of water molecules in each bin (about 3), the thermal fluctuations might not be significant from the statistical point of view.

There are then two indications showing that the PUT is adequate to regulate the temperature of the fluid. Nevertheless, using a Nosé Hoover algorithm close to the nanoparticle might change the local dynamics, such as for instance the local diffusional dynamics of water and ions. If one is interested in the surface conductivity or other surface phenomena, it is useful to see if another thermostatting strategy can be used which does not affect the area close to the solid-liquid interface.

Partial thermostat with a cylindrical geometry

We use in this part a partial thermostat with a cylindrical geometry. In parts of the simulation box far from the obstacle, outside of the cylinder, the flow field is mainly oriented towards the direction of the applied force. This subpart is thermostatted in the directions Figure 2: Velocity profiles of water molecules along the x-direction for different resolutions of the PUT grid. The reference velocity v o is the average velocity around a perfectly spherical obstacle periodically fixed within a cubic array. [START_REF] Hasimoto | On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres[END_REF] The center of the POM is located at x = 0. perpendicular to the flow. This requires to adequately choose the radius of the cylinder.

As it can be seen on Fig. 2, for x coordinates between -18 Å and -10 Å the velocity of the fluid in the x direction is almost constant, with variations below the uncertainty of our calculations (the center of the obstacle is in the position x = 0).

The tendency is the same for all regions of space further than 15.0 Å from the center of the obstacle, and this will be checked a posteriori. We choose thus to thermostat the solvent outside a cylinder of radius equal to 17.0 Å.

First, this cylindrical thermostat was used for equilibrium simulations, with a target temperature of 300 K. A perfect agreement was obtained between the time averaged temperature within this procedure and that obtained with a NH thermostat on the full simulation box.

We turned then to non-equilibrium simulations. A 2 ns non-equilibrium run was done to reach steady state before the production run started. At steady state, we expect that the Finally, we have checked that the density distribution functions around the spherical obstacle correspond to that of an equilibrium simulation at the right temperature (300 K), as the local temperature is expected to affect water structure. Equilibrium simulations for T = 350 K were also done. The significantly higher temperature of 350 K has been chosen to obtain a quantitative effect on distribution functions, one order of magnitude higher than noise. We compare the results obtained from these two sets of equilibrium simulations with non-equilibrium simulations using the Partial Thermostat and the Profile Unbiaised Thermostat. We show on Fig. 4 the distribution functions between hydrogen atoms of water and the phosphorus atom of the POM. They have been computed in a slice perpendicular to the x direction with a width of d bin = 3.1 Å. Within the slice, the rdf is computed in cylindrical bins of length d bin = 0.373 Å centered on the POM, as a function of r yz = y 2 + z 2 . As expected, increasing temperature makes water molecules more mobile and less structured: The maximum of the second peak of the distribution function decreases strongly with increasing temperature. The distribution function obtained from non-equilibrium simulations coincides very well with that obtained at equilibrium at the same temperature. In particular, the difference at the second peak (r yz around 7.6 Å) is much weaker than the difference with the distribution function obtained at 350 K. The agreement is observed whatever the size of the grid when the PUT is applied. This also confirms that there is no significant heating of the system at the solid/water interface with the partial thermostat, even if the thermostat only works far from this interface.

In In what follows, we have then only used the partial thermostat with a cylindrical symmetry. This thermostat only affects the components of water velocities in the direction perpendicular to the flow (T yz ), and far from the central nanoparticle. This also ensures that the balance of forces on the x-direction, which we use to compute the friction coefficient (see eq. 6), is not influenced by the thermostat. From equilibrium simulation at 300 K (black line), from equilibrium simulation at 350 K (red dotted line), and from non-equilibrium simulations with a partial thermostat with a cylindrical frontier (blue circles) at 300 K and with a PUT with 4 3 or 8 3 bins.

Friction coefficient of charged inorganic nanoparticles

Using the cylindrical partial thermostat, we have computed the friction coefficient of the POM and of inorganic nanoparticles analogous to the POM but with other charges, also in the presence of added salt, from non-equilibrium simulations. From this quantity, as described before, we have used Hasimoto's formula to deduce the hydrodynamic radius of the nanoparticle, assuming stick boundary conditions.

Ten different systems were studied, which are listed in Table 2. The exact number of water molecules in the simulation box and the precise size of the box, obtained after equilibration runs in the (N pT ) and (N V T ) ensembles are also given. The counterion of the POM is potassium K + in every case. In some cases, potassium chloride is added, with a total number of potassium ions always equal to 15, so that the electrostatic screening is the same for all systems. This corresponds to a concentration of c K + = 0.48 mol L -1 ≈ 0.5 mol L -1 .

The number of chloride ions is adapted to ensure electroneutrality. To obtain nanoparticles of different charges, the charge of the central atom is increased or reduced compared to that given in Table 1. The goal is to change as less as possible the surface charges of the POM so that differences emerge only as a function of the total structural charge. The structure of the POM is kept as described in Table 1. For the systems POM5, POM5-salt, POM0 and POM0-salt, an alternative method to change the total charge is also tested: The charge of the phosphorus and the charges of the tungsten were increased (or reduced) by the same amount. The radial distribution functions between the POM and water atoms for these two different charge distributions show no significant difference. It should be stressed that POM3 and POM4 ions exist, respectively with the structures [PW 12 O 40 ] 3-and [SiW 12 O 40 ] 4-. Similar systems with a radius of about 0.5 nm exist also for charges -5 e, -6 e and higher charges. [START_REF] Long | Towards polyoxometalate-integrated nanosystems[END_REF] The SPC/E water model was chosen with a time step of 1 fs in combination with the SHAKE algorithm. Several equilibration simulations were run, 100 ps in N V T , 1000 ps in N pT and another 100 ps in N V T . Subsequently the POM was fixed in the center of the simulation box setting all forces and its rotational velocity to zero at every time step.

Another equilibrium run was done to generate different starting configurations for the nonequilibrium simulations. For all systems, a flow was applied during a 2 ns simulation, prior to the production run, in order to establish steady state conditions. These systems were then simulated at least over 37 ns in total, if we add the simulation time starting from different configurations. To derive the friction coefficient, we have computed the averaged solvent velocity far from the POM, as described before.

The values of the friction coefficient obtained for the different nanoparticles are shown in Fig. 5 and compared to the values obtained in the presence of an added salt. The friction coefficient is divided by ξ 0 , which is the friction coefficient that would be obtained for an uncharged spherical obstacle of radius equal to 6.0 Å at infinite dilution in water with stick boundary conditions. The values of the hydrodynamic radius deduced from the friction coefficient using Hasimoto's result using the method detailed in the methodological section are given in Fig. 5.

The hydrodynamic radius of the neutral nanoparticle is found equal to 6.45±0.01 Å. The standard uncertainty is evaluated from the standard uncertainty of the velocity deduced from a block analysis. For the POM with charge -4 e, the hydrodynamic radius is found equal to 6.75 ± 0.01 Å. This value is close to our recent experimental study of silicotungstate ions. The electroacoustic signal of aqueous solutions of silicotungstate was indeed measured and analyzed with a novel theoretical treatment in terms of charge and hydrodynamic radius: [START_REF] Pusset | Nonideal effects in electroacoustics of solutions of charged particles: combined experimental and theoretical analysis from simple electrolytes to small nanoparticles[END_REF] The best fit of the experimental data corresponds to 6.30 Å.

As it is shown in Fig. 5, the friction coefficient of the nanoparticle increases with increasing charge, and as a consequence the hydrodynamic radius increases with increasing charge.

The increase of the computed hydrodynamic radius is significant, from 6.45 Å for the neutral nanoparticle to 7.24 Å for a charge equal to -5e without salt. It is rather surprising if one considers that the spatial coordinates of the atoms of the nanoparticle do not change at all. For such small nanoparticles, it is difficult to establish exactly what should be the value of the hydrodynamic radius from the microscopic structure. Nevertheless, it should be related to the average distance between the water molecules in the first solvation shell containing mobile water molecules and the center of the nanoparticle. Two effects can explain the variations of the hydrodynamic radius. First, water molecules are more attracted by the nanoparticle when its charge increases. This can be quantitatively observed in Fig. 6 which gives the number of water molecules as a function of the distance to the phosphorus atom of the POM. The number of water molecules as a function of the distance has been computed as an integral of the radial distribution function. The plateau in the coordination number plot as function of the distance for the systems POM4, POM5 indicates that the POM is solvated respectively by 8, and 10 water molecules when its charge increases from 4 to 5. On the contrary, no clear plateau appears for the POM3 which means that the solvation shell is less clearly defined than with larger charges. Therefore, the mobility of water molecules decreases with the charge of the nanoparticle, as it is shown in Fig. 7 that gives the mean velocity of water molecules as a function of the distance to the center of the POM nanoparticle. It appears clearly in Fig. 7 that the velocity of solvent molecules is almost zero (0.5 m.s -1 ) for systems POM4 and POM5 at a distance equal to 5 Å, with or without added salt, whereas it is about three times larger for POM3.

In order to compare the time spent by a water molecule around a nanoparticle depending Velocity v / v o q POM = -3 / no salt q POM = -4 / no salt q POM = -5 / no salt q POM = -3 / salt q POM = -4 / salt q POM = -5 /salt on its charge, residence times around POM3 and POM5 particles have been evaluated. The survival probability for a molecule to be found within a given distance d of the center of the nanoparticle can be calculated:

P (t) = S I (t) S I (0)
where S I (t) = 1 if the molecule is present in the sphere of radius d centered around the POM both at time 0 and t, and S I (t) = 0 otherwise. Such a definition allows for molecules to leave and come back into the sphere between 0 and t and characterizes an intermittent survival probability. d corresponds to the position of the first minimum of the radial distribution functions between the phosphorus atom and the oxygen atom of water molecules, i.e. 6.2 Å. P (t) can be fitted by a decreasing exponential function exp(-t/τ ), where τ is the residence time in the volume occupied by the sphere. More realistic ways to calculate residence times exist, which can be better compare with experiments. [START_REF] Laage | On the residence time for water in a solute hydration shell: Application to aqueous halide solutions[END_REF] However, in our case, the aim is to see how much the charge of the nanoparticle can slow down water molecules surrounding it. Therefore, the comparison between intermittent residence times calculated in the same manner already gives interesting insight. The probability distribution of residence times in a spherical volume of radius 6.2 Å are shown in Fig. 8. The average residence time of water around the nanoparticle is 28 ns in the case of the nanoparticle of charge -3e and 320 ns in the case of the nanoparticle of charge -5e. This confirms the structural behavior deduced from the coordination numbers: Water molecules stay a significantly much longer time close to the nanoparticle when it is more charged. As a consequence, it looks like the hydrodynamic radius is larger, accounting for water molecules that would be stuck at the surface of the nanoparticle. In recent simulations of electroosmosis in a slit geometry, such dependence of the flow at the surface of the solid on its charge has also been observed. [START_REF] Yoshida | Molecular dynamics simulation of electrokinetic flow of an aqueous electrolyte solution in nanochannels[END_REF] A second phenomenon which increases the hydrodynamic radius when the charge increases is the so-called electrostatic friction, arising from the interactions with both water molecules (dielectric friction) and small ions. The friction is influenced by the long-range electrostatic interactions, which can effectively increase the hydrodynamic radius. The Hubbard- Onsager theory [START_REF] Hubbard | Dielectric dispersion and dielectric friction in electrolyte solutions[END_REF][START_REF] Wolynes | Dynamics of electrolyte solutions[END_REF] yield an estimation of the increase of the friction due to the interaction with water dipoles. In a subsequent study, a detailed comparison of our results with this theory and a more exhaustive analysis of water relaxation around the POM will enable us to unravel whether the increase of the friction with charge can be quantitatively attributed to dielectric friction, or to the strong water adsorption. Please note that Hasimoto's treatment of the flow in an array of spherical obstacles does not take into account the influence of electrostatics on the local viscosity of water. We shall address this limitation in our forthcoming study.

The influence of salt on the friction coefficient is small but systematic: The friction in the presence of added salt is always a few percents larger than the friction without salt. As it can be observed in Fig. 7, the presence of added salt leads to a decrease of the velocity of water molecules compared to the case without salt. The decrease is of about 10%. This reflects a small enhancement of the effective viscosity of the solvent due to the presence of the added salt. This effect then slightly impacts the friction coefficient computed from the simulation (the friction increases slightly), which increases the effective hydrodynamic radius. This indirect effect of salt on the hydrodynamic radius is weak, but it is systematic.

In a previous work, Chowdhuri and Chandra 50 performed atomistic molecular dynamics simulations of aqueous solutions of potassium chloride with the same interaction potential as in the present work. They showed that the dynamics of water molecules slows down as the salt concentration increases. At a KCl concentration of 0.88 mol/L, they observe a decrease of the self-diffusion coefficient of water molecules of about 6%. We can thus argue that the presence of potassium and chloride ions (at a concentration of 0.5 mol/L in our simulations) impacts the dynamics of water molecules and yields a small increase of the effective viscosity of the solvent.

Conclusion

Our goal was to find a robust procedure to evaluate the hydrodynamic radius of a nanoparticle in water from numerical simulation. Indeed, this quantity is not always easy to deduce from experiments but must be known to perform mesoscopic numerical simulations of the dynamic properties of nanoparticle dispersions. Moreover, such study is very important to determine whether a specific nanoparticle is sufficiently large to behave as a typical colloidal particle from the dynamical point of view. For instance, the existence of a linear dependence between the size of the particle and its friction, or its diffusion coefficient, can be evaluated, i.

e. the relevance of the Stokes-Einstein relationship can be evaluated. [START_REF] Zhang | Hydrodynamic drag on diffusing nanoparticles for size determination[END_REF] Small nanoparticles are indeed in between colloids and electrolytes, and for electrolytes, the link between the structure (size, shape) and the friction is much less clear than for colloidal particles.

We have proposed to use non-equilibrium molecular dynamics simulations to study the flow of water around the nanoparticle, and then to compute the friction coefficient of the nanoparticle. One difficulty in non-equilibrium simulations is to find a good way to induce the flow without heating the system. As we have shown in the present study, this difficulty is more pronounced for the flow around a spherical nanoparticle because the geometry of the velocity field is less symmetric than in a channel. We compared two procedures to thermostat the simulation box, both available as options in common open-source molecular dynamics simulation packages. The first one, the PUT, divides the box volume into small bins, in which the average velocity of the solvent molecules is approximately uniform. The second one consists in thermostatting the system only in a subpart of the simulation box, where the velocity field is almost uniform. We have shown that for the flow around a roughly spherical object the choice of a cylindrical shape of the boundary between the thermostatted volume and the rest of the box allows one to obtain a constant temperature in all the simulation box.

We have then performed non-equilibrium simulations of the flow of water around nanoparticles of different charges using a partial thermostat of cylindrical shape. The atomic struc-ture of the nanoparticle represents the one of typical polyoxometalate ions (POM). The friction coefficient of the nanoparticle was deduced from the average solvent velocity with a good precision. This result enabled us to compute the hydrodynamic radius by using Hasimoto's approximate analytical result for the flow across a cubic array of spheres. We found that the hydrodynamic radius significantly increases with the charge of the nanoparticle, a phenomenon that had not been quantified so far using molecular dynamics. This result suggests that such nanoparticle should be seen as a big ion more than as a small colloidal particle. It should be of particular interest for the analysis of routine experiments on nanoparticles, such as dynamic light scattering or laser zetametry when the diffusion coefficient or the electrophoretic mobility are generally evaluated using the Stokes-Einstein estimation of friction, which only depends on the radius of the nanoparticle and not directly on its charge. While more sophisticated theories can be used, such as the Hubbard-Onsager theory of dielectric friction, their validity as a quantitative tool to analyse experimental data is questionable. Molecular dynamics provides a much more precise and convincing way to quantify the deviations from the Stokes-Einstein law. A forthcoming study should enable us to relate these quantitative results with existing analytical theories, and to get a clearer view of the relevant physical phenomena to take into account to predict the effective hydrodynamic radius or diffusion coefficient of a small charged nanoparticle.

  As stated in the introduction, we are interested in the flow of water molecules around a specific polyoxometalate ion (POM): the phosphotungstate anion [PW 12 O 40 ] 3-. Counterions are potassium ions K + , and in some cases chloride ions Cl -have also been added. POMs are inorganic anions consisting of a central atom surrounded by MO n polyhedrons, where M is a metal atom and O oxygen. The most known class of polyoxometalates is the Keggin anion that has the general formula [XM 12 O 40 ] n-, depicted in the right part of Fig. 1. In [PW 12 O 40 ] 3-, the central phosphorus is surrounded by oxygens forming a tetrahedron. Each of these oxygen atoms forms the corner of an octahedron, where all corners are formed by oxygen atoms, which enclose the tungsten. Therefore, the surface atoms accessible to the solvent are oxygen atoms shown in various blue shades on Fig. 1, which indicate distinct distances to the central atom. The POM is placed in the center of a cubic simulation box as depicted in the left part of Fig. 1. Periodic boundary conditions are applied in all three directions. The box is filled with water molecules. In NEQ simulations, an external force is applied on water molecules along the x-axis. Several previous theoretical investigations of polyoxometalates focused on their electronic properties using density functional theory (DFT). The optimized DFT structures for the Keggin anions [PW 12 O 40 ] 3-and [SiW 12 O 40 ] 4-well coincide with experimental data as cited in reference.

Figure 1 :

 1 Figure 1: Left: sketch of the simulation box, containing a POM ion in the center of the box. The flow is applied along the x-axis. Right: Structure of the POM ion ([PW 12 O 40 ] 3-). The central atom P is in red; Surrounding O atoms forming a tetrahedron are in yellow; Surface O are in blue, and the W atom is in green. Each W atom is surrounded by an octahedron of O atoms.

  For the same parameters, we define as a reference the fluid velocity v o = | Fx | ξ o , where | F x | is the total external force acting on solvent molecules in the x direction. In what follows, there are about 1700 water molecules around a POM in the simulation box, with an external force on each water molecule of 0.1 × 10 -13 N: v o = 12.31 m.s -1 .

  molecules inside the cylinder keep a constant temperature through heat transfer with the thermostatted region of the box. Indeed, no heating was observed over 20 ns: The mean temperature averaged over time and over the whole simulation box have been computed a posteriori using the simulation trajectory. It is 299.5 K, and thus it agrees perfectly with the target temperature of the T yz thermostat. Moreover, the homogeneity and the time invariance of the temperature inside the simulation box have been checked: The temperature was computed in subregions of the simulation box, within hollow cylinders with the same axis as the thermostat boundaries, over durations of 0.2 ns along the production run. All cylinders keep a temperature of ≈ 300 K. As with PUT, the velocity profile averaged over successive durations of 4 ns fluctuates by 1 m s -1 . As shown in Fig.3, the averaged velocity of water molecules obtained when using the Partial Thermostat is very close to (i) the one obtained with the analytical formula for an array of spherical obstacles, i.e. it is close to v 0 , and (ii) the one obtained with the Profile Unbiased Thermostat. This agreement is a strong evidence showing the ability of both thermostats to regulate temperature in non-equilibrium simulations around a quasi-spherical obstacle.

Figure 3 :

 3 Figure3: Velocity profiles of water molecules along the x-direction obtained with the partial thermostat with a cylindrical frontier and compared to those obtained with the PUT for a resolution of 4 3 bins. The center of the POM is located at x = 0.

  conclusion, we have identified two ways of properly thermostatting the simulation box in non-equilibrium simulations of liquid flow around a spherical obstacle: (i) The Profile Unbiased Thermostat, and (ii) the use of a cylindrical Partial Thermostat. Both methods allow us to obtain the same velocity profiles of water molecules at steady state, and the same density distribution functions around the central nanoparticle. Still, in principle the Nosé Hoover algorithm may affect the local diffusional dynamics of molecules even if the average properties (velocity profiles, distribution functions) are correct. The friction coefficient we are interested in is especially influenced by the structure and dynamics of water molecules close to the central nanoparticle.

Figure 4 :

 4 Figure 4: Distribution functions between the center of mass of water molecules and the phosphorus atom of the POM computed in a cylindrical bin centered on the POM, of length d bin = 0.373 Å, as a function of r yz = y 2 + z 2 . Comparison between results obtained:From equilibrium simulation at 300 K (black line), from equilibrium simulation at 350 K (red dotted line), and from non-equilibrium simulations with a partial thermostat with a cylindrical frontier (blue circles) at 300 K and with a PUT with 4 3 or 8 3 bins.

Figure 5 :

 5 Figure 5: Friction coefficient (top) and hydrodynamic radius (bottom) of the central nanoparticle as a function of its charge, with and without added salt, computed from non-equilibrium simulations.The error bars correspond to the standard uncertainty. 24

5 Figure 6 :

 56 Figure 6: Number of water molecules as function of the distance to the center of the nanoparticle in the systems with no added salt, computed from the radial distribution function between the phosphorus central atom of the nanoparticle and the oxygen atom of the water molecule.

Figure 7 :

 7 Figure 7: Mean velocity of water molecules as function of the distance to the center of the nanoparticle.

Figure 8 :

 8 Figure 8: Residence time of water molecules within a spherical volume of radius around the nanoparticle.

Table 1 :

 1 Force field parameters for the [PW 12 O 40 ] -3 ion and distances of the different atoms to the center of the POM.

	q [e]

i d P -i [ Å] i-i [kcal.mol -1 ] σ i-i [ Å]

Table 2 :

 2 Simulated systems with abbreviations. The box length L box is given in Å. Concentrations c POM , c K + and c Cl -are given in mol.L -1 , q P is the charge on the phosphorus in multiples of the elementary charge e.

	simulation	L box	n H 2 O c POM c K +	c Cl -	q P
	POM0	37.32456 1705 0.032	0.0	0.0	4.51
	POM0-salt 37.35643 1676 0.032 0.478 0.478 4.51
	POM3	37.21562 1702 0.032 0.097	0.0	1.51
	POM3-salt 37.42509 1681 0.032 0.475 0.380 1.51
	POM4	37.35305 1701 0.032 0.127	0.0	0.51
	POM4-salt 37.27600 1681 0.032 0.481 0.353 0.51
	POM5	37.18164 1701 0.032 0.162	0.0	-0.49
	POM5-salt 37.27939 1682 0.032 0.481 0.321 -0.49

Acknowledgement

The work has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 674979-NANOTRANS.