Mael Bosson
email: mael.bosson@inria.fr

Sergei Grudinin
email: sergei.grudinin@inria.fr

Xavier Bouju
email: xavier.bouju@cemes.fr

Stephane Redon
email: stephane.redon@inria.fr

Interactive physically-based structural modeling of hydrocarbon systems

Keywords: Interactive modeling, Adaptive algorithms, Hydrocarbon builder, Brenner potential

Hydrocarbon systems have been intensively studied via numerical methods, including electronic structure computations, molecular dynamics and Monte Carlo simulations. Typically, these methods require an initial structural model (atomic positions and types, topology, etc.) that may be produced using scripts and/or modeling tools. For many systems, however, these building methods may be ineffective, as the user may have to specify the positions of numerous atoms while maintaining structural plausibility.

In this paper, we present an interactive physically-based modeling tool to construct structural models of hydrocarbon systems. As the user edits the geometry of the system, atomic positions are also influenced by the Brenner potential, a well-known bond-order reactive potential. In order to be able to interactively edit systems containing numerous atoms, we introduce a new adaptive simulation algorithm, as well as a novel algorithm to incrementally update the forces and the total potential energy based on the list of updated

Introduction

Scientific and technological advances are increasingly enabling the design of objects at the molecular scale. [START_REF] Ample | A Morse manipulator molecule for the modulation of metallic Shockley surface states[END_REF][START_REF] Dietz | Folding DNA into twisted and curved nanoscale shapes[END_REF][START_REF] Grill | Rolling a single molecular wheel at the atomic scale[END_REF][START_REF] Joachim | The design of a nanoscale molecular barrow[END_REF][START_REF] Shirai | Surface-rolling molecules[END_REF] Similar to the digital prototyping tools that have been developed for macroscopic engineering, powerful tools are more and more in demand to study, design and prototype nanoobjects. The complexity of the underlying physics makes it difficult to manually create and edit virtual nano-objects, though. Hence, if atomic positions are not guided by the digital prototyping tool to plausible configurations, most users will often create unrealistic, high-energy structures that need to be minimized prior to start simulation steps. The final structure might then be very different from the initial plan of the user. An effective way to help the user overcome this problem is to continuously minimize the system's energy during editing.

Such physically-based interactive modeling is a challenging problem, though, since it requires a very low computational cost at each modeling step. Semi-empirical electronic structure calculations and first principle calculations may not be appropriate for such a modeling because of their high computational cost. Moreover, most classical force fields essentially model valence bonds using a harmonic potential, thus preventing topology changes during simulation, which reduces their applicability in interactive builders. An effective approach consists in simulating bond formation and destruction with so-called reactive force fields, such as the ReaxFF [START_REF] Van Duin | ReaxFF: a reactive force field for hydrocarbons[END_REF] method or the Brenner potential.

The Brenner potential is an empirical many-body potential-energy expression with a chemical bond description. [START_REF] Brenner | Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films[END_REF][START_REF] Brenner | The art and science of an analytic potential[END_REF] It has been used to model e.g. carbon nanotube growth [START_REF] Sinnott | Model of carbon nanotube growth through chemical vapor deposition[END_REF] or reconstructed diamond growth [START_REF] Harris | Growth on the reconstructed diamond (100) surface[END_REF] through chemical vapor deposition, nanotube fracture, [START_REF] Belytschko | Atomistic simulations of nanotube fracture[END_REF] buckminsterfullerene dynamics, [START_REF] Brenner | Molecular dynamics simulations of the nanometer-scale mechanical properties of compressed Buckminsterfullerene[END_REF] etc. The first proposed version of the Brenner potential had some limitations that were later removed. [START_REF] Brenner | A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons[END_REF] Several other extensions or corrections have been proposed, e.g. to handle other atom types [START_REF] Dyson | Extension of the Brenner empirical interatomic potential to C-Si-H systems[END_REF] or weak long-range interactions. [START_REF] Los | Intrinsic long-range bond-order potential for carbon: Performance in Monte Carlo simulations of graphitization[END_REF][START_REF] Stuart | A reactive potential for hydrocarbons with intermolecular interactions[END_REF] Brenner et al. have described the use of this potential to simulate nano-engineering, which involves the creation and destruction of chemical bonds. [START_REF] Brenner | Simulated engineering of nanostructures[END_REF] We show in this paper that this potential is also particularly adapted to virtual nano-engineering, i.e. interactive digital modeling of carbon and hydrocarbon structures.

We have recently proposed an approach to speed up a modeling time step in the framework of a classical force field, and demonstrated it on the CHARMM19 and CHARMM22 force fields. [START_REF] Rossi | Adaptive torsion-angle quasi-statics: a general simulation method with applications to protein structure analysis and design[END_REF] The method relies on partial updates of neighbor lists and interaction potentials and forces, which are made possible by incremental updates of the molecular state description (precisely, only a subset of the torsional degrees of freedom are moved at each time step). The other key feature is an efficient algorithm able to predict which degrees of freedom are the most important. We propose a similar approach in cartesian coordinates to dynamically focus the available power on the most active region of the system in editing. Precisely, the algorithm chooses a subset of mobile atoms and, taking advantage of the fact that Brenner's reactive force-field function only depends on atomic relative positions, performs an incremental update of the Brenner potential. This enables interactivity and efficient guide for the user during the editing of a system.

The paper is organized as follows. An overview is proposed in Sect. 2. It provides a brief introduction to our interactive modeling methodology and to the Brenner potential. Sect. 3 analyzes the structure of the computations involved in updating the Brenner potential, which makes an incremental approach possible, and deduces an algorithm to incrementally update the total potential energy and interatomic forces. Sect. 4 describes the new adaptive approach. In Sect. 5 we present the applications and results. The computational cost of the incremental update is analyzed. A new minimizer is proposed as an application of the new adaptive approach. This leads to the interactive modeler that may be used for digital prototyping of hydrocarbon systems. Finally, Sect. 6 concludes the paper and discusses possible extensions and other applications of our work.

Overview

Interactive physically-based structural modeling

The main idea behind interactive physically-based structural modeling is to continuously minimize the system energy during editing to help the user build realistic structures. In practice, our builder alternates a time step of minimization with a time step of user action. In order to maintain interactive update rates even when the model contains numerous atoms, though, we propose to use an adaptive, partial update of the system's state at each time step.

Precisely, the structure of interactive modeling algorithm time step consists in alternating two main steps: When the user performs a local action on the system of atoms, the potential energy and interatomic forces are incrementally updated: thanks to the structure of the Brenner potential, only some forces need to be updated, and the total potential energy may be updated through recomputation of some terms only. This helps us to control the computational cost of a time step.

We describe in Sect. 3 a new algorithm to incrementally update energies and forces.

Once all forces and the potential energy are known, the adaptive minimization step chooses a set of "most mobile atoms" and updates the positions of these atoms only (since moving all atoms might prevent interactive rates when constructing large systems). The choice of the set of active degrees of freedom (i.e. the set of mobile atoms) is performed automatically at each time step, based on user-defined precision constraints (e.g. the total number of mobile atoms at each time step, or the maximum allowable error on the acceleration of the molecular system), as well as the current forces being applied on each atom. Controlling the number of mobile atoms, and thus the cost of the system's update, allows us to decouple the cost of the time step from the number of atoms in the system, and helps us maintain interactive editing rates. We describe in Sect. 4 a new adaptive algorithm to determine the set of most mobile atoms.

Note that detecting neighboring atoms is independent of the underlying force-field as well, and we may use for example the adaptive neighbor list algorithm described in Rossi et al.. [START_REF] Rossi | Adaptive torsion-angle quasi-statics: a general simulation method with applications to protein structure analysis and design[END_REF] As a result, we assume in this paper that, when we need to update forces and energies at each time step, we are given the list of pairs of neighboring atoms whose relative position has changed due to user actions or to adaptive minimization.

Brenner potential

For completeness, we now introduce the second generation Brenner potential. [START_REF] Brenner | A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons[END_REF] This potential belongs to the class of bond-order potentials, and thus has the ability to evolve smoothly between different bonding states. Shortly, it al-lows the calculation of the interaction energy of chemically covalent bond parametrized by a bond-order function.

Potential expression

The Brenner potential has the particularity to be only focused on covalent bonds (no long-range interaction), i.e. the total potential energy is a sum over interacting atoms (two atoms are said to be interacting if the distance between them is smaller than 2 Å):

E b = i j>i [V R (r ij ) -b ij V A (r ij )] , (1)
where b ij is the bond order. The functions V R and V A are respectively repulsive and attractive terms. A complete description of these terms and their parameters is provided in the original publication [START_REF] Brenner | A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons[END_REF]. The complete analytic form of the function is recalled in appendix.

Forces

In order to proceed, we need to determine the force terms, i.e. the gradient of the potential, from the description above. Because the force applied on one atom only depends on relative positions, we can express it as a sum over the interacting pairs in which this atom is involved:

∂V dx i = j,(i,j)∈B ∂V dr ji ∂r ji dx i , (2)
where r ji is the distance between atoms i and j, and B is a set of all the pairs of atoms involved in interaction:

B = (i, j), r ij < D max ij , (3)
where D max ij depends on atom types.

For future use, we also define:

A = {(i, j, k), (i, j) and (j, k) ∈ B, i = j} , (4)
D = {(i, j, k, l), (i, j), (j, k) and (k, l) ∈ B, i = k = l} , (5)
and we decompose the potential as follows:

V = i j>i [V R (r ij ) -b ij V A (r ij )] = (i,j)∈B, j>i V B ij . (6)
Then, a direct analysis of the influence of the relative position of every pair in interaction on the potential is sufficient to know the forces on each atom. Because of the bond order term, the relative position of atoms i and j can influence several V B kl terms. Indeed, the total force on one atom can be split into three contributions.

∂V dx i = j,(i,j)∈B   (k,l)∈B ∂V B kl dr ji   ∂r ji dx i , (7)
(k,l)∈B ∂V B kl dr ji = ∂V B ij dr ji + ∂V A ij dr ji + ∂V D ij dr ji , (8)
where V A ji is the sum of bond terms that involve either i or j:

V A ij = k,(j,i,k)∈A V B ik + k,(i,j,k)∈A V B jk , (9)
and V D ji is the sum of bond terms that involve neighbors of i or j:

V D ij = (k,l),(i,j,k,l)∈D V B kl + (k,l),(j,i,k,l)∈D V B kl . (10)
3. An incremental algorithm to update forces and potential

We may now deduce an algorithm able to incrementally update both the potential and atomic forces.

Influence of relative motions

In the previous section, we have shown that relative positions of interacting atoms may influence several terms in the potential. Figure 1 presents these dependencies in a schematic case. Adaptive modeling also requires us to know the dependency of forces on the relative position of two atoms, but this may be directly deduced from the previous analysis. Indeed, equation 7 shows that the force on atom i needs to be updated when either one of the three terms in the sum in equation 8 has been updated. These terms, in turn, need to be updated when the bond terms involved in equation 9 and equation 10 are updated. This dependency is illustrated in Figure 2. Figure 3 shows the complete dependency between the terms involved in both the potential and the atomic forces.

The algorithm

As noted in Sect. 2, the incremental collision detection algorithm provides us with a set of pairs of atoms whose relative position has changed since the last time step. Let B updated denote this set. Note that this set may include pairs of newly interacting atoms, as well as pairs of atoms which are now separated by a distance longer than the force-field cutoff (2 Å in the case of the Brenner potential), but which were interacting at the previous time step.

At each time step, the algorithm starts with this information to update the sets B, A and D. This update is incremental in the sense that the sets are not rebuild from scratch, but elements are inserted or deleted. In practice, the algorithm stores all terms which contribute to the total potential and to the atomic forces. To achieve this, we introduce four sets of tuples:

I Atom = I Atom i = (N H i , N C i , N t i , N aux i ), i ∈ Atoms , I Bond = I Bond ij = (V R ij , V A ij , b σ-π ij , S σ-π ij , b σ-π ji , S σ-π ji , Π RC ij , b DH ji , S DH ji , N conj ij ) , i < j, (i, j) ∈ B} , I Angle = I Angle ijk = (f ik (r ik )G(cos(θ ijk ))e λ ijk ), i < k, (i, j, k) ∈ A , I Dihedral = I Dihedral ijkl = {(1 -cos(Θ ijkl ))f ik (r ik )f jl (r jl )}, i < l, (i, j, k, l) ∈ D . (11)
These sets contain the terms that need to be updated to compute the total potential for the adaptive simulation. They contain respectively the information attached to one, two, three or four atoms (please refer to the Appendix for a detailed description of the potential function).

Similarly, we define four more sets to allow for an incremental update of the forces, corresponding to the decomposition proposed in the previous section:

F atom = F atom i = ∂V dr i , i ∈ Atoms , F B = F B ij = ∂V B ij dr ji , i < j, (i, j) ∈ B . F A = F A ijk = ∂V B ij dr ik , i < k, (i, j, k) ∈ A , F D = F D ijkl = ∂V B ij dr kl , i < l, (i, j, k, l) ∈ D . (12)
F atom contains the forces acting on each atom, and is updated based on the intermediate force terms contained in the sets F B , F A , and F D .

Pseudo-code

The structure of the algorithm is deduced from the hierarchical dependency presented in Figure 3. The algorithm works with auxiliary sets Atoms, B, A and D containing respectively singletons, pairs, triplets and quadruplets of atomic indices which correspond to obsolete terms. The computation of the potential and forces is organized in four main steps:

1. Initialization.

Based on B updated , the output of the collision detection algorithm, incrementally update the sets B, A and D as indicated in Sect. 3.2. Set Atoms, B, A and D to ∅.

First level potential update.

Detect terms that have to be updated directly from the list of bonds with a relative motion, i.e. ∀(i, j) ∈ B updated : Extend the list of potential terms that have to be updated to include bonds without a relative motion, but with a bond order change:

(a) insert (i, j) in B (b) insert i and j in Atoms (c) ∀k, (i, j, k) ∈ A, insert (i, j, k) in A (d) ∀(k, l), (i, j, k, l) ∈ D, insert (i, j, k, l) in D (e) ∀(k, l), (k, i, j, l) ∈ D, k < l, insert (k, i, j, l) in D Update terms: (a) ∀ i ∈ Atoms, update the four corresponding values in I Atom i (b) ∀(i, j) ∈ B, update V R ij , V A ij and N conj ij in I Bond ij (c) ∀(i, j, k) ∈ A,
(a) ∀(i, j) ∈ B updated , ∀(k), (i, j, k) ∈ A, insert (j, k) in B (b) ∀(i, j) ∈ B updated , ∀(k, l), (i, j, k, l) ∈ D, insert (k, l) in B Update the total potential (a) ∀(i, j) ∈ B, update all the bond order terms (i.e. b σ-π ij , b σ-π ji ,Π RC ij and b DH ji ). Update V B ij based on V R ij , V A ij
and b ij , and update its contribution to the total sum V .

Force update.

Extend invalid sets to invalid force terms: We do not provide the details of the update of the F B , F A and F D terms, as their expression is straightforward -it can be directly retrieved form the expression of the potential -but tedious.

∀(i, j) ∈ B, (a) ∀k, (i, j, k) ∈ A, insert (i, j, k) in A (b) ∀(k, l), (i, j, k, l) ∈ D, insert (i, j, k, l) in D
Note on sums updates: one important point in our algorithm is how we update the sums involved in the Brenner potential and forces. Indeed, sums are never recomputed from scratch, but are incrementally updated. Precisely, when we update a term and have to update its contribution to a sum (e.g. one term V B kl has been recomputed, and the sum in Equation 10has to be updated), we subtract the old value of the term from the sum, and add the new value to it. This allows us to update obsolete terms only. Numerical error may be increased by this approach, however, due to the accumulation of subtractions and additions. To avoid this issue, a periodic re-initialization may be performed every e.g. million steps, where energies and forces are computed for the whole molecular system. Periodic re-initialization is very rarely needed, though, since the sums involve only a limited number of terms, and the loss of precision due to one incremental sum update (one addition and one subtraction) is around machine precision. Thus, the computational cost is amortized over the whole interactive session.

Extension of the Brenner potential

For some applications, a notable limitation of the Brenner potential is the absence of long-range interactions. To overcome this problem, some extensions have been proposed. [START_REF] Los | Monte Carlo simulations of carbon-based structures based on an extended Brenner potential[END_REF][START_REF] Los | Improved longrange reactive bond-order potential for carbon. I. Construction[END_REF][START_REF] Ghiringhelli | Improved longrange reactive bond-order potential for carbon. II. Molecular simulation of liquid carbon[END_REF][START_REF] Karssemeijer | Phonons of graphene and graphitic materials derived from the empirical potential LCBOPII[END_REF][START_REF] Stuart | A reactive potential for hydrocarbons with intermolecular interactions[END_REF] The general idea is to add a non-bonded interaction potential term that is a sum of pairwise potential contributions:

E N B = 1 2 i,j V N B (r ij ). ( 13
)
The difficulty is to correctly treat the intermolecular interactions without biasing the intramolecular part. This correction has been done carefully in e.g. the popular AIREBO potential [START_REF] Stuart | A reactive potential for hydrocarbons with intermolecular interactions[END_REF]. In this potential, the long-range term is progressively activated and combined with the accurate short-range part of the Brenner potential.

Long-range interactions are particularly important when one deals with carbon structures above a surface. The adsorption of fullerenes and nanotubes has been extensively studied with a large range of surfaces, including metallic, semiconducting and insulating ones. Here, we focus on the graphite surface and graphene as well. A previous study was based on a self-consistent real-space scheme for calculating the van der Waals interaction energy between a fullerene molecule and a graphitic substrate with atomic surface corrugation. [START_REF] Gravil | Adsorption of C 60 molecules[END_REF] By using linear response theory to evaluate the dipole-dipole interactions between the molecule and the substrate, the van der Waals energy was accurately computed and the optimized conformations of the adsorbed molecule were found.

For simplicity, we use the approach of Los and Fasolino [START_REF] Los | Monte Carlo simulations of carbon-based structures based on an extended Brenner potential[END_REF], and we modulate the exp-6 van der Waals potential parameters:

V N B(r ij ) = b exp(-c 0 r) -ǫ σ r 6 -V shif t , (14)
to retrieve optimum adsorption energies close to the already calculated ones and the corresponding distances between the center of the molecule and the surface. [START_REF] Gravil | Adsorption of C 60 molecules[END_REF] Precisely, with parameters b = 3224.9 eV, c 0 = 3.5995 Å-1 , ǫ = 0.01396 eV, σ = 3.44 Å, the energy errors are 0.0162eV and 0.0108eV and the distance errors are 0.08 Å and 0.09 Å for the top and hollow configurations respectively. In this case, the two-body potential expression can be introduced efficiently in our adaptive approach thanks to the notion of partial forces. [START_REF] Rossi | Adaptive torsion-angle quasi-statics: a general simulation method with applications to protein structure analysis and design[END_REF] As illustrated in Figure 4, the preferred orientation of the fullerene is always with a six-membered ring parallel to the surface. The agreement with other calculations and experiments is very satisfactory. [START_REF] Dappe | Weak chemical interaction and van der Waals forces between graphene layers: A combined density functional and intermolecular perturbation theory approach[END_REF][START_REF] Dappe | Intermolecular interaction in density functional theory: Application to carbon nanotubes and fullerenes[END_REF][START_REF] Seydou | Atomic force microscope measurements and LCAO-S 2 + vdW calculations of contact length between carbon nanotube and graphene surface[END_REF] This stacking is somehow similar to the adsorption of aromatic molecules on graphene, where π-π interactions dominate. [START_REF] Dappe | Weak chemical interaction and van der Waals forces between graphene layers: A combined density functional and intermolecular perturbation theory approach[END_REF][START_REF] Dappe | Intermolecular interaction in density functional theory: Application to carbon nanotubes and fullerenes[END_REF][START_REF] Shtogun | Many-body van der Waals interaction between graphitic nanostructures[END_REF][START_REF] Björk | Adsorption of aromatic and anti-aromatic systems on graphene through π-π stacking[END_REF] Once a confident set of parameters is extracted, one can easily deal with larger systems where long-range interaction are important. For instance, it is well known that large fullerenes and nanotubes are deformed when adsorbed on a surface. [START_REF] Hertel | Deformation of carbon nanotubes by surface van der Waals forces[END_REF] One can handle such phenomena with this new parametrization. The approach we propose in this paper should also deserve attention to investigate structural properties of graphite [START_REF] Lebègue | Cohesive properties and asymptotics of the dispersion interaction in graphite by the random phase approximation[END_REF] or graphene sheets in which one is able to design structures like ribbons or Y-shaped junctions for nano-electronics.

Finally, the interaction of a carbon nanotube with a graphene surface is an important topic in the context of ultimate atomic force microscopy (AFM). [START_REF] Wilson | Carbon nanotube tips for atomic force microscopy[END_REF] Indeed, it has been shown that a single-walled carbon nanotube anchored at the apex of an AFM tip increases the sensitivity of the measured forces. By cycling approaches and retracts, one is able to extract physical information of the surface, such as energy dissipation. [START_REF] Seydou | Molecular mechanics investigations of carbon nanotube and graphene sheet interaction[END_REF][START_REF] Seydou | Atomic force microscope measurements and LCAO-S 2 + vdW calculations of contact length between carbon nanotube and graphene surface[END_REF][START_REF] Delmas | Contact angle hysteresis at the nanometer scale[END_REF] Nevertheless, it is crucial to understand the mechanical behavior of the nanotube, by mixing strong deformations when it vertically approaches the surface with long-range van der Waals contributions when it is deposited on it. We believe that all these developments can be addressed by the on-the-fly calculations using the algorithm and the numerical tool described in this paper.

A new adaptive algorithm

The incremental algorithm presented in the previous section is able to take advantage of the freezing of groups of atoms. We now focus on the choice of the active degrees of freedom. An adaptive algorithm has been developed in the framework of internal coordinates and articulated bodies simulation [START_REF] Redon | Adaptive dynamics of articulated bodies[END_REF]. Here, we extend this approach to cartesian coordinates in order to address the problem of providing interactive and efficient structure optimization. This is the essential feature that guides the user in our modeler.

The adaptive approach in cartesian coordinates is to decide for each atom if it might move or be frozen in space. Thus, the algorithm has to decide to activate or freeze an atom depending on the norm of the force applied to the atom. This decision is made by comparing the norm with a certain threshold value. This threshold value is defined either automatically based on the system state (e.g. average acceleration, maximum acceleration) or can be simply predefined by the user. This cartesian adaptive approach is illustrated in Figure 5. In this paper, we implement the cartesian adaptive approach using the steepest descent minimization algorithm, [START_REF] Press | Numerical recipes in FORTRAN: The art of scientific computing[END_REF] which serves to guide the model towards an optimized geometry. In practice, the adaptive minimizer does not use the components of the energy gradient (the partial derivatives) that are below a given threshold. This amounts to providing a perturbed minimization direction. This is detailed in Sect. 5 as an application of this new adaptive approach coupled with the incremental update of the Brenner potential.

An important feature in our algorithm is that the threshold criteria and the active atoms set are not computed from scratch. At the end of each point energy calculation a list of the atoms with force change and the corresponding new threshold are incrementally computed as explained in Sect. 3. In order to avoid the linear cost of determining the new set of active atoms, a binary tree is used to represent the system. Each leaf node represents an individual atom, while each internal node represents a set of atoms. Each leaf node stores the norm of the force applied to the corresponding atom. Each nonleaf node stores the maximum of the two force norms of its children, as illustrated in Figure 6. We use two tree passes in order to update tree nodes' values and to determine the new active atoms. In the first, bottom-up pass, force norms are updated in a sub-tree of the binary tree, starting from the leaves with modified norms, in O(k old (log( n k old ) + 1)) times where k old is the number of active atoms and n the total number of atoms. In the second, top-down pass, the new active atoms (i.e. the atoms with the force norms which are now the largest), are determined in O(k new (log( n k new ) + 1)) times

where k new is the new number of active atoms. This process is illustrated in

Results and applications

Computational complexity of the Brenner potential incremental update

Our algorithm is linear with the size of the working sets Atoms, B, A and D. Based on the assumption that each atom may only have a bounded number of neighbors (which is true for realistic systems), these sets have a linear size in the number of bond updates |B updated |. Thus, the overall complexity is linear in the number of relative positions updates indexed in B updated , since each relative position update has an impact on a bounded number of terms.

In practice, the complexity of the algorithm depends on the local topology of the system and the bond update pattern. To demonstrate this, we generated a carbon nanotube of 8000 atoms thanks to TubeGen, [START_REF] Frey | Tubegen 3.3[END_REF] and perturbed the positions of a set of N atoms, 1 < N < 8000. This benchmark is presented in Figure 7 where the solid curves represent the computational cost of the forces and potential update. To illustrate the dependency of this cost on the activation pattern (i.e. the evolving set of active atoms), both random and continuous activations are compared. The curve with a continuous activation pattern shows a linear behavior. Indeed, continuous activation is equivalent to consider a nanotube system of length equal to the length of the active region plus a buffer zone. Thus, the linear behavior is a direct consequence of the linear dependency of the Brenner potential on the number of bonds. The computational cost of the random activation step is more important, since updated relative positions are isolated in space with a high probability, which results in a larger number of terms that have to be recomputed. In general, for any molecular system, there exists two constants C min and C max such that the total time T U to update the potential and forces after updating N b relative positions satisfies:

C min N b ≤ T U ≤ C max N b . ( 15
)
Figure 7 shows the update time T U as a function of the number of bonds with a relative motion.

! "!!! #!!! $!!! %!!! &!!!! &"!!! ! '! &!! &'! ()*

Efficient and interactive energy minimization of local deformation

Energy minimization is a classical tool for structure analysis and has already been used with the Brenner potential. [START_REF] Hod | Carbon nanotube closed-ring structures[END_REF][START_REF] Huang | Thickness of graphene and single-wall carbon nanotubes[END_REF][START_REF] Reddy | Equilibrium configuration and continuum elastic properties of finite sized graphene[END_REF] We apply the adaptive algorithm presented in Sect. 4 with the incremental update to the steepest descent minimization. Classical minimization algorithms typically update the positions of all atoms at each iteration. As a result, all energies and forces are updated at each iteration. In many cases, however, the displacement of many atoms will be small, and will not contribute much to a decrease in potential energy. We thus propose to focus the computational resources by freezing in space the least mobile atoms at each iteration. As discussed in the previous section, our incremental algorithm is able to take advantage of the freezing of groups of atoms. As a result, computational resources will be automatically focused on the most mobile areas of the system. This approach is expected to speed up energy minimization of systems with non-uniform force distributions.

We demonstrate this approach on a locally deformed carbon nanotube with a few broken cycles, which should converge to the reference structure (Figure 8).

We minimize the deformed structure using an automatic, relative threshold (only the atoms that have an acceleration larger than half the maximum atom acceleration at the current time step will move). We stop the minimization procedure when the root-mean-square-deviation (RMSD) [START_REF] Theobald | Rapid calculation of RMSDs using a quaternion-based characteristic polynomial[END_REF] between the structure being minimized and the reference one becomes smaller than 0.01 Å. The Brenner potential accuracy is in the order of 0.01 Å, [START_REF] Petukhov | Reconstructions of diamond (100) and (111) surfaces: Accuracy of the Brenner potential[END_REF] and thus a better convergence would not be significant. Figure 9 plots the RMSD as a function of wall-clock time. The adaptive minimizer converges 3.5 times faster than the original minimizer to the RMSD of 0.1 Å and more than two times faster to the RMSD of 0.01 Å. Indeed, in a first period, the deformation is very local and an adaptive approach is very efficient. Then, when the force distribution is more uniform, the non-adaptive approach becomes faster due to the overhead of the adaptive algorithm. This is illustrated by the second curve in Figure 9, which shows the performance of the adaptive steepest descent when all degrees of freedom are active. In this case, detecting which terms have to be updated is useless and time consuming. In our implementation, the overhead is around 6 %. Figure 10 presents the snapshot of a video showing the computational benefits of the incremental approach in a similar situation of a nanotube locally deformed with a few broken cycles.

The video is accessible on the web: link.

To understand why our algorithm can accelerate energy minimization, recall that the adaptive steepest descent algorithm does not follow the classical gradient direction on the energy surface but a close, truncated one, to reduce the cost of a minimization step. In fact, this is equivalent to looking for a more efficient direction in the sense of the ratio between the energy decrease and the effective computational cost of the position update. The classical steepest descent algorithm is based on the idea that the most efficient way to reduce the energy is to follow the steepest direction on the energy surface.

However, in this argument, the update of the energy potential is seen as a black box and thus the algorithm does not consider the (potentially variable) update cost. We believe this update cost should be taken into account in order to choose the minimization direction, and that it might be beneficial to try and maximize the energy decrease per computational unit at each time step. The overall benefit of such a strategy may depend on the general shape of the potential energy surface. Therefore, it is still difficult to predict the advantage of such an approach in the general case. We have shown that our adaptive approach allows us to finely tune between the update cost and the energy decrease, which speeds up the minimization when the force distribution is localized (when the force distribution is uniform, moving all atoms whose force is larger than half the maximum force will lead to activating all atoms, so that our truncated gradient will be equal to the actual gradient).

Interactive modeling

The classical approach to study hydrocarbon systems is typically to first build a structural model, and then to perform geometry optimization to prepare the model for simulation. Because of the complexity of the physics at this scale, however, this approach is not efficient in many cases. Indeed, when the user creates the initial structural model, he can easily produce high-energy, unstable structures, in particular for large systems. As a result, the geometry optimization step may bring on very different structures from the ones targeted by the user, so that the user would have to go back to the construction step. Overall, creating a complex structural model by decoupling the construction step from the geometry optimization step may be particularly ineffective.

Thus, we believe an effective way to help the user to construct a structural model of an hydrocarbon system is to perform geometry optimization during construction. In practice, at each time step, we combine user actions with energy minimization to produce structures that have lower potential energies. This physically-based interactive modeling requires a very low computational cost at each simulation step. One possibility for the user is to act on a small area of the system, creating local deformations only. Another possibility is to perform deformations which require physics simulation in local areas only (such as when two large groups of atoms move rigidly and interact with each other through a few atoms only). The minimizer presented in the previous section has been designed for such cases, since the force distribution should reach a peak in a local area. Furthermore, the adaptive minimizer can guarantee interactivity independently of the size of the system (by directly controlling the threshold value and thus the number of allowed components in the truncated gradient -a choice specifically based on available computational resources). We also recall that the reactive Brenner potential is particularly adapted, since a fundamental step in structural modeling is the creation and destruction of chemical bonds. This makes coupling an adaptive minimization algorithm with an incremental reactive potential update algorithm an effective way to achieve interactive modeling.

Remark that a local deformation can have a local impact as well as a global impact. The locality of the active zone does not prevent the minimizer to produce delocalized deformation. Indeed in the spirit of the adaptive paradigm, the active region is automatically determined thanks to the information of the atomic forces and it might go through the whole system in a linear number of time steps when needed. An important benefit is that a local user action might lead to a non-local, coordinated displacement of a group of atoms (for example when a user pulls on an atom that belongs to a chain), which enables the user to induce large, physically-based deformations of the nanosystem. More generally, such an approach allows the user to perform even non-local editing (e.g. rotation and translation of entire groups of atoms), while ensuring that the underlying physics is taken into account.

We illustrate this possibility in Figure 14, where a complex, original nanosystem is built in a few minutes thanks to the interactivity of the physically-based modeler. In Step 0, the user loads two equilibrated graphene sheets in the modeling tool. In Step 1, the user connects both sheets (left part in Figure 14), simply by pulling on atoms. As the two sheets are being attached, the minimizer progressively induces a curvature in the structure, which actually helps the user to attach the sheets. In Step 2, the right part is connected, again by pulling on atoms. In the final step, the two remaining sides have been attached, and a few atoms has been removed from the corners to form a "nano-vesicle". A closed bilayer graphene has already been observed experimentally by high-resolution transmission electron microscopy [START_REF] Liu | Open and Closed Edges of Graphene Layers[END_REF]. This shows that the present numerical tool should be efficient for experimentalists to design and to build graphitic nanostructures. Similarly, carbon nanotubes with a large variety of defects have already been observed experimentally by high-resolution transmission electron microscopy [START_REF] Chamberlain | Reactions of the inner surface of carbon nanotubes and nanoprotrusion processes imaged at the atomic scale[END_REF] and our methodology allows to intuitively prototype such complex geometries. Figure 11 shows how our modeler may be used to interactively close a nanotube. The user has chosen to add atoms and put them close together to create cycles of 5, 6 or 7 carbon atoms. The tool automatically induces the corresponding curvature on the system leading to the design of the apex of the nanotube. As discussed above, the process could have been done with a nanotube of any size as the computational cost of each minimization step only depends on the number of active atoms. We also present two videos to illustrate the efficiency of the interactive editor.

Figure 12 presents a snapshot of a video in which a user designs the apex of a nanotube. Figure 13 presents a snapshot of a video in which the user adds structural defects to a carbon nanotube model. The videos are also accessible on the web: link.

Conclusion

In this paper, we have presented two algorithms enabling digital prototyping of hydrocarbon structures. Interactivity is guaranteed by the adaptive minimization algorithm coupled with the incremental update algorithm of the Brenner potential. This allows us to focus computational resources on the regions that are the most affected by user actions. Furthermore, when user actions have a local impact, the adaptive approach appears to be an effective way to rapidly reach neighboring energy minima, which helps the user build realistic structures. The proposed modified steepest descent algorithm may allow for an important speed up when the energy gradient is non-uniform. This happens when the user stresses the system locally. This might happen in other contexts as well (e.g. docking problems). We would like to extend our approach in several directions. First we have mentioned that the overhead resulting from the incremental update and the marginal cost of a relative motion may make our approach slower than the classical one in some cases. This suggests the need for hybrid minimization algorithms which would be able to switch between the classical and the adaptive approach at runtime based on the distribution of forces in the system. Also, we would like to investigate whether our new cartesian adaptive approach to compute a truncated gradient could be useful in combination with other minimization algorithms.

!"#$%& !"#$%' !"#$%( !"#$%) *+,-.%/"#$%0%!+1#%2+#3 *+,-.%/"#$%0%45$%2+#3
Finally we want to determine whether the cartesian adaptive approach may be combined with an adaptive internal degrees-of-freedom simulation method [START_REF] Rossi | Adaptive torsion-angle quasi-statics: a general simulation method with applications to protein structure analysis and design[END_REF], i.e. group atoms into clusters which might deform while creating and breaking bonds, which should provide great functionality both for simulation and digital prototyping.

Despite the computational efficiency of the Brenner potential, the cost of molecular dynamics or Monte Carlo simulations may still be too high when the size of the system increases, or when long simulations are required. A classical approach to address this problems is to use parallel computers. [START_REF] Rosenblum | Multi-processor molecular dynamics using the Brenner potential:. parallelization of an implicit multi-body potential[END_REF][START_REF] Omeltchenko | Crack front propagation and fracture in a graphite sheet: A molecular-dynamics study on parallel computers[END_REF][START_REF] Srivastava | Molecular dynamics simulation of largescale carbon nanotubes on a shared-memory architecture[END_REF] We believe that we may extend our work in the future as an alternative answer. Indeed, our incremental algorithm will make it possible to recompute changing terms only, based on the list of moving atoms, and can thus be used in any Monte Carlo methods based on the Brenner potential [START_REF] Los | Monte Carlo simulations of carbon-based structures based on an extended Brenner potential[END_REF][START_REF] Pikunic | Structural modeling of porous carbons: Constrained reverse Monte Carlo method[END_REF][START_REF] Los | Intrinsic long-range bond-order potential for carbon: Performance in Monte Carlo simulations of graphitization[END_REF].

We are also exploring the possibility of adaptive dynamic simulation, with non-zero temperature. Indeed, the proposed tool helps the user produce molecular models that have optimized geometries. However, these structures should thus be analyzed through simulations, to compute e.g. free energy barriers, statistical properties, etc.

Finally, another important challenge in interactive molecular modeling that we would now like to address is to extend the adaptive approach to semiempirical quantum mechanical models, which would provide better accuracy as well as interactive electronic structure computation.

where D max ij and D min ij are type-dependent parameters, and r is the distance between atoms i and j. Originally, the parameters have been determined for carbon and hydrogen only [START_REF] Brenner | Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films[END_REF]. Later they have been extended to Silicon [START_REF] Dyson | Extension of the Brenner empirical interatomic potential to C-Si-H systems[END_REF][START_REF] Que | Application of the extended Brenner potential to the Si(111)7×7:H system I: cluster calculations[END_REF][START_REF] Que | Application of the extended Brenner potential to the Si(111)7×7:H system II: periodic calculations[END_REF][START_REF] Schall | Elastic constants of silicon materials calculated as a function of temperature using a parametrization of the second-generation reactive empirical bond-order potential[END_REF]. Bonds are defined dynamically via the f function. As a result, the potential has the ability to describe chemical reactions.

Bond order

The benefit of a reactive potential mainly depends on the bond order. In the Brenner potential case, we can decompose it into four terms depending on environment parameters including angles, dihedral angles and number of neighbors.

b ij = 1 2 [b σ-π ij + b σ-π ji ] + Π RC ij + b DH ji (20)

Angle influence

The first term in Equation 20 is given as a sum over the angle θ ijk in which the bond (i, j) is involved.

b σ-π ij = 1 + k =i,j
f ik (r ik )G(cos(θ ijk ))e λ ijk + P ij (N C i , N H i )

-1/2 (21)
where G is a spline and P a bicubic spline. λ ijk is a constant. N C i and N H i are respectively the number of carbon and hydrogen neighbors. The last numbers are continuous values as they are computed by a sum of the f function.

N C i = carbon atoms k =i,j f ik (r ik )

N h i = hydrogen atoms k =i,j f ik (r ik ) (22)
For the description of our incremental update algorithm we define :

S σ-π ij = k =i,j
f ik (r ik )G(cos(θ ijk ))e λ ijk [START_REF] Gravil | Adsorption of C 60 molecules[END_REF]

1 .

 1 User action step (a) let the user move / create / delete one or more atoms (b) incrementally update forces and the potential energy 2. Adaptive minimization step (a) choose a new set of mobile atoms (b) perform a minimization step with the corresponding set (c) incrementally update forces and the potential energy

!Figure 1 :

 1 Figure 1: Dependency of the Brenner potential terms on the relative motion of atoms (see Sect. 3.1). Atoms belonging to "Type 1" bonds have a relative motion, which induces a modification in every term of the bond potential contribution. "Type 2" bonds are rigid, but with a bond order change due to its dependency on the number of neighbors, as well as on angles and on dihedrals. "Type 3" bonds are rigid, but with a bond order change caused by the change on the conjugate number N conj . The relative motion of atoms in "Type 1" bonds does not impact the potential contribution of "Type 4" bonds.

Figure 2 : 5 ,Figure 3 :

 253 Figure 2: Dependency of atomic forces on relative positions updates (see Sect. 3.1).

3 .

 3 update the angle-dependent value I Angle ijk and its influence in I Bond ij and I Bond jk (i.e. the terms S σ-π ij and S σ-π ji ) (d) ∀(i, j, k, l) ∈ D, update the dihedral-dependent value I Dihedral ijkl and its influence in I Bond jk (i.e. the term S DH ji ) Second level potential update.

  Update forces: (a) ∀(i, j) ∈ B, update the corresponding F B terms and their influence in the sum defining F Atom i and F Atom j (b) ∀(i, j, k) ∈ A, update the corresponding F A terms and their influence in the sum defining F Atom i , F Atom j and F Atom k (c) ∀(i, j, k, l) ∈ D, update the corresponding F D terms and their influence in the sum defining F Atom i , F Atom j , F Atom k and F Atom l

Figure 4 :

 4 Figure 4: Buckminsterfullerene in interaction with a graphite surface thanks to the fitted pairwise van der Waals term. The fitting procedure has been performed to respect the fullerene-surface stable distance.

Figure 6 Figure 5 :

 65 Figure 6 as well.

20

 20

Figure 6 :

 6 Figure 6: This figure illustrates the procedure to determine the active zone, when the threshold is automatically set to half the largest atomic force norm. The algorithm uses abinary tree to represent the atomic system, in which each leaf corresponds to one atom. In this example, the four leaves correspond, from left to right to atoms 1 to 4. Internal tree nodes represent groups of atoms. The value indicated in each leaf node is the norm of the force applied to its corresponding atom. For internal nodes, this value is the maximum of the norms of the forces applied to atoms in the corresponding group. In step 0, the threshold is automatically set to 10. As a result, only atom 1 moves. In step 1, the potential is incrementally updated, and the norms of the forces applied to atoms 1 and 2 are updated. In step 2, the values associated to the tree nodes are incrementally updated through a bottom-up pass that starts from the modified leaf nodes values. Because of this bottom-up update, the adaptive threshold becomes equal to 4. In step 3, the new active atoms are determined through a top-down pass, by visiting only the nodes that have a value larger than the adaptive threshold.

Figure 7 :

 7 Figure7: Forces and potential incremental update complexity of a carbon nanotube composed of 8000 atoms and 11980 bonds. In this benchmark, a number of pairs of atoms have a relative motion and the solid curves represent the time cost of the forces and potential update. The curve of the random activation pattern is steep at the beginning as the relative motions are isolated in space with a high probability and cause a maximum number of terms to recompute. Here, the slope is similar to the slope C max . For a large number or relative motions, the curve becomes flat since a new relative motion does not add many new terms to recompute. Thus the marginal cost of a new relative motion is minimal. The curve with a continuous activation pattern illustrates the linear dependency of the Brenner potential on the number of bonds, since this continuous activation is equivalent to consider a nanotube system with an increasingly longer length. Its slope is also approximately equal to the slope of the minimal cost C min .

Figure 8 : 4 Figure 9 :

 849 Figure 8: Different snapshots of an experiment where a carbon nanotube is being relaxed after a local deformation. The adaptive steepest descent minimizer automatically focuses the computational resources on the most mobile atoms. This leads to a more efficient minimization scheme in such a situation.

Figure 10 :

 10 Figure 10: Snapshot of a video showing the computational benefits of the incremental approach when minimizing large structures that have been locally deformed.

Figure 11 :

 11 Figure 11: Different steps to prototype a "nano-vesicle" with the adaptive interactive modeler. Starting with two graphene sheets, the adaptive methodology enables the user to prototype a new structure by applying forces on atoms with the mouse. We demonstrate the efficiency of the modeler by the creation of this new original carbon structure.

Figure 12 :

 12 Figure 12: Snapshots of a nanotube capping process with the adaptive interactive modeler.A classical nanotube is loaded in the modeler and the closure is performed by applying forces on bordering atoms. The curvature at the tip appears automatically thanks to the adaptive minimization running while the user edits the structure's topology. Thanks to the adaptive methodology, this operation can be done in a few minutes.

Figure 13 :

 13 Figure 13: Snapshot of a video showing a user interactively prototyping the closure of a carbon nanotube.

Figure 14 :

 14 Figure 14: Snapshot of a video demonstrating how our algorithm may be used to interactively add structural defects to a carbon nanotube model.

Software availability

The new algorithms will be made available through SAMSON (Software for Adaptive Modeling and Simulation of Nanosystems) [START_REF]NANO-D[END_REF], being developed in the NANO-D group .

Acknowledgement

We wish to thank Brigitte Bidégaray-Fesquet for careful reading of the manuscript and helpful remarks. This work was supported in part by the French National Research Agency (ANR) through the COSINUS program (project SAMSON ANR-08-COSINUS-015).

Appendix

Potential expression

Recall that the Brenner potential is a sum over the chemical bonds.

where b ij is the bond order. The functions V R and V A are respectively repulsive and attractive terms.

These two terms are limited to the nearest neighbors by a simple smooth cut-off function f . The values of Q, A and B n have been determined for each atomic type.

Radical energetics and π bonds influence

The third term in Equation 20 depends on whether the bond has a radical character and is part of a conjugate system.

where F is a tricubic spline. N t i is the total number of neighbors of atom i. N conj ij determines whether the bond is part of a conjugate system from the information of the number of neighbors of each neighbor. This term is decomposed in parts that depend either on atom i, or on atom j :

Dihedral term

The last term in Equation. 20 is a sum over the dihedral angles Θ ijkl

where T is a tricubic spline. This term favors a planar configuration in the case of double C-C bonds. For the description of our incremental update we define :

(1cos(Θ ijkl ))f ik (r ik )f jl (r jl )

For a complete description of the potential, please refer to [START_REF] Brenner | A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons[END_REF].