Lithuanian Journal of Physics, Vol. 62, No. 4, pp. 235–242 (2022)
© Lietuvos mokslų akademija, 2022
Received 12 September 2022; accepted 3 October 2022
A model, based on the pairwise intermolecular halogen–hydrogen (Br–H) and halogen–halogen (Br–Br) bonding, is proposed to describe the self-assembly of Br4Py molecules into two different planar structures (Phase I and Phase II). The pair bonding interactions are calculated by the density functional theory for the two-molecule and four-molecule clusters. It is shown that about 60% of bonding strength is due to the electrostatic Br(top)–Br(belt) interactions, while the remaining originates from Br–H interactions. The obtained values of pair interactions are further used for Monte Carlo calculations. The model for these calculations is proposed on a square lattice. The two main pair interactions are needed for the emergence of the Phase I ordering, while the Phase II ordering is obtained using a single interaction. The obtained results explain the emergence of both phases.
Keywords: molecular self-assembly, statistical models of phase transitions, Monte Carlo simulations, density functional theory calculations
Bottom-up nanotechnology exploiting weak intermolecular interactions has been established as an efficient approach for molecular engineering on different solid surfaces . 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 . As a result, the X-bonding is highly directional, hydrophobic and tunable by changing the halogen element . It is recognized  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)  and Bi(111)  surfaces. In particular, 1,3,6,8-tetrabromopyrene (Br4Py) was shown to assemble into two coexisting molecular arrangements on Au(111) at room temperature under ultrahigh vacuum (UHV) conditions . 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 . Further annealing at 473 K resulted in cleaving of some C–Br bonds and the emergence of disordered metal-coordinated molecular networks . 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 . In comparison, the deposition of Br4Py on Cu(111) leads to partial polymerization into Cu-coordinated tetramers already at room temperature .
Phase I was also obtained for the related 1,6-dibromo-3,8-diiodopyrene (Br2I2Py) compound, then it was deposited on Au(111) at room temperature . 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 Br2I2Py molecules, yet without debromination, thus molecular chains were joined together into one monolayer by X···H bonds. On Ag(111) at room temperature, the Br2I2Py molecules undergo a partial dehalogenation .
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 Br4Py 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.
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 Br4Py by DFT with the ORCA 5.0.3 program package . The clusters represent the typical arrangement of molecules in Phases I and II in the STM images of Ref. .
DFT calculations were carried out in the gas phase using a ωB97X-D4  range-separated hybrid functional with the DFT-D4 dispersion correction . Relativistic effects were accounted for by the zeroth-order regular approximation (ZORA) . 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 En minus energies of its relaxed one-molecule constituents E1,
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 e1x and e1y, we obtained values –4.68 and –5.83 kcal/mol, respectively. The value of the pair interaction in Phase II was e2 = –4.69 kcal/mol (Fig. 2(c)), which is surprisingly very close to e1x.
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(e1x + e1y) for Phase I and 4e2 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 , that would lead to two slightly different interactions e2 (along both axes of the square lattice). Our results also demonstrate that the interaction between nearby molecules 1 and 4 (e14) 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 e1x, e1y, e14 (Fig. 2(d)), and e2 (Fig. 2(e)) constrained to their exact optimized positions in the 4-molecule clusters. The results were the following: e1x = –4.29, e1y = –5.43, e14 = –1.52 and e2 = –4.67 kcal/mol. Despite geometric constraints, the results agree with the main trend given by the two-molecule clusters. They also show that e14 is small, most likely due to the top–top electrostatic interaction of bromine atoms. It is known  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 e1x, e1y and e2 interactions, while the Br(top)···Br(top) contact leads to a rather weak e14 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 e1x 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 e1x 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 e1x 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 e1x interaction might increase to more than 60%.
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 Br4Py 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 . In the previous section, we identified three intermolecular bonding interactions, e1x, e1y and e2, 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 e1x and e1y. 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 e2 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: e1y = 1 and e1x = 0.8, while e2 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 Ei 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 Ef is calculated. The new state and position of a molecule is accepted if ∆E = Ef – Ei < 0, or accepted with a probability ~ exp(–∆E/kBT) 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 105 MC sweeps at each temperature point, where a MC sweep corresponds to L2 Metropolis steps.
It should be noted that the ground state interaction energies of Phases I and II are 2(e1x + e1y) and 4e2, respectively. The coexistence of the phases is possible when these values are equal. In our gas phase DFT calculations we obtain e1x + e1y = –10.51 and 2e2 = –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.  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 e2 is less than 0.92, and Phase 2 (Fig. 4(b)) dominates for e2 > 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 (e2 ≈ 0.8), and it is also slightly larger than anticipated (e2 = 0.9) from the balance of the phase energies.
A stable coexistence of both phases in the equilibrated system could be expected at e2 ≈ 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, Tc ≈ 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 Tc 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 Tc, where the movement of molecules is hindered. In our calculation at T = 0.08 after 105 MC steps we found large domains of both phases (Fig. 4(d)), resembling the molecular arrangement observed in experimental STM images.
In this paper, we proposed a model and studied the formation of self-assembled structures of Br4Py 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.
 D.P. Goronzy, M. Ebrahimi, F. Rosei, Arramel, Y. Fang, S. De Feyter, S.L. Tait, C. Wang, P.H. Beton, A.T.S. Wee, et al., Supramolecular assemblies on surfaces: nanopatterning, functionality, and reactivity, ACS Nano 12, 7445–7481 (2018).
 G.R. Desiraju and T. Steiner, The Weak Hydrogen Bond in Structural Chemistry and Biology (Oxford, New York, 1999).
 A. Priimagi, G. Cavallo, P. Metrangolo, and G. Resnati, The halogen bond in the design of functional supramolecular materials: recent advances, Acc. Chem. Res. 46, 2686–2695 (2013).
 T. Clark, M. Hennemann, J.S. Murray, and P. Politzer, Halogen bonding: the sigma-hole, J. Mol. Model. 13, 291–296 (2007).
 G. Cavallo, P. Metrangolo, R. Milani, T. Pilati, A. Priimagi, G. Resnati, and G. Terraneo, The halogen bond, Chem. Rev. 116, 2478–2601 (2016).
 A.M.S. Riel, R.K. Rowe, E.N. Ho, A.C. Carlsson, A.K. Rappe, O.B. Berryman, and P.S. Ho, Hydrogen bond enhanced halogen bonds: a synergistic interaction in chemistry and biochemistry, Acc. Chem. Res. 52, 2870–2880 (2019).
 R. Gatti, J.M. MacLeod, J.A. Lipton-Duffin, A.G. Moiseev, D.F. Perepichka, and F. Rosei, Substrate, molecular structure, and solvent effects in 2D self-assembly via hydrogen and halogen bonding, J. Phys. Chem. C 118, 25505–25516 (2014).
 A. Okada, S. Hasui, K. Mino, S. Hara, M. Yoshimura, and K. Kadono, Structural variations of two-dimensional networks of 2,4,6-tris(4-bromophenyl)-1,3,5-triazine on Au(111) surface from solutions, Surf. Sci. 702, 121718 (2020).
 F. De Marchi, G. Galeotti, M. Simenas, M.S. Gallagher, E. Hamzehpoor, O. MacLean, R.M. Rao, Y. Chen, D. Dettmann, G. Contini, et al., Temperature-induced molecular reorganization on Au(111) driven by oligomeric defects, Nanoscale 11, 19468–19476 (2019).
 W. Zhao, L. Dong, C. Huang, Z.M. Win, and N. Lin, Cu- and Pd-catalyzed Ullmann reaction on a hexagonal boron nitride layer, Chem. Commun. 52, 13225–13228 (2016).
 Z. Yang, L. Fromm, T. Sander, J. Gebhardt, T.A. Schaub, A. Görling, M. Kivala, and S. Maier, On-surface assembly of hydrogen and halogen-bonded supramolecular graphyne-like networks, Angew. Chem. Int. Ed. 59, 9549–9555 (2020).
 T.A. Pham, F. Song, M.-T. Nguyen, and M. Stöhr, Self-assembly of pyrene derivatives on Au(111): substituent effects on intermolecular interactions, Chem. Commun. 50, 14089–14092 (2014).
 T.A. Pham, F. Song, M.-T. Nguyen, Z. Li, F. Studener, and M. Stöhr, Comparing Ullmann coupling on noble metal surfaces: on-surface polymerization of 1,3,6,8-tetrabromopyrene on Cu(111) and Au(111), Chem. Eur. J. 22, 5937–5944 (2016).
 M. Lischka, M. Fritton, J. Eichhorn, V.S. Vyas, T. Strunskus, B.V. Lotsch, J. Björk, W.M. Heckl, and M. Lackinger, On-surface polymerization of 1,6-dibromo-3,8-diiodpyrene – comparative study on Au(111) versus Ag(111) by STM, XPS, and NEXAFS, J. Phys. Chem. C 122, 5967–5977 (2018).
 M. Lischka, G.S. Michelitsch, N. Martsinovich, J. Eichhorn, A. Rastgoo-Lahrood, T. Strunskus, R. Breuer, K. Reuter, M. Schmittel, and M. Lackinger, Remote functionalization in surface-assisted dehalogenation by conformational mechanics: organometallic self-assembly of 3,3’,5,5’-tetrabromo-2,2’,4,4’,6,6’-hexafluorobiphenyl on Ag(111), Nanoscale 10, 12035–12044 (2018).
 J. Hu, J. Hu, H. Wang, K. Shen, H. Zhang, C. Huang, L. Xie, Q. Tian, H. Huang, Z. Jiang, and F. Song, Initiating Ullmannlike coupling of Br2Py by a semimetal surface, Sci. Rep. 11, 3414 (2021).
 A. Ibenskas, M. Šimėnas, and E.E. Tornau, Multiorientation model for planar ordering of trimesic acid molecules, J. Phys. Chem. C 122, 7344–7352 (2018).
 A. Ibenskas and E.E. Tornau, Modeling of high-temperature ordered structures with weak intermolecular C–H···F and C–H···N bonds, J. Phys. Chem. C 125, 19560–19569 (2021).
 J. Lisiecki and P. Szabelski, Surface-confined metal-organic precursors comprising naphthalene-like derivatives with differently distributed halogen substituents: A Monte Carlo model, J. Phys. Chem. C 124, 20280–20293 (2020).
 K. Nieckarz, P. Szabelski, and D. Nieckarz, Monte Carlo simulations of the self-assembly of hierarchically organized metal-organic networks on solid surfaces, Surf. Sci. 719, 122041 (2022).
 D.K. Jacquelín, F.A. Soria, P.A. Paredes-Olivera, and E.M. Patrito, Reactive force field-based molecular dynamics simulations on the thermal stability of trimesic acid on graphene: implications for the design of supramolecular networks, ACS Appl. Nano Mater. 4, 9241–9253 (2021).
 N. Kalashnyk, M.D. Mansour, J. Pijeat, R. Plamont, X. Bouju, T.S. Balaban, S. Campidelli, L. Masson, and S. Clair, Edge-on self-assembly of tetra-bromoanthracenyl-porphyrin on silver surfaces, J. Phys. Chem. C 124, 22137–22142 (2020).
 A. Ibenskas and E.E. Tornau, Modeling of different ordering schemes for halogen-functionalized molecules with triazine and benzene core, J. Phys. Chem. C 126, 8079–8089 (2022).
 F. De Marchi, G. Galeotti, M. Simenas, E.E. Tornau, A. Pezzella, J. MacLeod, M. Ebrahimi, and F. Rosei, Room-temperature surface-assisted reactivity of a melanin precursor: silver metal–organic coordination versus covalent dimerization on gold, Nanoscale 10, 16721–16729 (2018).
 F. Neese, Software update: the ORCA program system – Version 5.0, WIREs Comput. Mol. Sci. 8, e1606 (2022).
 A. Najibi and L. Goerigk, DFT-D4 counterparts of leading meta-generalized-gradient approximation and hybrid density functionals for energetics and geometries, J. Comput. Chem. 41, 2562–2572 (2020).
 E. Caldeweyher, C. Bannwarth, and S. Grimme, Extension of the D3 dispersion coefficient model, J. Chem. Phys. 147, 034112 (2017).
 E. van Lenthe, J.G. Snijders, and E.J. Baerends, The zero‐order regular approximation for relativistic effects: the effect of spin-orbit coupling in closed shell molecules, J. Chem. Phys. 105, 6505–6516 (1996).
1,3,6,8-tetrabromopireno (Br4Py) 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.