Marc-André Dubois

Xavier Bouju

Alain Rochefort

Toward interactive scanning tunneling microscopy simulations of large-scale molecular systems in real time

Keywords: numbers: 73.20.-r, 71.15.Nc, 82.65.+r STM, surface adsorption, large-scale simulation, molecular systems, parallel computing

published or not. The documents may come

Toward interactive scanning tunneling microscopy simulations of large-scale molecular systems in real time

I. INTRODUCTION

Scanning tunneling microscopy (STM) is an imaging technique that exploits the quantum nature of electrons mainly to study the electronic states of a substrate or adsorbed molecular systems, from a single molecule to a molecular layer on a surface. Given the atomic resolution and the nature of the interactions between the sample and the STM tip, this technique can also be used to manipulate atoms and molecules in a controlled manner [1][2][3][4] to build complex nanostructures in a bottom-up fashion, for instance.

Despite STM has a very good resolution, a direct interpretation of experimental images can be difficult considering the complexity of the electron propagation into the media by tunneling. For example, the STM signature of an adatom on a metal surface can be quite different depending on its adsorption site. [5][6][7] A similar behavior is observed with molecules where the various adsorbed sites i.e. the different molecule-surface interactions) play a role as important on the resulting STM contrasts as the variations of internal geometry of the molecules. 8,9 In addition, the non-trivial relationship between tunneling probability and molecule-surface interactions cannot be easily a) Electronic mail: alain.rochefort@polymtl.ca generalized. Thus, the development of simulation tools that support fundamental aspects related to surface science is highly relevant [9][10][11][12][13][14][15] to support and settle definitive conclusions on the structural and electronic properties of the probed species. Two steps are necessary to perform such simulations: first, the most stable geometry of the adsorbate has to be determined, and second, the STM image needs to be evaluated. Although each step consumes a significant amount of computational time that depends of the level of accuracy needed and the size of the system investigated, the evaluation of the electronic structure always remains the bottleneck of STM simulations.

The quality of simulated STM images mostly depends on the tip-molecule-surface junction model from which the electronic structure is evaluated. Most relevant models follow rigorous theoretical frameworks in which the search of a good accuracy in the determination of the electronic structure and the electron transport properties strongly impacts the computational time consumption. For example, high level theory techniques such as Density Functional Theory (DFT) methods usually provide excellent agreement with experimental results for molecular systems containing a few hundred of atoms. An interesting theoretical approach has been recently developed by Jelínek and co-workers with a CO (or Xe atom) functionalized tips to provide high-resolution images. [16][17][18][19] During the tip scanning, the CO tip apex is relaxed at each lateral position of the probe. The flexibility of such a COterminated tip has to be taken into account to obtain calculated images are in good agreement with experiments. The procedure of tip relaxation is only considered for low tip-surface distances, and is not relevant for large tip heights or with a clean STM tip. Another case where relaxation has to be considered is single-molecule STM manipulations. Whatever the manipulation mode (pushing or pulling), a mechanical modification of the adsorbed system is revealed by calculated STM signals that can be obtained with a numerical STM combining Elastic Scattering Quantum Chemistry (ESQC) technique and molecular dynamics. Such a numerical tool is able to mimic as close as possible an actual STM setup comprising feedback loops and molecular dynamics coupled to the tunnel current calculations. At each position of the tip during the constant-current imaging process, the adsorbate is free to relax within an adiabatic approximation. Such a procedure delivers variation of the junction resistance as a function of the tip-surface distance 20 or tunnel current signatures during a controlled manipulation. 21,22 Unfortunately, such approaches do not allow large-scale calculations and on-the-fly STM simulations.

Here, we propose a strategy to reach real-time imaging in which the image production time is extremely fast. Staying in the imaging regime, we considered that the mechanical action of the tip remains minimal. We also consider the possibility of a local structural or chemical modification of surface species called intrusion in the following. Such an intrusion on an adsorbed system could be induced by an external user, 23 by a STM tip through a mechanical action, 1 by a local electric field generated at the tip apex 24,25 or by using inelastic-electron tunneling induced manipulation processes. 26 A general approach in STM simulation that includes intrusion within the limits of Tersoff-Hamann and Bardeen theories was already developed by our group over the last years. 23 However, despite considering an explicit intrusion followed by a very rapid calculation of the STM image, this approach does not consider any molecular relaxation after the intrusion, and this essentially leads to a static STM imaging mode. The objective of the present computational development is to produce a numerical toolm that can rapidly and accurately provide a reasonable adsorbate geometry and more importantly, a STM image that can be compared directly to experimental results.

In the present work, we have developed an STM simulation framework called DyFlex that considers the concept of intrusion, that accounts for subsequent molecular relaxation in the presence of a metallic surface, and produces nearly real-time STM images. We propose a strategy that combines simple but adequate models within a computational implementation that also benefits from modern architectures involving graphical processing units (GPUs). Section II presents a summary of the formalism we used to evaluate the energy, to consider molecular relaxation and then to obtain STM im-ages. We first validate our integrated tool for molecular systems by checking the optimization for two different datasets: a subset of the Cambridge Structural Database (CSD) 27 containing 500 molecules with geometries highly dictated by strong intramolecular interactions, and the s22 dataset 28 that contains 22 complexes bonded through non-covalent interactions. We verify the robustness of our optimization scheme by considering also the complexes of the s22x5 dataset, 29 which is an extension of s22 for non-equilibrium geometries. Next, we show an example of intrusion in which we have studied the first step of the Ullmann reaction that involves the breaking of a halide-carbon bond, followed by the relaxation of the fragments on a metal (Cu) surface. The resulting optimized geometries for the chemical intermediates species are in very good agreement with the experimental and DFT results. Finally, we show how the use of GPUs drastically improves the computation of STM images associated to a known supramolecular aggregations of 1,4-di(4',4"-pyridyl) benzene (NN) on a highly doped Si(111)-B surface. 30

II. METHODS

This section introduces the different theoretical components we have combined to consider molecular relaxation. First, we present how we use the Atom Superposition and Electron Delocalization Molecular Orbital theory (ASED-MO) combined to an empirical van der Waals functional form to compute the electronic structure and total energy of the molecular systems. We then describe how we use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) non-linear optimization algorithm along the molecular relaxation process. Finally, we present the formalism to evaluate the tunneling current, and how it was implemented over a hybrid CPU (central processing unit)/GPU infrastructure.

A. Electronic structure calculations

The ASED-MO theory

The well-known extended Hückel Molecular Orbital theory (EHMO), developed by Hoffmann, 31 is a simple semi-empirical quantum theory where the Hamiltonian elements (H ij ) are proportional to the overlap matrix elements (S ij ) such:

H ij = K H ii + H jj 2 S ij , (1)
where K is the Wolfsberg-Helmholtz constant. The values of H ii(jj) are defined by the so-called valence state ionization potential (VSIP) 31 of the atomic orbital centered on atom i (j). As in non-orthogonal tight-binding approaches, overlap matrix elements are explicitly evaluated with EHMO. The total energy E EHMO is simply obtained by summing the energies over the occupied orbitals. This method is well-known to give relevant qualitative eigenvalues and eigenfunctions. [31][START_REF] Albright | Orbital Interactions in Chemistry[END_REF][START_REF] Ballhausen | Molecular orbital theory: an introductory lecture note and reprint volume[END_REF] However, since it does not explicitly include two-body electrostatic interactions, EHMO generally fails to give realistic total energy (geometry) when used with non-linear optimization techniques. [START_REF] Calzaferri | [END_REF] To overcome this deficiency, Anderson proposed the Atom Superposition and Electron Delocal-ization Molecular Orbital theory (ASED-MO). 35,36 Following the ASED-MO scheme, the electronic density is split into two contributions: a perfectly following term where the atomic electron density is centered on atom (ρ PF ) and follows perfectly the nuclei. The nonperfectly following term (ρ NPF ) is associated to the formation of chemical bonds where the density is shared across the molecule (then not follows perfectly the nuclei). The energy associated to ρ PF is the sum of repulsive energy between atomic cores and the attractive energy between the nuclei and the electronic density:

E PF (R) = i i-1 j Z i Z j R ij - 1 2 i i-1 j Z i ρ j |R ij -r| dr + Z j ρ i |R ij -r| dr , (2)
where Z i(j) and ρ i(j) are the atomic charge and the electronic density of the atom i(j), and R ij is the distance between atoms i and j. The energy associated to ρ NPF is attractive and can be simply represented by the energy difference between the whole molecule and the isolated atomic components:

E NPF (R) ≈ E EHMO (R) - i I i , (3)
where I i is the ionization potential of the atom i and E EHMO is the molecular energy obtained from the EHMO theory. In the present work, the Wolfsberg-Helmholtz constant (K) in off-diagonal Hamiltonian terms is replaced by a distance dependent term: 34,37

H ij = {1 + κ ij exp [-δ ij (R ij -R cov )]} H ii + H jj 2 S ij , (4)
where κ ij and δ ij are empirical parameters while R cov is the covalent radius between the orbitals i and j. Following this theoretical framework, the molecular orbitals (ψ) are represented by a linear combination of atomic orbitals (φ) with the form of two Slater Type Orbitals (STO) centered over an atom:

38 φ i (r) = βY l m (r)r n-1 c 1i e -ζ1r + c 2i e -ζ2r , (5)
where β is a normalization constant, Y l m is the normalized spherical harmonics and ζ µ and c µi are parameters related to the primitives representing the atomic orbital. The ASED-MO method has already been advantageously used to simulate rather large systems with a divideand-conquer technique as implemented in the SAMSON code. 39,40 Additionally, it was used to obtain optimized geometries of adsorbed molecules on metallic [41][42][43][44] and semiconducting surfaces. 45,46

Contribution from van der Waals interactions

Dispersive interactions play an important role in numerous physical process, 47,48 they are of primary importance for nearly all organized molecular systems that are physisorbed [START_REF] Oura | Surface Science[END_REF] on surfaces. To consider these long-range and weak interactions, we have used an empirical potential term similar to the -1/r 6 term of a Lennard-Jones potential. Such energy correction was recently introduced by Grimme 47 and already often used in conjunction to Density Functional Theory (DFT) calculations to correctly describe weakly bonded molecular systems such as van der Waals (vdW) complexes: 50

E disp (R) = -s 6 N -1 i=1 N j=i+1 C ij 6 R 6 ij f dmp (R ij ) , (6)
where C ij 6 et R ij are the coefficients obtained by Grimme, and s 6 is a scaling factor. In order to avoid an important attractive contribution at small R ij value, a damping function f dmp is used:

f dmp (R ij ) = 1 1 + exp -d Rij R ij 0 -1 , (7)
where d and R ij 0 are, respectively, an empirical constant and the sum of van der Waals radius for the atoms i and j. Finally, by considering the various energy contributions described above, E NPF , E PF and E disp , we obtain the total energy of the system:

E TOT = E NPF + E PF + E disp . (8)

B. Geometry relaxation

Since the molecular relaxation that follows an intrusion is generally localized, it is then reasonable to use a local optimization algorithm to explore the potential energy surface for finding the minimum of the energy function defined in (8). The more general approach in optimizing a non-linear function consists to follow a descending direction. At least, one needs to know the function f (x) to optimize and a vector of coordinates x 0 that acts as the starting point of the algorithm:

x k+1 = x k + α k d k , (9)
where α k and d k are respectively, the step length and the descent direction. We choose to consider an adaptive step-length based on the Wolfe conditions 51,52 for the determination of the optimal step-length while using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method 53 to evaluate the descent direction. This method is frequently used in non-linear optimization problems [START_REF] Nocedal | Numerical Optimization[END_REF] and is known to perform well even with non-smooth functions. [START_REF] Lewis | [END_REF] The BFGS method was especially chosen for its efficiency to converge toward a local minimum. This should work fine for systems where the starting point is not drastically different from the relaxed geometry.

Force computation

As described previously, within a local minimization scheme, the forces acting on each degree of freedom need to be determined. In the ASED-MO theory, these can be obtained analytically 56 over the primitives (µ, ν) of orbital i with respect to the coordinate x of a given atom:

∇ x E = µ ν P µν ∇ x H µν - µ ν W µν ∇ x S µν +∇ x E PF + ∇ x E disp , (10)
where P µν and W µν are respectively, the following density and the energy-weighted density matrices:

P µν = i≤N/2 2V µi V νi , W µν = i≤N/2 2 i V µi V νi , ( 11
)
and i is the energy of the orbital i. This analytic form has the advantage of reducing the quartic complexity of the numerical gradient to a cubic power. 57

C. Evaluating the tunneling current

The tunneling current I between a tip and a surface under a bias voltage has an exponential relation to the tip-surface distance, and depends strongly on the applied voltage and the local density of states (LDOS) of the sample: 58

I ∝ E F En=E F -eV T |ψ n | 2 e -2βz , (12)
where ψ n is the state n of the substrate, β is a parameter related to the tunnel barrier and z is the distance between the tip and the substrate. From Eq. ( 12), we can see that the current is proportional to the propagation states (ψ n ) contained within an energy window defined by the Fermi level (E F ) and the applied bias voltage (V T ).

In order to compute an entire STM image, one needs to describe the electronic structure of the tip and the substrate, and to consider the interactions between them. [START_REF] Paz | [END_REF] Most rigorous methods explicitly considered these interactions in the theoretical tunneling framework. [60][61][62] For short tip-surface distances, inelastic electrons scattering becomes more significant, and the tunneling process needs to be more carefully described. [START_REF] Paz | [END_REF] In normal operation conditions of STM, the tip-substrate distance is relatively larger (∼5-10 Å), and it is reasonable to adopt a perturbative approach to obtain the main STM image contrasts. 63,64 This is the central idea in the Tersoff-Hamann approximation (TH), 64 where the tip is approximated by a spherical potential. The interpretation of the tunneling current then becomes exclusively dependent of the calculated substrate states ψ s lying in the energy window (δ) such as:

I(r) ∝ V i |ψ s i (r)| 2 δ(E i -E F ) . ( 13
)
Since a tunneling current is calculated for each pixel of the STM image, the algorithmic complexity of such approach is O(nmij) for a (n × m) image with an energy window containing i molecular states, constructed from a linear combination of j atomic orbitals. Fortunately, each pixel can be calculated independently. Tersoff-Hamann technique remains an approximation of the Bardeen theory where the tip is more explicitly considered. [START_REF] Chen | Introduction to Scanning Tunneling Microscopy[END_REF]63,[START_REF] Tersoff | Scanning Tunneling Microscopy[END_REF] Although we can also produce STM images at the Bardeen theory level with our computational framework, 23 we are focusing on the TH approximation here because it generally reveals the most important features of STM images, especially for two-dimensional molecular systems, and it is coherent with our real-time STM simulation objective. In this special context, a strategy implying the parallelization of the pixels computation becomes even more highly suitable.

D. Parallel computing of STM images on GPUs

An STM microscope can operate in two imaging modes:

1. The constant height mode consists in scanning the surface with a fixed tip-surface distance and recording the fluctuation of the current directly. This mode has the advantage to be fast, but is limited to the study of non-corrugated surfaces to avoid the risk of collision between the tip and the surface.

2. The constant current mode is the most used mode. It consists in scanning the surface and adjusting the tip-surface distance z at each pixel with the help of a feedback loop to re-insure a constant current. This mode is used to extract the topography of the surface.

Parallel computation of the constant-height mode is straightforward and can be done for each pixel in a single step with Eq. ( 13). It is more complex for the constantcurrent mode where the current needs to be evaluated at multiple heights (z) to localize the targeted current. In these steps, we define a domain along the z axis and we gradually converge to a given height through a dichotomic search algorithm at each pixel. While the current I is computed using the GPU, the dichotomic search is performed on the CPU, inducing an additional computational cost related to the communication between the CPU and the GPU memories. In hybrid CPU/GPU implementations, the memory exchanges between the host and the device(s) must be minimized to lower the communication latencies that limits the gain performance of using GPUs. In the present work, the memory initialization takes the form shown in Table I. For example, a metallic surface containing 1 000 atoms involves [(1 000s) + (3 000p) + (5 000d)] = 9 000 atomic orbitals. If we consider an image of 100 Å by 100 Å with a pixel resolution of 0.3 Å, 111 111 pixels need to be calculated. By combining these information, the total memory required is: Mem = (16 × 111 111 px) + ((8 + 40) × 9 000 orb.) + (8 × 9 000 2 orb.) = 650 MB .

Despite the data transfer needed for the initialization, the memory associated with the storage of the pixels represents a small proportion (1.8 Mo) of the total memory needed (650 Mo). The additional cost of using a dichotomic search on the CPU in the constant-current mode is then negligible compared to the performance gain of using GPUs to compute the image. The computational flow chart of DyFlex is shown in Fig. 1. It summarizes the different decision pathways among each modules of the entire program, including the interactions between modules and the information shared.

III. RESULTS AND DISCUSSION

We start to validate our computational approach that includes ASED-MO and vdW corrections by determining the geometries of molecular systems in gas phase with high covalent bonding character but also of complexes where vdW interactions play a major role. The first group is a subset of the CSD database of about 500 crystallographic structures, 27 and the second is determined by considering the same s22 dataset while starting from a well-defined off-equilibrium geometries. We complete this validation scheme by studying the on-surface 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18

DRMSD ( Å)

Probability (A.U.)

FIG. 2: Distribution of the root-mean-square deviation (DRMSD) in bond lengths obtained with ASED-MO+vdW with respect to the molecules contained in the CSD database subset.

Ullmann reaction in which we are comparing the calculated geometries of physisorbed and chemisorbed species on a Cu(111) surface with DFT results, and simulated STM images with experiments. The general quality of the calculated distances is made by evaluating the rootmean-square deviation (DRMSD) over all pairs of atoms between the reference geometry (experimental data or DFT calculations) (X ref ) and the calculated geometry with ASED-MO+vdW (X) of the targeted molecule:

DRMSD X ref , X = 1 N (N -1) × i =j δ x ref i , x ref j -d (x i , x j ) 2 ,
where N corresponds to the number of atoms in the molecular systems. Finally, we look at the performance improvement obtained from GPUs in the computation of STM image of a supramolecular arrangement of (1,4-di(4',4"-pyridyl)-benzene, NN) molecule on a highly doped Si(111)-B surface. 30 A. Highly covalent bonded molecules Figure 2 reveals that the average deviation of intramolecular bond length is about 0.022 Å. In fact, 95 % of the molecules contained in the CSD subset, have a deviation that is less than 0.048 Å. Such deviations that our optimized geometries are mainly indistinguishable from the reference geometries, and are more than satisfying. Our results are also comparable to typical cases of geometry optimizations with much more sophisticated methods such as Hartree-Fock (HF) or DFT, 66 but at a much lower computational cost.

B. Weakly bonded complexes

In order to check the ability of our ASED-MO+vdW approach to deal with van der Waals interactions, we have compared the optimized geometries from the s22 dataset with and without the E disp term. Figure 3 shows a major improvement in the determination of bond length when adding the dispersion term (E disp ) into the original ASED-MO formalism, where the typical DRMSD decreases from 0.13 Å to 0.027 Å. The complexes of the s22 dataset can be classified in three categories: 1) Hbonding complexes (1-7), 2) dispersive or van der Waals complexes (8-15) and, 3) a mixture of 1) and 2) complexes (16-22). From Fig. 3, the improvement observed for complexes characterized by van der Waals interactions by including E disp is spectacular, and more especially when the complex involves only dispersion interactions such in CH 4 -CH 4 (8), C 2 H 4 -C 2 H 4 (9) and CH 4 -C 2 H 6 (10) complexes. As an example, Fig. 4 shows a superposition of the optimized structures of a methane-benzene complex obtained with (green) and without (blue) the E disp term with respect to the reference structure (gray: carbon atoms, white: hydrogen atoms). The different colors helps to evaluate the influence of both energetic terms. The geometry obtained without dispersion correction in blue is quite far from reference geometry, while giving appropriate internal molecular geometry. The geometry obtained with correction in green shows very good methane-benzene distance, but tends to slightly contract the CH bonds in methane and benzene. The ASED-MO scheme contains ρ i(j) terms representing partial charges that allows a partial stabilization of weakly bonded complexes (8-10) through electrostatic interactions.

C. Robustness of relaxation

The robustness of our optimization scheme was tested by optimizing the s22x5 dataset, a variant of the s22 FIG. 4: Effect of dispersion on the more stable geometry of a weakly bonded methane-benzene complex (complex 10). The complex without correction is in blue, with the correction is in green, and the reference complex is grey (C) and white (H). DRMSD is 0.13 Å and 0.027 Å without and with the correction, respectively. dataset in which the intermolecular distances of the reference complexes are stretched or squeezed by a geometrical factor S that varies from 0.9 to 2.0. This approach allowed us to more quantitatively characterize the influence of the starting geometry on the overall optimization process to obtain reasonable geometries from which a STM image could be calculated. Figure 5 illustrates that the DRMSD obtained after optimization remains reasonable for a factor S < 1.5. At higher S values, the final optimized geometries diverge from the expected final geometries. As for the s22 dataset, the high dependence of S on the final geometries can be explained by the nature of the intermolecular interactions. The observed results suggest that our approach to obtain reliable geometries can be successfully used when the atomic relaxation is limited to an average radius of 0.6 Å around the real optimized geometry.

D. Molecular relaxation on a surface: the Ullmann reaction

The surface-assisted Ullmann coupling of aryl halides represents a classical STM example where the tip is used to break chemical bonds, to mechanically manipulate molecular species and to image the associated chemical reaction. Hla et al. first realized the synthesis of a biphenyl molecule with a STM from the Ullmann coupling of two adsorbed iodobenzene molecules on a Cu(111) surface. 67 The complete reaction on the surface consists in three steps: (1) the breaking of the C-I bond by a voltage pulse, (2) the migration of the intermediate chemical species by pushing mode with the STM tip and, FIG. 5: Comparison of DRMSD obtained with non-equilibrium geometries of s22 dataset where the distance between molecules is multiplied by a S factor of 0.9 (blue), 1.0 (black), 1.2 (green), 1.5 (orange) and 2.0 (red).

(3) the formation of a C-C bond between adsorbed phenyl radicals with another voltage pulse to form a biphenyl molecule. This reaction was recently revisited with the help of DFT calculations by Björk et al. 68 where they obtained geometries and adsorption conformations of the chemical intermediates.

We have used DyFlex in the intrusive mode to investigate the first step of the Ullmann reaction in which we put a halide-benzene molecule on the surface, we artificially break the C-halide bond and then, we let the phenyl fragment to relax on the surface and determine the most stable chemical intermediate. STM images can also be rapidly calculated along the molecular relaxation process. We start first by placing a benzene molecule on the Cu(111) surface model containing 118 Cu atoms shown in Fig. 6. The size of the Cu cluster model used here to simulate the Cu(111) surface is sufficiently large to minimize the effects related to open boundary conditions of finite systems and to reproduce the main electronic structure features of the extended surface. We then substitute a hydrogen by a bromine atom (chemical intrusion), and we allow the bromobenzene molecule to relax on the fixed geometry Cu surface. Figure 6 shows the optimized geometry of bromobenzene (C 6 H 5 Br) on Cu(111) obtained at the ASED-MO+vdW level. The molecule is slightly tilted by 6 • and the mean adsorption height of the central benzene ring is about 2.94 Å, in very good agreement with the value (2.93 Å) obtained by Björk et al. 68 with DFT method. In addition, Fig. 6(b) shows that the bromobenzene molecule is physisorbed on a three-fold hollow site on Cu(111), also found by DFT, 68 and confirmed the ability of the present approach in reproducing relevant molecular geometries of weakly adsorbed species.

To provoke the dehalogenation of bromobenzene that can be experimentally induced by a voltage pulse from the STM tip or by a thermal activation on the surface, we 7 where the phenyl moiety becomes chemisorbed on the surface. As a result, the molecule is bonded to a single surface atom at a C-Cu distance of 2.05 Å and the phenyl fragment is strongly tilted on the surface. Similar adsorption configuration was observed with DFT, where a C-Cu distance of 1.99 Å was reported. 68 Such small differences in the geometry of an adsorbed species could not be distinguished within the STM technique. Since our goal is to rapidly produce relevant STM images, the structure calculated with ASED-MO+vdW is quite satisfying. Concerning the simulation of STM images, we have limited our comparison with experiments to the physisorption of bromobenzene and the chemisorption of a phenyl fragment on the Cu(111) surface as shown in Fig. 8. As observed in the experimental STM image [Fig. 8(b)], the bromobenzene molecule appears as an asymmetric protrusion where a maximal contrast is centered over the halogen atom [Fig. 8(a)] and witha mirror symmetry along the Br-H line that is slightly shadowed by the rendering. In addition, the molecular fragment resulting of the dehalogenation of bromobenzene can be easily recognizable through its pseudo-

E. Performances

Finally, we have benchmarked DyFlex in generating STM images using computer nodes built from CPU (2.5 GHZ 10 cores Intel Xeon Ivy Bridge E5-2670 V2, 128 GB RAM) and GPU (NVidia K20 GPU, 2496 cores at 706 MHz, 5 GB GDDR5 memory) components. The benchmark consists of a comparison of computing time on CPU only and then on GPU to generate various image sizes of a recently studied supramolecular arrangement of 1,4-di(4',4"-pyridyl) benzene (NN) molecules that form nanolines on a highly-doped Si(111)-B surface. 30 To give a fair comparison between the two architectures (CPU vs GPU), we compare the GPU version to a parallelized and vectorized CPU version.

Figure 9 reveals a drastic difference in computing time, where the GPU version is faster than the CPU version by about 30 for large images. For comparison, the slopes for the CPU and GPU are respectively, 0.30 and 0.009 min/px. The use of GPU over CPU clearly becomes beneficiary when large molecular models are used. Moreover, the strong linearity in computing time over the and the GPU version (120 nm 2 ). For both images, the computing time is the same. These results give the user the ability to observe not only the STM contrasts obtained from an isolated molecule, but also give access to the supramolecular patterns that appear in extended systems. As it was recently reported, 30 the ability of our solver to consider relaxation while computing efficiently STM images of supramolecular system helps to identify the key aspects of the conformational arrangement of the molecules such as the formation of molecular dimers.

IV. CONCLUSION

We have presented a computational approach related to chemical and structural modifications of molecular systems with subsequent relaxation. Our approach gives reasonably good geometries for a rather large set of molecules, and where the consideration of a dispersion energy term considerably improves the calculated geometries of a large set of molecules and van der Waals complexes. In addition, we have demonstrated that our local optimization algorithm can be used with confidence even for geometries that are not necessarily close to the minimum of the potential energy surface, especially for adsorbed molecules. Finally, we developed a procedure to speed-up the calculations of STM image for comparison with experimental results. We showed that the use of GPUs improves the performance of our solver by a factor of 30 to generate simulated STM images. The present development in the simulation of nearly real-time STM images exploiting new computing architectures opens ac-

FIG. 3 :

 3 FIG.3: Distance root-mean-square deviations (DRMSD) obtained for the s22 dataset without (red) and with (blue) the dispersion energy corrections to the ASED-MO energy.

FIG. 6 :FIG. 7 :

 67 FIG. 6: (a) Side view and (b) top view of the optimized geometry of a bromobenzene molecule adsorbed on a Cu(111) surface. The tilt angle of the molecule is 6 deg and the adsorption height of the central benzene ring is 2.94 Å.

FIG. 8 :

 8 FIG. 8: Comparison of the simulated STM images (a, c) with DyFlex and measured (b, d) experimentally of benzene-halide (C 6 H 5 X) (upper panels) and phenyl fragment (C 6 H 5 ) (lower panels) on Cu(111). (Both (b) and (d) images were adapted with permission from Hla et al. 67 . Copyright 2000 American Physical Society.)

FIG. 9 :

 9 FIG. 9: Computing time to generate STM images as a function of domain size with CPU (blue) and GPU (red) architectures.

FIG. 10 :

 10 FIG. 10: STM simulations of NN nanolines on a Si(111)-B surface obtained with CPU and GPU architectures. The experimental image is also provided for reference. (Adapted with permission from Zhan et al. 30 . Copyright 2017 American Chemical Society.)

TABLE I :

 I Summary of the memory elements transfered between the CPU to the GPU when simulating an STM image in DyFlex. n is the number of pixels, m the nuber of orbitals of the system. The size of an arbitrary element is defined in bytes.

			Molecular model
			Formation of the matrices
			Diagonalization
			S, H, V ,
	yes	Convergence?	Energy evaluation Forces computation	New geometry
		no
			Evaluation of α
			Acceptable?	no
			yes
		Geometry
			Energy window evaluation
		V ,
			Evaluation of Fermi factors
		Constant current		Constant height
	Computation of the current		Computation of the current
	for each pixel		for each pixel
			no
	Current comparaison	Acceptable?
			yes
				Image
	FIG. 1: Flow chart of the different steps realized by the
	DyFlex solver to consider molecular relaxation following
	an intrusion in the computation of a STM image. The
	steps identified in red are entirely computed on the
			GPU.

V. ACKNOWLEDGMENTS

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). This work would not be possible without the computational resources provided by Calcul Québec and Compute Canada. We acknowledge Programme Investissements d'Avenir Grant No. ANR-11-IDEX-0002-02, reference ANR-10-LABX-0037-NEXT for support.