MODEL FOR SELF-ASSEMBLY OF Br–H AND Br–Br BONDED Br₄Py MOLECULES

1. Introduction

Bottom-up nanotechnology exploiting weak intermolecular interactions has been established as an efficient approach for molecular engineering on different solid surfaces [1]. Among various types of non-covalent interactions, the hydrogen and halogen (X) bonding has been extensively exploited for engineering of novel nanostructures [2, 3].

The X-bonding primarily reveals itself as the electrostatic interaction between the electropositive region (usually the top of a halogen atom) of one molecule and the electronegative region (usually the belt of a halogen atom) of a neighbouring molecule along the  covalent C–X bond  [4]. As a  result, the  X-bonding is highly directional, hydrophobic and tunable by changing the halogen element [5]. It is recognized [6] that a polarized halogen encircled by electronegative charge can engage in both halogen–halogen (X···X) and halogen–hydrogen (X···H) bonding simultaneously. Therefore, in a large variety of chemical and biochemical molecular systems  [7–11], the  X···X interactions are equally present and perhaps even stronger than the X···H ones.

Both C–X···X and C–X···H bondings play an important role in stabilizing the self-assembled supramolecular structures of Br- and I-containing pyrene derivatives on Au(111) [12–14], Ag(111) [14, 15], Cu(111)  [13] and Bi(111)  [16] surfaces. In particular, 1,3,6,8-tetrabromopyrene (Br₄Py) was shown to assemble into two coexisting molecular arrangements on Au(111) at room temperature under ultrahigh vacuum (UHV) conditions [12]. We will name these structures Phase I (Fig. 1(a)) and Phase II (Fig. 1(b)) for consistency with the original paper. Phase I was reported as the dominant one (90% of surface coverage), but both structures consisted of fully intact molecules, and the molecule–substrate interactions were considered to be weak, because initial surface reconstruction remained unhindered [12]. Further annealing at 473 K resulted in cleaving of some C–Br bonds and the emergence of disordered metal-coordinated molecular networks [13]. At 523 K, the monolayer contained some 1D polymer chains connected by C–Au–C bonds. Similar molecular chains with C–Ag–C linkage were observed after the room temperature deposition of the fluorine-rich tetrabromopyrene derivative, via partial debromination catalyzed by the Ag(111) surface [15]. In comparison, the deposition of Br₄Py on Cu(111) leads to partial polymerization into Cu-coordinated tetramers already at room temperature [13].

Phase I was also obtained for the related 1,6-dibromo-3,8-diiodopyrene (Br₂I₂Py) compound, then it was deposited on Au(111) at room temperature [14]. It was found in coexistence with the organometallic phase consisting of dimers with one C–Au–C linkage which were formed due to partial deiodination of molecular pairs. The annealing of this structure at 373 K resulted in the assembly of extended 1D organometallic chains with C–Au–C intra-chain interlinks, driven by full deiodination of Br₂I₂Py molecules, yet without debromination, thus molecular chains were joined together into one monolayer by X···H bonds. On Ag(111) at room temperature, the Br₂I₂Py molecules undergo a partial dehalogenation [14].

Computational modelling of similar supramolecular systems is often employed to aid the experiments, especially for molecules engaging in multiple weak bonds. The nature of mixed X···X and X···H interactions between molecules in experimentally observed molecular arrangements can be revealed by density functional theory (DFT) calculations, which also provide the  geometry associated with the ground state energy [17, 18]. The  collective behaviour of large molecular ensembles is usually studied using less accurate, but more efficient computational methods, such as Monte Carlo (MC) or molecular mechanics modelling. Such studies based on force fields or utilizing purposely built interaction models for treating the assembly of various small organic molecules have been published recently [9, 19–24].

In this work, the  DFT calculations were performed to estimate the  interaction energies for several configurations of Br₄Py molecular pairs, as well as four-molecule clusters representing Phases I and II. The obtained values were used as input in MC simulations of large ensembles of molecules representing the formation process of Phases I and II. Our DFT and MC results suggest that the phase transitions and stability of phases might be explained by the interplay between Br···H and Br···Br contributions to the intermolecular interactions.

2. DFT calculations

To determine the  magnitudes of pair intermolecular interactions, which are the  main driving force for the assembly of Phases I and II, we performed geometry optimizations of two- and four-molecule clusters of Br₄Py by DFT with the ORCA 5.0.3 program package [25]. The clusters represent the typical arrangement of molecules in Phases I and II in the STM images of Ref. [12].

DFT calculations were carried out in the  gas phase using a  ωB97X-D4  [26] range-separated hybrid functional with the  DFT-D4 dispersion correction  [27]. Relativistic effects were accounted for by the  zeroth-order regular approximation (ZORA) [28]. The all-electron relativistically recontracted ZORA-def2-TZVPP basis set with the corresponding SARC/J auxiliary basis was used on all atoms. Integration grid was increased to ‘DefGrid3’, and tight convergence thresholds (‘tightopt’) were applied in all calculations. Some of the larger molecular clusters were pre-optimized with an in-plane constraint, but all final calculations were completed without any geometrical restrictions. The interaction energy of the n-molecule cluster e was calculated as the total energy of the cluster E_n minus energies of its relaxed one-molecule constituents E₁,

First, we performed the calculation of pair interaction energies for two-molecule clusters. The optimized pairs of molecules are shown in Fig.2(a, b). For two interactions of Phase I, which we denote as e_1x and e_1y, we obtained values –4.68 and –5.83 kcal/mol, respectively. The  value of the  pair interaction in Phase II was e₂ = –4.69 kcal/mol (Fig. 2(c)), which is surprisingly very close to e_1x.

In order to check whether we have chosen the  appropriate pair interactions for further calculations, we also performed the  optimization of the  four-molecule clusters representing Phases  I and II. The intermolecular interaction energies of these clusters were –20.75 and –19.20  kcal/mol, respectively, as shown in Fig. 2(d, e). We assume that these interaction energies comprise four pair interactions in both phases in a form 2(e_1x + e_1y) for Phase I and 4e₂ for Phase II. These results quite well correspond to the  results that were obtained for the two-molecule clusters. We did not find the deformation of the Phase II structure, observed with the PBE density functional in the plane-wave DFT calculations [12], that would lead to two slightly different interactions e₂ (along both axes of the square lattice). Our results also demonstrate that the interaction between nearby molecules 1 and 4 (e₁₄) in the Phase I cluster (see Fig. 2(d)) should be much weaker than other pair interactions.

Thus, additionally, we calculated the molecular pair clusters, corresponding to e_1x, e_1y, e₁₄ (Fig. 2(d)), and e₂ (Fig. 2(e)) constrained to their exact optimized positions in the 4-molecule clusters. The results were the  following: e_1x  =  –4.29, e_1y  =  –5.43, e₁₄ = –1.52 and e₂ = –4.67 kcal/mol. Despite geometric constraints, the results agree with the main trend given by the two-molecule clusters. They also show that e₁₄ is small, most likely due to the top–top electrostatic interaction of bromine atoms. It is known [4] that in typical molecules with halogen functional groups the top of halogen is positively charged, while the belt is charged negatively. Therefore, the Br(top)–Br(belt) attraction brings a substantial contribution to the  e_1x, e_1y and e₂ interactions, while the Br(top)···Br(top) contact leads to a rather weak e₁₄ intermolecular interaction.

While trying to determine the fraction of electrostatic Br(top)–Br(belt) attraction in the overall pair interaction, we performed the following DFT calculation. Originally, the e_1x interaction is comprised of two Br(top)–Br(belt) interactions (distance 3.86) and two Br–H interactions (distance 2.87). In the optimized configuration of Fig. 2(a), we artificially substituted two Br atoms by H atoms, thus eliminating the top–belt bonding between Br atoms in the e_1x interaction (Fig. 2(f)). Correspondingly, the  energy of the  remaining interaction was reduced to –2.24  kcal/mol. Neglecting the Br···H bonding, that arises instead of the Br···Br bond, the value of Br(top)–Br(belt) interaction is –4.68 + 2.24 = –2.44 kcal/mol, i.e. more than 50% of the total e_1x energy. The new longer (>4) Br···H intermolecular bonds can also be roughly evaluated as –0.68 kcal/mol by substituting the two other Br atoms by H atoms as demonstrated in Fig. 2(g). Thus, the  share of Br(top)–Br(belt) bonding in the e_1x interaction might increase to more than 60%.

3. MC calculations

In our model, we assume that the molecules move over the sites of a square lattice. A site can be occupied by the  centre of the  molecule only. Since the Br₄Py molecule is symmetric and therefore can act as a binding entity in the same manner at both sides, we assume that it has two states (orientations) A and B on the lattice differing by 90 degrees rotation of the molecule (Fig. 3(a)). The position of a molecule with respect to neighbouring molecules on the lattice is always chosen by taking into account its size and shape to prevent their physical overlap. In Fig. 3(b) the sites around the molecule that cannot be occupied by any part of other molecules are shown by red circles.

In our model, mutual distances between molecules in Phases  I and II are chosen as close as possible to the experimental data [12]. In the previous section, we identified three intermolecular bonding interactions, e_1x, e_1y and e₂, responsible for building the ordered Phases I and II. These interactions occur at specific distances given in (x, y) coordinates of the  lattice. In Phase  I (Fig.  3(c)), molecules of the same state are bound together by the interactions e_1x and e_1y. The latter correspond to intermolecular square lattice distances (5,1) and (1,4), respectively. In Phase II, the ordering of molecules in a checkerboard manner is governed by the  interaction e₂ which arises at the  distance (1,5) between the molecules of different states (see Fig.  3(d)). The  values of interaction energies are taken from our DFT results and scaled with respect to the strongest interaction: e_1y = 1 and e_1x = 0.8, while e₂ was varied within reasonable limits.

MC simulations are performed on a  L  ×  L square lattice of the sizes L = 100 or 200 with periodic boundary conditions. The lattice is initially randomly populated by 280 or 1000 molecules, respectively, with two possible orientations. The larger system was mainly used for testing the  finite size effects. We employ the Metropolis scheme with the  non-local version of the  Kawasaki dynamics algorithm as implemented in previous studies [17, 23]. At each Metropolis step, we first calculate the  initial pair interaction energy E_i of a  randomly selected molecule. After finding an empty site, we attempt to transfer the  molecule with a randomly generated orientational state (A or B) to this new site. Here, the final energy E_f is calculated. The new state and position of a molecule is accepted if ∆E = E_f  – E_i < 0, or accepted with a probability ~ exp(–∆E/k_BT) if ∆E > 0. Simulation is performed by gradually cooling (∆T = 0.05 and 0.02) the system until the minimum energy structure is obtained. We make 10⁵ MC sweeps at each temperature point, where a MC sweep corresponds to L² Metropolis steps.

It should be noted that the ground state interaction energies of Phases I and II are 2(e_1x  +  e_1y) and 4e₂, respectively. The coexistence of the phases is possible when these values are equal. In our gas phase DFT calculations we obtain e_1x + e_1y = –10.51 and 2e₂ = –9.38 kcal/mol using two-molecule clusters (12% deviation from the  coexistence limit), but –9.72 and –9.34 kcal/mol using four-molecule clusters (4% deviation). The similarity of the energies of phases indicates a possibility of coexistence. The fact that the energy of Phase I is always lower corroborates the  experimental finding of Pham et al. [12] which demonstrated that 90% of surface area in the STM image was covered by Phase I and only 10% by Phase II and estimated the difference in the energies of phases as 7%.

The results of our MC simulations show that Phase I (Fig. 4(a)) completely prevails when e₂ is less than 0.92, and Phase 2 (Fig. 4(b)) dominates for e₂ > 0.92, owing to the balance between the sum of interactions in both phases. The threshold value of 0.92 is somewhat larger than the corresponding DFT result (e₂ ≈ 0.8), and it is also slightly larger than anticipated (e₂  =  0.9) from the  balance of the phase energies.

A stable coexistence of both phases in the equilibrated system could be expected at e₂ ≈ 0.92, but we did not observe such a situation at low temperatures when following the typical cooling procedure. It seems that different symmetry of the two structures might be preventing their physical linkage and simultaneous occurrence in simulations. Other effects may also have some influence, such as minor intermolecular interactions, molecule–surface interaction, and possibly rapid quenching process in experiments. However, smaller fragments of Phase I and Phase II were sometimes seen together in snapshots (Fig. 4(c)) at the temperatures close to the phase transition (molecular ordering) point, T_c ≈ 0.38. In some cases, we also noticed a repetitive forward and backward transition from Phase I to Phase II at successive temperature points near T_c when the temperature step was very small.

The MC procedure tends to determine the structure with the lowest energy, but it cannot take into account the  freezing of other structures due to quenching of the  system or the  substrate effect. Such effects might lead to the coexistence of phases seen in experiments [12, 13] that were performed at room temperature, certainly lower than the ordering temperature. Thus, we decided to perform the MC simulation at one temperature point well below T_c, where the movement of molecules is hindered. In our calculation at T = 0.08 after 10⁵ MC steps we found large domains of both phases (Fig. 4(d)), resembling the molecular arrangement observed in experimental STM images.

4. Conclusions

In this paper, we proposed a  model and studied the formation of self-assembled structures of Br₄Py molecules. By performing the DFT calculations, we evaluated the pairwise intermolecular interactions and chose those which played the most important role in molecular ordering. For MC calculations, we developed a two-state model on a square lattice with three interactions chosen from the  DFT results. In our MC simulations, we observed the pure Phases I and II, as well as their transient coexistence at higher temperature.

References

Br–H IR Br–Br RYŠIAIS SUSIJUNGIANČIŲ Br₄Py MOLEKULIŲ SAVAIMINIO SUSITVARKYMO MODELIS

^a Fizinių ir technologijos mokslų centras, Vilnius, Lietuva

^b Vilniaus universiteto Fizikos fakultetas, Vilnius, Lietuva

Santrauka

1,3,6,8-tetrabromopireno (Br₄Py) molekulių susitvarkymui į dvi skirtingas plokščias struktūras (I ir II fazes) aprašyti yra siūlomas gardelės modelis, grindžiamas porinėmis tarpmolekulinėmis halogeno–vandenilio (Br–H) ir halogeno–halogeno (Br–Br) sąveikomis ant kvadratinės gardelės. Porinių sąveikų energijos įvertinamos skaičiuojant tankio funkcionalą dviejų molekulių ir keturių molekulių klasteriams. Rezultatai rodo, kad apie 60 % sąveikų energijos sudaro elektrostatiniai Br (viršus) – Br (šonas) ryšiai, o likusi dalis atsiranda dėl Br–H ryšių. Gautos sąveikų energijų vertės toliau naudojamos Monte Karlo skaičiavimams atlikti, remiantis pasiūlytu gardelės modeliu. I fazei  susiformuoti reikalingos dvi pagrindinės porinės sąveikos, o išsidėstyti į II fazę užtenka vieno sąveikos tipo. Gauti rezultatai paaiškina abiejų fazių susidarymo ypatumus.