Coronene Molecules in Hexagonal Pores of Tricarboxylic Acids: A

Aug 14, 2015 - We find two structures formed by dimeric and trimeric interactions which can contain coronene in their pores and mix in different propo...
0 downloads 12 Views 3MB Size
Subscriber access provided by TEXAS A&M INTL UNIV

Article

Coronene Molecules in Hexagonal Pores of Tricarboxylic Acids: A Monte Carlo Study Mantas Šim#nas, Andrius Ibenskas, and Evaldas E. Tornau J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.5b06690 • Publication Date (Web): 14 Aug 2015 Downloaded from http://pubs.acs.org on August 20, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Coronene Molecules in Hexagonal Pores of Tricarboxylic Acids: a Monte Carlo Study ˇ enas,† Andrius Ibenskas,‡ and Evaldas E. Tornau∗,‡ Mantas Sim˙ Faculty of Physics, Vilnius University, Sauletekio 9, LT-10222 Vilnius, Lithuania, and Semiconductor Physics Institute, Center for Physical Sciences and Technology, A. Goˇstauto 11, LT-01108 Vilnius, Lithuania E-mail: [email protected] Phone: +370 5 2616915. Fax: +370 5 2627123

Abstract We propose a 7-state (orientation) model with some exclusions to describe the ordering of symmetry-reduced biphenyl-3,4’,5-tricarboxylic acid (BHTC) molecules into planar ordered structures on a substrate lattice. The model is based on the 3-state Bell-Lavis model, which is used to describe the symmetric triangular molecules. The reduction of molecular symmetry allows us to predict different honeycomb-like structures H-bonded by dimeric and trimeric interactions. The model is supplemented by the guest-host interaction term which enables to analyze the effect of coronene on ordering of molecular structures. We find two structures formed by dimeric and trimeric interactions which can contain coronene in their pores and mix in different proportions depending on coronene concentration. The effect of host-coronene interaction on ordering of trimesic (TMA) and benzene-tribenzoic (BTB) acids is also studied using the ∗

To whom correspondence should be addressed Faculty of Physics, Vilnius University, Sauletekio 9, LT-10222 Vilnius, Lithuania ‡ Semiconductor Physics Institute, Center for Physical Sciences and Technology, A. Goˇstauto 11, LT-01108 Vilnius, Lithuania †

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3-state Bell-Lavis model. We demonstrate that the interaction with coronene might induce the co-existence of coronene-filled honeycomb phase and coronene-free superflower phase in the TMA molecular system. We also predict that this interaction might induce the occurrence of superflower phase in the system of BTB molecules.

Introduction Supramolecular chemistry is extensively used for targeted design of planar self-assembled supramolecular formations on solid surfaces. In recent years many diverse surface structures have been prepared and quite a large part of them comprises nicely ordered H-bonded molecular networks. The directionality, chemical selectivity and flexibility of H-bonds allow to organize open molecular networks with well-defined pores used for confinement of guest molecules with different functional properties. 1–4 Many important prerequisites have to be taken into consideration to organize a planar molecular structure. The electronic 5 and symmetry (surface commensurability with molecular network) properties of substrate template are important. The choice of correct properties and proportions of solute and solvent (see 6 and Refs therein) plays the role when the selfassemblies occur at the solid-liquid interface. Thermodynamic conditions are important as well as concentration of solute (or molecular packing density in vacuum experiments 7 ), since by tuning the molecular density the same molecular building block might be used to build different ordered planar molecular structures. Nevertheless, the leading role in building supramolecular entities is played by a detailed knowledge of intermolecular and substrate-molecule interactions. 1–4 In natural conditions a dimeric intermolecular H-bond interaction, featuring a high degree of directionality (with preferred bond angles being 170◦ -180◦) and strong covalency , 8 is established. The symmetric cyclic trimeric arrangement is also quite often observed in ordering of triangular molecules. Along with experimental (especially STM) measurements, the theoretical modeling is necessary to better understand the molecular ordering processes and to explain and predict 2

ACS Paragon Plus Environment

Page 2 of 34

Page 3 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

new surface structures. For such a calculation the experimental findings by the STM are used as a starting material to construct the main model with an appropriate set of chosen interactions. As a result, the possible patterns might be found without performing actual experiments. Moreover, the modeling might be really helpful for interpretation of superstructure geometry of STM snapshots obtained with submolecular resolution. Theoretical studies of planar H-bonded molecular orderings are performed using ab initio, 9–11 molecular mechanics 12,13 and molecular dynamics 14,15 calculations or combination of these studies with Monte Carlo (MC) calculations. 14,16–18 Molecular ordering is also modeled using statistical phase transition models, most frequently - lattice models. 16–23 In such calculations a molecule is represented by a simple geometrical form with a minimal number of molecular orientations. In addition, only the most characteristic molecular interactions are chosen which have to reveal the main orderings of the model. The models are sometimes supplemented by some exclusion rules. 23,24 Here we present a 7-state model and MC calculation for ordering of triangular molecules of reduced symmetry. This model can be easily reduced to the Bell-Lavis model 25–27 which is used to describe 21 the self-assembly of symmetric triangular molecules. This paper is mostly devoted to reveal the role played by guest (coronene) molecules in ordering of tricarboxylic acid molecules. We consider both symmetry-reduced molecules, such as biphenyl-3,4’,5tricarboxylic acid (BHTC), and symmetric molecules, such as trimesic acid (TMA) and 1,3,5-benzene-tribenzoic acid (BTB) (Fig. 1). The results are presented here not in a form of rigorous phase diagrams, but in a simple form of structure snapshots more appropriate for prediction of results and analysis of experimental data. Both TMA and BTB molecules are three-fold-symmetric hydrogen-bonding units. The TMA molecule has one (central) phenyl ring and three carboxylic groups along each of the three bonds of the molecule. The BTB molecule has four phenyl rings, central and three peripheral ones, with carboxylic group originating from each perypheral phenyl ring. BHTC has two carboxylic groups originating from the central phenyl ring as in TMA molecule.

3

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1: Schematic representation of (a) TMA, (b) BTB, (c) BHTC and (d) coronene molecules. Large (small) grey circles denote carbon (hydrogen) atoms, red circles - oxygen atoms. The third carboxylic group is on a bond elongated by a peripheral phenyl ring as in BTB molecule. The symmetry of BHTC is lower than their better known counterparts, therefore the ordering possibilities are more diverse. The BHTC molecules demonstrate very interesting self-assemblies at the solid-liquid interface. 28 The dimeric H-bond interactions, analogous to those observed in honeycomb orderings of TMA and BTB molecules, play a crucial role in formation of ordered porous BHTC structures. Along with the dimeric bonds, a trimeric H-bonding motif similar to that found for TMA flower structure 7,29–31 is experimentally observed in BHTC system 28 when the coronene molecules are added. The modeling of self-assembly of TMA 21,22 and melamine 16,24 molecules (both possessing C3 symmetry) was performed by mapping the shape of these molecules onto an equilateral triangle with interacting vertices (TMA) or sides (melamine). The BTB molecule might be mapped onto a larger equilateral triangle with interacting vertices. Following the same logic, the BHTC molecule, which has a reduced C2 symmetry, might be mapped onto an isosceles triangle with interacting vertices. The fact that the radius of the BTB molecule is twice that of the TMA molecule can be used as an advantage: the BHTC molecule might be nicely incorporated on a substrate with a central phenyl ring located on a site of a substrate lattice in such a way that interactions by carboxyl groups are possible on both short and long bonds 4

ACS Paragon Plus Environment

Page 4 of 34

Page 5 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

of the molecule. The decrease of molecular symmetry from equilateral to isosceles triangle leads to a very rich phase diagram with a large number of ordered structures. All these phases are porous honeycomb-like formations, some of them demonstrate the nets of slightly irregular hexagonal pores. The occurrence of structures with trimeric bonds is one of the most interesting problems encountered in studies of the triangular molecules interacting by their vertices. The trimeric bonds are found either at higher molecular density 7,29–31 or when foreign molecules are incorporated in pores of a host molecules matrix. 28 In a former case they support the formation of flower phases of TMA molecules. In a latter case the guest (coronene molecule) stimulates the occurrence of hexagonal porous formations in BHTC system which must possess the trimeric bonds in order to build a symmetric structure with similar pore sizes. Actually, it is not clear if the increase of a density due to a guest molecule insertion or the host-guest interaction itself is the main reason of a trimeric arrangement. The ratio of a trimeric to dimeric interaction (per molecule) in DFT calculations for tricarboxylic acids is more or less the same (0.8-0.9

28,31

). As shown by our studies here, such a value of trimeric interaction is

sufficient to independently create the phases with trimeric bonds at high molecular densities for TMA molecular system, while for molecules with lower conformational flexibility (BHTC and BTB) some support of host-guest interaction is needed. It should be noted, that our modeling has some similarity with a well-known ”tiling problem“ realized experimentally by Blunt et al 32 for a two-dimensional arrangement of pterphenyl-3,5,3”,5” -tetracarboxylic acid (TPTC) on graphite. The authors demonstrated that the unique honeycomb ordering might be achieved by 5 different assemblies of TPTC molecule mapped onto a rhombic tile. The STM interpretation becomes cumbersome, kind of a ”bird’s eye view”, and the final structure is a likely mixture of possible assemblies or a spin-glass (from the point of view of phase transition problem). Here, in our simulations for BHTC molecules, we also find several honeycomb structures demonstrating same density and symmetry of pores which might be assembled in different ways, the details of the assem-

5

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

bly being, probably, hardly distinguishable by the STM. Especially clearly this problem is manifested while studying experimentally observed 28 low density porous coronene-stabilized phase. We found that the pores of the same density might be formed by two similar, but yet slightly different molecular orderings. Our results imply that an interpretation of the experimental structure as a mixture of the two phases should be more feasible. The paper is organized in a following way. The model for ordering of symmetry-reduced molecules, main interaction parameters and methods of calculation are presented in Sec. 2. The ground state calculations for the BHTC system without coronene are given in Sec. 3.1. The effect of coronene on lower symmetry (BHTC) and higher symmetry (TMA and BTB) molecular systems is presented in Sec. 3.2 and 3.3, respectively.

Figure 2: States (orientations), interactions and exclusions of the model: (a) six molecular orientations, (b) main pair and triangular H-bond interactions (LL stands for long-long; SL short-long; SS1 (SS2)-short-short interactions of the first (second)-type; T -trio interactions) of the BHTC molecules. Lower right drawing depicts the exclusions for the molecule in state 1. The black solid circles show the exclusion for the center of the molecule in any state; the upper and lower crosses show the exclusions for molecules in neighboring states 2 and 6, respectively (the corresponding exclusion rule is used for other molecular states). (c) Two examples of the BHTC-coronene interaction. Both interactions are considered to be equal in magnitude.

6

ACS Paragon Plus Environment

Page 6 of 34

Page 7 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Model and methods of simulation To describe the ordering of the BHTC molecules, we use the lattice model on a triangular lattice as an underlayer. The lattice is chosen for convenience and simplicity of calculations. The molecule is assumed to be “anchored” to the site of a substrate lattice in such a way that the central phenyl ring, which is bonded to peripheral ring and two carboxylic groups, has the on-site location. While moving from site to site, the molecule might have any of six orientations on each site with the longer phenyl arm being kind of a “vector” pointing along each of the lattice directions. The length of such a vector is ∼ 2a, where a is the triangular lattice constant. The orientations of the molecule are shown in Fig. 2a. Our model comprises 6 molecular states and a vacancy state which occurs when the site is not occupied by the molecule. Thus, in general, the model is a 7-state model with interacting vertices (carboxylic groups) and some exclusion rules dictated by the size and form of the molecule. Since the length of a longer arm (distance from the center of the central phenyl ring to the carboxylic group) is roughly twice as that of a shorter arm, the main H-bond interactions existing in such a system might be conveniently taken into account. Here we employ the Bell-Lavis model 25–27 previously used to describe the ordering in twodimensional bonded fluid. The lattice-gas variable ni = 1(0) describes the existence (nonexistence) of the center of the BHTC molecule on a site i, while the variable τiij characterizes the orientational ordering: if the H-bond exists between the molecules at sites i and j, τiij = 1, if - not, τiij = 0. Thus, the Hamiltonian (per molecule) for BHTC molecules has the form

HBL = −

1X 1X ni nj nk eijk τiijk τjjik τkkji , ni nj eij τiij τjji − 2 i,j 3 i,j,k

(1)

where eij and eijk are dimeric and trimeric interactions shown in Fig. 2b and coefficients 1/2 and 1/3 are taken to avoid the double and triple counting in first and second sum, respectively.

7

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The dimeric H-bond interactions eij might be of a different type: the long-long interaction, ell , representing the interaction by longer arms of the molecules when the centers of interacting molecules are 4a apart; the short-long interaction esl when the molecules are 3a apart and the two short-short interactions, ess1 and ess2 , when the molecules are 2a apart and their “vectors” make the mutual angles of 180 and 60 degrees, respectively. There are no other distances and angles appropriate for the dimeric H-bond to be established. There is also a trimeric interaction in the model, eijk = et , which is possible, if three √ molecules exist on vertices of an equilateral triangle with a side length 2a 3 and their longer arms are oriented towards mutual center creating a trimeric H-bond (Fig. 2b), i.e. all variables τiijk = τjjik = τkkji = 1, where τiijk = τiij τiik , etc. Note, that there might be other trimeric interactions in the model. However, they are not considered here, since the structures with such molecular configurations are not experimentally observed and addition of these interactions would significantly complicate our model. To avoid the overlap of molecules and occurrence of non-realistic structures, we also use some exclusion rules which are shown for a molecule in the state 1 in a bottom right drawing of Fig. 2b. Here the black dots correspond to sites which cannot be occupied by centers of other molecules, while the crosses show the sites which are forbidden for centers of molecules in states 2 (upper cross) and 6 (lower cross) neighboring to the state 1 (to avoid the structure with overlapping carboxylic groups on a longer arm of three neighboring molecules). For molecules in other states we use the same kind of exclusions. It should be noted that, most likely, some other exclusions could be taken in the model, but they are important only at high densities of molecules which can hardly be achieved in real experiments. When guest (coronene) molecules are inserted into the BHTC system, we should account for the interaction, ecor , stimulating the occurrence of the ordered BHTC-coronene systems. We assume that this interaction acts between the coronene molecule and BHTC molecule (in any of its six orientations) when the distance between their centers is 2a (Fig. 2c). The resulting energy with this interaction taken into account is

8

ACS Paragon Plus Environment

Page 8 of 34

Page 9 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

HBLc = HBL −

X

ni tj ecor ,

(2)

i,j

where tj is the lattice-gas variable which represents the existence (tj = 1) or absence (tj = 0) of the center of coronene molecule on the lattice site j. It should be noted, that TMA and BTB (Fig. 1) molecules, which can be mapped on equilateral triangles, might also be described by the Bell-Lavis model in the form (1) and (2). However, for these high-symmetry molecules the appropriate model is simpler: a molecule in a site of a triangular lattice might have any of the three states - one of two orientations, which differ by 60 degrees rotation of the molecule, or a vacancy state. The dimeric H-bond interaction ed in the 3-state Bell-Lavis model can exist (τij = 1 in (1)) at a distance 2a, if the ordering of TMA molecules is considered. For BTB molecules this distance is 4a. The trimeric interaction exists if three molecules are arranged on vertices of an equilateral √ √ triangle with a side length a 3 (for TMA molecules) and 2a 3 (for BTB molecules). We performed the MC simulation of the proposed 7-and 3-state models using Metropolis algorithm and Kawasaki and Glauber dynamics. 33 When we used the latter dynamics, we considered the energy in the form

H = HBL + µ

X

ni ,

(3)

i

where µ is the chemical potential of studied molecules. The Glauber calculation starts when the lattice is randomly filled with molecules in different states corresponding to some fixed value of µ. The interaction energy Ei on a randomly chosen site is calculated. Then the initial state on that site is changed (with equal probability) to one of the six remaining states, and the final energy Ef is calculated. Then the Metropolis procedure is performed: the new state is accepted, if ∆E = Ef − Ei < 0, or accepted with the probability ∼ exp(−∆E/kB T ), if ∆E > 0. Note that if the forbidden state is chosen then the Metropolis procedure is interrupted and a new molecule is randomly selected. In a process of Glauber dynamics, the

9

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

molecular concentration is changed, while the chemical potential is fixed. For calculations with coronene we used the Kawasaki dynamics and Hamiltonian (2). Employing this method, the diffusion of molecules is simulated for fixed value of molecular concentration. To decrease the computation time and avoid freezing of the system in metastable states, we employed the version in which the molecule can “jump” to any nonoccupied site at arbitrary distance. Our simulations start from disordered random mixture of tricarboxylic acids and coronene. The algorithm is realized by randomly choosing the molecule at initial site (with energy Ei ) and the empty site, where the molecule is supposed to jump. One of the six states of the molecule is randomly chosen for the empty site, and corresponding final energy Ef is calculated. Then the Metropolis procedure is performed: the new state is accepted, if ∆E < 0, or accepted with the probability ∼ exp(−∆E/kB T ), if ∆E > 0. For MC calculations with Kawasaki dynamics we used the triangular lattices of sizes L × L commensurate with expected molecular structures. Mostly, we used L = 72 with periodic boundary conditions. We probed our model for higher values of L, but the obtained results did not differ. Our simulations were performed starting from higher temperature and using random initial particle configuration. Then the temperature was lowered below the order-disorder phase transition point which depends on molecular concentration. Then the temperature was further slightly decreased to reduce thermal fluctuations and obtain fully ordered structure. The total number of MC steps at lowest temperature was at least 107 . Our calculations of the H-bond dimeric interactions using density functional theory (DFT) and accounting for the van der Waals (vdW) interactions and BSSE correction showed that ess1 & esl & ell ≈ ess2 , and the magnitudes of the dimeric interactions differ by about 1-1.5% (while the absolute values are in between 9.62 and 9.49 kcal/mol per molecule). Our calculation demonstrates some differences from the former calculations, 28 in which the vdW contribution and BSSE correction were neglected. However, both calculations were performed for the gas phase and did not account for some very important effects (direct effects

10

ACS Paragon Plus Environment

Page 10 of 34

Page 11 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

of the substrate lattice and indirect effects due to free atoms of the substrate 34 ) which might change the value of interactions. We assume that the difference (∼ 5%) in values obtained by two DFT calculations for a gas phase might be smaller than that given by accounting for the substrate lattice. 9 Since in our MC calculations the molecule-substrate interaction is neglected, we left a larger frame for variation of dimeric interaction constants. This allows us, by slight tuning of the intermolecular interactions, to predict the possibility of some orderings not seen in experiments. We fixed the values ell = esl as a constant ed to which other constants and temperature are rescaled to dimensionless form. Note, that the ratio of trimeric interaction to the largest dimeric interaction in our DFT calculations was 0.86 (per molecule) in good accordance with previous results. 28 We approximately kept this ratio in our MC calculations.

Figure 3: Arrangement of molecules in phases F2 and F4 (c = 3/32), F1 and F3 (c = 2/21) and HG (c = 1/10) and their representation in MC simulations. The contours of occurring hexagonal pores are highlighted in red.

11

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 34

Results Ground state calculations, ecor = 0 We used MC calculations to determine possible ordered structures occurring in our model for different values of interaction parameters and calculated the ground state (GS) energies of these structures. For further studies we have chosen the phases which emerged for different BHTC molecular concentrations when the difference in magnitude of dimeric interaction parameters was within ∼ 20%. The first group of phases consists of four flower-type hexagonal structures shown in Fig. 3 which might be observed at lower molecular densities. To describe the density, here we use the concentration, c. The concentration of any particular phase is obtained by calculating the number of occupied by molecules lattice sites and dividing this number by L2 . In our calculations we assume that a molecule occupies only one lattice point. For example, a BHTC molecule occupies a site which is under the central phenyl ring. The GS energies of the flower phases are

EF1 = −cF1 (esl + ess1 /2 − µ),

(4)

EF2 = −cF2 (esl + ess1 /2 − µ), EF3 = −cF3 (esl + ess2 /2 − µ), EF4 = −cF4 (ess2 + ell /2 − µ), cF1 = cF3 = 2/21 ≈ 0.095, cF2 = cF4 = 3/32 ≈ 0.094. It is interesting to note that the structures F1 and F2 have the same energy and just slightly different concentrations. The F1-phase features one type of hexagon with two side

12

ACS Paragon Plus Environment

Page 13 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

lengths equal to 2a and four side lengths - 3a. The F2-phase demonstrates two types of hexagons: one - regular (side length equal to 3a) and the other - irregular (three sides have length 2a and three - 3a). In GS phase diagrams only F1 phase is emerging due to denser packing. However, in finite temperature MC calculations both F1 and F2 phases use to occur, most likely, due to higher entropy of the F2 structure. The low-density flower structure F3 has the same geometry and concentration as the F1-phase. The pores of both these phases consist of four longer sides and two shorter sides of a schematic BHTC triangle, just the shorter sides of the F1-phase are on opposite sides of a hexagonal pore, while in F3-phase they are close to each other. The most exotic is the F4-phase which consists of central regular hexagonal pore, created by shorter sides of six BHTC molecules. This pore is surrounded by six irregular hexagons containing three long (4a) and three short (2a) sides (Fig. 3). There is one additional structure which has low molecular density and therefore might be attributed to the flower phases. We call it the hour-glass (HG) phase (cHG = 0.1, see Fig. 3). The fragments of the HG are experimentally found 28 at higher molecular densities in between large domains of a high-density zigzag structure with cZ = 0.143. The relation between the HG and zigzag phases is a subject of a separate topic which is planned to be published elsewhere. The range of concentration makes the HG phase a likely competitor of the flower phases in the phase diagram. The GS energy of the HG structure is

EHG = −cHG (ell /2 + ess1 − µ).

(5)

The second group of phases consists of two higher-density phases (Fig. 4) with trio H-bond interactions. One of these structures, C1 (cC1 = 0.125)), was suggested 28 when interpreting the STM images of the BHTC-coronene ordering - in this structure the coronene molecules fill the hexagonal pores formed by the BHTC molecules. The C1 phase is a sequence of two hexagons (both have side length equal to 2a). The sides of a first hexagon are assembled of six shorter sides of a schematic BHTC triangle. The sides of a second 13

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 34

hexagon are made of four longer sides of such a triangle. Further we call these hexagons “six-molecule” and “four molecule” pores. The ratio of number of six- to four-molecule pores in the C1 structure is 1:3. Performing MC simulations we noticed another very similar structure possessing the same molecular density, as well as form and size of hexagonal pores, as C1. We called this phase C2 (cC2 = 0.125). It also demonstrates a sequence of two hexagons (both have side length equal to 2a): one is formed of four longer sides of a schematic BHTC triangle as a four-molecule pore of C1, but the other hexagon is assembled of five sides of the triangle (two longer sides and then three shorter sides in a row) comprising “five-molecule” pore. There is the same number of five- and four-molecule pores in the C2 structure. The GS energies of these two phases are

EC1 = −cC1 (ess2 + et /3 − µ),

(6)

EC2 = −cC2 (ess1 /3 + 2ess2 /3 + et /3 − µ), i.e. the energies of these two phases are equal when both short-short interactions are equal, ess1 = ess2 . Despite diversity of obtained phases and dimeric interactions, both low-density flower (F) phases (4) including the HG phase (5), i.e the phases with c ≤ 0.1, and trimer-organizing phases C1 and C2 (further denoted as C) with c = 0.125 can be analyzed using GS phase diagrams. Let us initially neglect a small difference in magnitude of dimeric bond interactions, ell = esl = ess1 = ess2 = e. Then the GS energies of the F and C phases are EF = −cF (3e/2 − µ) and EC = −cC (e + et /3 − µ), respectively. Here as cF we take cHG as the densest of the first group of phases. The resulting phase diagram is shown in Fig. 5, where the line 5et = 3e + 3µ separates the F and C phases. It should be noted, that in this diagram only the interval of 1.2 < µ/e < 1.5 is interesting, because its upper bound is the gas phase limit, while at µ/e < 1.2 the finite temperature MC calculations demonstrate the molecular

14

ACS Paragon Plus Environment

Page 15 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 4: Arrangement of molecules in phases C1 and C2 (c = 1/8) and their representation in MC simulations. Two hexagons with the six- and four-molecule pores in phase C1 and five- and four-molecule pores in phase C2 are highlighted in red.

15

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

density which is too large for any of the C phases to be obtained. We have chosen two points in this phase diagram at optimal value of chemical potential, µ/e = 1.3, with rather small (et /e = 1) and rather high (et /e = 1.8) trio interaction to compare with the DFT result. Note, that the DFT ratio of trimeric- to dimeric interaction (0.86 per molecule) is equal to et /e ≈ 1.3 in our MC calculations. At et /e = 1 and 1.8 we perform further calculation of the GS phase diagrams using formulas (4)-(6), assuming ell = esl = ed and varying ess1 /ed and ess2 /ed in small limits (0.9 ÷ 1.1). The obtained diagrams are shown as insets to Fig. 5. The main conclusions are the following: (1) the magnitude of et has to be rather large for the C phases to occur: the C1 phase exists at ess1 < ess2 and the C2 - at ess1 > ess2 ; (2) roughly the same inequality separate the F4 and F3 phases on one side and F1 and HG phases on the other for not too large et ; (3) the DFT result for et implies that this value is too small to obtain the C1 or C2 phase without additional interactions. Note, that the structures, C1 and C2, are not experimentally observed without coronene. 28 Formally, the phases C1 and C2 cannot be formed simultaneously - it is the ratio of interaction constants, ess2 /ess1 which predetermines which phase would be formed. The gas phase DFT calculations provide us with very similar magnitudes of dimer interactions, ess2 /ess1 = 1.01 28 and 0.97 (our calculations). Given the similarity of dimeric interactions in DFT calculations, the effect of substrate and free surface atoms might easily shift this ratio in favour of one of the phases. Therefore, the conclusions of the GS analysis, certainly, do not deny the situation when both phases exist as domains in different regions of the lattice or when they mix. The situation with mixing of these two phases is very similar to that found for a rhombic tiling 32 when one type of honeycomb arrangement could be achieved by 5 different assemblies of TPTC molecules on graphite, and the resulting structure is a kind of molecular “spinglass“. In our case the structure (which might seem unique in STM measurement) might be obtained by two different assemblies, and be either pure phase or their “spin-glass“ mixture.

16

ACS Paragon Plus Environment

Page 16 of 34

Page 17 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 5: Ground state (µ, et ) phase diagram. Insets are (ess1 , ess2 ) diagrams at two points of main diagram µ/ed = 1.3 and et /ed = 1.0 and 1.8. The cross and dashed-line mark the DFT result. 28

Effect of coronene on occurrence of C1 and C2 phases Our GS estimate and following MC calculations demonstrate that the trio interaction has to be at least of the same magnitude (per molecule) as the main dimer interactions in order to form the C1 or C2 phases. The DFT calculations 28,31 show that the magnitude of et is lower. This is the most probable reason why the C1 or C2 phases do not occur in the absence of coronene. 28 Meanwhile, it was experimentally shown 35 that the coronene can be easily incorporated into the honeycomb (six-molecule pore) phase of TMA molecules. It increases the binding energy of the system by 3.9-4.6 kcal/mol (DFT result 36 ). In BHTC-coronene system the coronene plays the same role of a binding element which effectively strengthens the trio interaction and, consequently, leads to formation of the C1 and C2 phases. Here we demonstrate how the coronene-BHTC interaction, ecor , affects the occurrence of these phases when the concentration of coronene ccor is varied. It is very interesting that the variation of ccor might lead to switching from one phase to another: for low values of ccor and ess1 ≈ ess2 the C1 phase is preferred due to six-molecule surrounding of the coronene, while at high value of ccor which phase would appear depends on the ratio of interaction constants ess1 /ess2 . 17

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 34

Here we perform the Kawasaki-type treatment and therefore neglect the effect of chemical potential. We also introduce the parameter α = ccor /cC1 which is the ratio of coronene coverage to BHTC coverage in any of the C1 or C2 phases (cC1 = cC2 ). When hexagonal pores of BHTC are fully occupied by the coronene in both these phases (α = 2/3), their GS energy might be written

2/3

EC1 = EC1 − 3cC1 ecor ,

(7)

2/3

EC2 = EC2 − 3cC2 ecor . It is seen that at full occupancy the term related to interaction ecor brings the same energy contribution to the C1 and C2 phases, despite the difference of coronene surroundings. This difference is important when the parameter α (concentration of coronene ccor ) decreases obviously the coronene is best bonded in six-molecule pores and worse - in four-molecule pores. This preference leads to occurrence of the C1 phase with coronene at the six-molecule centers at α = 1/6 and the C2 phase with coronene chains in five-molecule pores at higher values of α. Taking into account that the decrease of α leads to gradual reduction of coronene firstly in four-molecule pores, then - in five-molecule-pores and finally - in six-molecule pores, we can write the energy in more general form

1 α EC1 = EC1 − ( + 4α)cC1 ecor , 3 ≤1/6

EC1

= EC1 − 6αcC1 ecor ,

1 α EC2 = EC2 − ( + 4α)cC1 ecor , 3 ≤1/3

EC2

= EC2 − 5αcC1 ecor ,

1/6 ≤ α ≤ 2/3

(8)

α ≤ 1/6 1/3 ≤ α ≤ 2/3 α ≤ 1/3.

At 1/3 ≤ α ≤ 2/3, i.e. until all coronene vanishes from the four-molecule pores of the C2 phase, the energies of both phases are equal, if the main dimer interactions in C1 and

18

ACS Paragon Plus Environment

Page 19 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

C2 phases are equal (ess1 = ess2 and, consequently, EC1 = EC2 ). With further decrease of α ≤ 1/3 the coronene leaves four-molecule pores in phase C1, and five molecule-pores in phase C2. As a result, the energy of C2 phase relatively increases, and the C1 phase becomes more preferred at α ≤ 1/3 (see Fig. 6). At α = 1/6 the C1 phase completely dominates, because the number of coronene molecules and the number of six-molecule pores coincides. It is interesting that the C1 phase is more favourable (just for slightly decreased interval of α < 1/3) also for ess1 & ess2 (see line 1 in Fig. 6 and its inset), i.e. for such a ratio of dimeric interactions which favours the C2 phase both at high coronene concentrations (7) and in absence of coronene molecules (6), but large et .

α α Figure 6: (a) Ground state energies of phases C1 and C2 and their difference ∆E = EC1 −EC2 (inset) as a funcion of α = ccor /cC1 at ess1 /ed = 1.2 > ess2/ed = 1.05 (1), ess1/ed = ess2 /ed = 1/6 1.1 (2) and ess1 /ed = 1.05 < ess2 /ed = 1.17 (3). Red circles mark the energy of EC1 and 1/3 EC2 mixing at α = 1/4 (see text). For convenience, the lines 1 are shifted up and the lines 3 are shifted down along the E-axis by 0.04. (b) The relative number of six-molecule N6 (solid curves) and five-molecule N5 (dashed curves) pores as a function of α obtained in MC calculations. The denotation of curves and values of interaction parameters are the same as in (a). The data is rescaled to dimensionless form with respect to their values of in ideal phases C1(N6 ) and C2 (N5 ) at α = 2/3. The data points are averages over 4-6 runs. Lines are guides to the eye. Here et /ed = 1.4 and ecor /ed = 0.6.

It should be noted that we seldom obtain pure structures in our MC calculations. The analysis of the C1 and C2 phases demonstrates that they can mix very easily due to a common building block of both phases - zigzag chains of four-molecule pores (Fig. 4). These chains are the main ”sticking“ element which easily combines pieces of domains of both 19

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 34

phases in any proportions without energy losses due to domain walls. These different pieces of domains are six-molecule stars in C1 phase and zigzag five-molecule chains in C2 phase . There are several formal reasons for mixed structures to be seen in our numerical experiment. The entropy of a mixed phase at finite temperature can be higher than the entropy of any of the pure phases. Besides, it might be also shown from purely energetic consideration that the mix of the C1 and C2 phases for some values of α can be energetically more favourable than any of the pure phases, C1 or C2. Consider e. g. the coronene concentration at α = 1/4 for different values of ess1 and ess2 . Note, that the most appropriate combination of these two phases at α = 1/4 can be obtained by mixing C1 phase at stoichiometry αC1 = 1/6 (coronene in six-molecule pores) and C2 phase at stoichiometry αC2 = 1/3 (zigzag chains of fully coronene-occupied five-molecule pores) in equal proportions. For ess1 . ess2 the energy of such a mixing is higher than that of the lower-energy C1 phase; for ess1 = ess2 1/4

1/4

it is equal to EC1 , but for ess1 & ess2 the mix has lower energy than EC1 . This comparison is shown by red circles at α = 1/4 in Fig. 6. Our MC calculations confirm this result: we rather obtain the mixing of coronene-occupied six-molecule pores and coronene-occupied five-molecule zigzag chains than domains of pure phases C1 or C2 at α ≈ 1/4. In Fig. 7 we collected the snapshots at α = 2/3, 1/4 and 1/6 reflecting the most characteristic features of our calculations. To improve the visual perception of our snapshots, we marked the coronene molecules in six- and five-molecule pores by different color. Simultaneously, we calculated the share of each of the phases C1 and C2 in a mixture which might be found calculating the number of six- and five-molecules pores in each snapshot for different values of 1/6 ≤ α ≤ 2/3 (Fig. 6b). The curves were obtained by averaging over 4-6 runs. In Fig. 7 we show two snapshots of ideal C2 (Fig. 7a) and C1 (Fig. 7d) phases at α = 2/3. They were obtained for ess1 > ess2 and ess1 < ess2 , respectively. It should be noted, that the pure phases at α = 2/3 occur in about 30% of our simulations for such a choice of interaction parameters as in Figs. 6 and 7. We also demonstrate that with decrease of α the share of six-molecule pores increases both for ess1 > ess2 and ess1 ≤ ess2 , and phase mixing 20

ACS Paragon Plus Environment

Page 21 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

is rather prominent at α = 1/4 (see Fig. 6b and Fig 7b and e). Finally, we demonstrate the result of a competition between values of dimeric interactions (ess1 > ess2 ) favouring the C2 domain and value of coronene concentration (α = 1/6) favouring the C1 domain in Fig. 7c. In Fig. 7f we show the structure which is obtained when the former effect is neutral and the latter promotes the C1 domain. The structure is dominated by the six-molecule coronene-filled pores. For ess1 < ess2 and α = 1/6 the C1 domain is almost ideal.

Figure 7: Snapshots of MC simulations for different values of α and ess1 /ess2 ratio: (a) α = 2/3 (ideal C2 phase), (b) α = 1/4, (c) α = 1/6 for ess1 /ed = 1.2 > ess2 /ed = 1.05, (d) α = 2/3 (ideal C1 phase) for ess1 /ed = 1.05 < ess2 /ed = 1.17 and (e) α = 1/4, (c) α = 1/6 for ess1/ed = ess2 /ed = 1.1. For better visibility the coronene molecules are colored red, blue and gray in six-, five- and any-molecule pores, respectively. Here L = 72, et /ed = 1.4, and ecor /ed = 0.6. At α < 1/6 there is a lack of coronene to maintain the C1 phase. Therefore, some part of coronene-filled six-molecule pores are substituted by empty, but five-molecule pores. The C1+C2 structure is maintained until very small values of α. For very small coronene coverage (α ≈ 0.02) the system cannot maintain the C1+C2 structure anymore, but still rather large domains of this structure might be found around each coronene molecule. These domains are connected by rather disordered flower-phase and HG hexagons. The fact that so small a number of coronene molecules maintain so large domains of the C1+C2 phases 21

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

might also be related with our Kawasaki dynamics, where the jump event occurs between non-neighbouring sites. This treatment might create small patches of the C1+C2 phases over a whole lattice, thus artificially maintaining the overall C1+C2 phase. On the other hand, it cannot be done if the energy of the C1+C2 mixture and the corresponding F phase would not be rather similar. At α → 0 the C1+C2 domains are replaced either by the F-phases or a system of disordered hexagons. To facilitate the ordering of our system, we used rather high values of BHTC-coronene interaction, ecor /ed = 0.6. To find, how the ordering depends on this interaction, we performed an extra analysis varying ecor . It turned out that the results do not change up to ecor /ed = 0.3 for all values of α. For lower values of ecor /ed (0.2 and 0.1) the system still maintains the same type of ordering as in Fig. 7 at high values of α ≈ 2/3, but for lower values of α a strong competition of the C1 and C2 domains with the flower-type and HG phases starts to occur.

Coronene in pores of symmetric molecules The importance of the host-coronene interaction in a system of BHTC molecules raises the question of the effect this interaction might have on higher symmetry molecules, TMA and BTB. For these molecules the model (2) is much simpler and might be reduced to the 3-state Bell-Lavis model. 21,27 The honeycomb (HON) phase, which has the dimeric bonds only, is experimentally observed in both these systems (see, e.g. TMA, 7,29,30 BTB 38,39 and mixed TMA-BTB 40 ). In TMA system the HON phase is the first member (n = 1) of homologous series of flower phases. Further members of this series (n = 2, ..., ∞) emerge in TMA molecular system when the coverage of TMA molecules exceeds the coverage of the HON phase, chon tma = 1/6. Increase of ctma > chon tma is characterized by decrease in number of dimeric and increase in number of trimeric bonds between the molecules. The last member (n = ∞) of the series is called the superflower (SF) phase, and is characterized by trimeric bonds only. Alternatively, one can 22

ACS Paragon Plus Environment

Page 22 of 34

Page 23 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

say that with increase of molecular density the decrease in number of larger (honeycomb or six-molecule) pores and increase in number of smaller (triangular or three-molecule and rectangular or four-molecule) pores takes place. 7,22 The six-molecule pore (Fig. 8a) is the same-type of pore we encountered studying the ordering of the C1 phase of BHTC molecules, and this pore can incorporate a coronene molecule. The three- and four-molecule pores in the TMA flower phases are too small to guest a coronene molecule. Therefore, in calculations for TMA-coronene system with ecor > 0 one can expect the coronene in honeycomb (six-molecule) pores of the HON phase for any concentration of 36 coronene, if 2ccor ≤ chon In flower tma , and this is corroborated by the DFT calculations.

phases the coronene could also be incorporated in the honeycomb pores, but the ccor should be rather small, since any increase of molecular concentration, both TMA or coronene, might violate the molecular density (and a subtle balance between dimeric and trimeric bonds) at which the flower phases exist. The most interesting results for TMA-coronene system, performed by Kawasaki dynamics, are presented in Fig. 8. For optimal ratio of coronene/TMA density equal to 0.5 in a perfect HON phase, the six-molecule pores are fully occupied by the coronene molecules (Fig. 8b). For slightly higher TMA concentration (ctma = 2/9) either the flower (n = 2) phase or the HON and n = 3 two phase co-existence can be found without coronene. In our calculations here we obtain the two-phase co-existence. Introducing rather small amount of coronene interacting with TMA molecules (ecor /ed = 0.4) we obtain very interesting phase separation into three phases: HON, n = 3 and SF (Fig. 8c). Further increase of coronene concentration leads to co-existence of HON and SF phases (Fig. 8d). It should be noted, that this result agrees with the experimental fact that the SF phase might be found only as high density fragments in co-existence with other molecules (e. g. alcohols 41 ). The BTB molecules also can form the simplest HON phase. 15,38 In contrast to the sixmolecule honeycomb pore of the TMA molecules, such a pore of the BTB molecules is too large to appropriately host a coronene molecule, since the diameter of BTB molecule is

23

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 8: (a) Coronene molecule in a six-molecule pore of TMA HON phase. Snapshots of TMA-coronene system for different TMA coverages and different ratios of coronene coverage to TMA coverage: (b) ctma = Ntma /L2 = 1/9 and ccor /ctma = 0.5, (c) ctma = 2/9 and ccor /ctma = 0.05, (d) ctma = 2/9 and ccor /ctma = 0.17. Here L = 60, et /ed = 1.4, ecor /ed = 0.4 (b),(c) and 0.2 (d). roughly two times larger than that of the TMA (Fig. 1). The exact number of coronene guest molecules adsorbed in each pore of the BTB HON phase was not determined by the STM measurements, but it was predicted, 15 on a basis of molecular mechanics and molecular dynamics calculations and entropic arguments, that this pore can hold three coronene molecules. It is also known that pure BTB molecular system does not form any flower phases for higher BTB concentrations, because the low conformational flexibility of BTB is opposed to the tendency for dense packing. 38 Using our simulations, we tried to determine if the interaction with coronene might give extra stability to the BTB molecular system. We explored two cases: (i) when the coverage of BTB molecules corresponds to the stoichiometric coverage of the HON phase, cbtb ≈ chon btb , but the coronene coverage is even higher. Here we tried to find if the stability of the HON phase with three coronene molecules in a pore is plausible; (ii) when the coverage of BTB molecules exceeds the stoichiometric coverage of the HON phase, cbtb > chon btb . Here we studied if the coronene can induce extra flower-type phases at higher BTB concentrations. Like for BHTC molecules, we raised the question, if ecor might support the trio interactions which are not strong enough to independently form the trimeric bonds. The snapshot of our calculations for case (i) is presented in Fig. 9b. The prediction 15 24

ACS Paragon Plus Environment

Page 24 of 34

Page 25 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

of three molecules in a pore is corroborated also by our results - we find only one stable structure of coronene molecules in a pore of the BTB HON phase with three coronene molecules arranged symmetrically in an equilateral triangle (see Fig. 9a).

Figure 9: (a) Three coronene molecules in a six-molecule pore of BTB HON phase. Snapshots of BTB-coronene system for different BTB coverages and different ratios of coronene coverage to BTB coverage: (b) cbtb = Nbtb /L2 = 7/180 and ccor /cbtb = 1.5, (c) cbtb = 1/15 and ccor /cbtb = 0.58, (d) cbtb = 3/40 and ccor /cbtb = 1. Here L = 60, et /ed = 1.2 (b) and 1.4 (c),(d) and ecor /ed = 0.2. (e) Coronene molecule in a three-molecule pore of BTB molecules (superflower phase). While studying case (ii), we have found that the hexagonal pore organized by three BTB molecules (Fig. 9e) has an ideal size to host a coronene molecule. Following the analogy with the TMA orderings, we consider the BTB molecular system of such three-molecule pores as the SF phase of BTB molecules. Further simulations confirmed our prediction. For ideal ratio of coronene coverage to BTB coverage we obtain the SF structure with coronene in centers of three-molecule pores (Fig. 9d). For lower value of this ratio we obtain a very interesting two-phase co-existence of coronene-free HON structure and coronene-filled SF phase (Fig. 9c). It should be noted, that the coronene molecule in a six-molecule pore of BHTC-coronene or TMA-coronene system should be more stable than that in a three-molecule pore of BTBcoronene SF system. Actually, the reason for occurrence of the host-coronene interaction is the weakening of main dimeric double H-bond, since one carboxyl group starts to participate in two H-bonds (host-host dimeric and host-guest, ecor ). The negative charge which can be transferred from six H-bonds of six-molecule pore to coronene is twice less than might be given by three H-bonds in three-molecule BTB pore. This is seen comparing the schematic 25

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figs. 8a and 9a. We explored very small BTB-coronene interactions, and coronene-supported BTB SF phase still emerged in our simulations.

Conclusions Coronene molecule is an important element strengthening the ordering of tricarboxylic acid molecules. Owing to similarity of hexagonal lattices of coronene and graphite, an almost perfect π − π stacking is established between the coronene and substrate. As a result, the TMA HON structure with coronene located in honeycomb pores is almost commensurate with the substrate lattice, while the same structure without coronene demonstrates some incommensuration. 37 Thus, it is more likely that the coronene rather than the BHTC molecule is the main condensation center when the coronene is added to the BHTC system. Especially as it is known 36 that the coronene gives extra stability to the hexagonal pore of the TMA HON phase. 32 Thus the occurrence 28 of a coronene-filled six-molecule pore could be expected. However, in addition to dimerically bonded six-molecule pore, further ordering of the BHTC-coronene system requires the pores of similar size which also could contain coronene molecules. The trimeric interactions are inevitable to form such pores, since it is impossible to enclose the coronene by BHTC molecules bonded by the dimeric interactions only. In such a way, the trimeric interactions, which are not strong enough to independently organize the C1 structure of BHTC molecules, gain some aid from neighboring coroneneBHTC interaction. Along with the C1 structure, the MC calculations predict the C2 structure which also creates the hexagonal four-molecule pores with dimeric and trimeric interactions. In addition to those pores, the C2 structure has five-molecule pores and is characterized by alternating zigzag chains of four-molecule and five-molecule pores. The size of a pore is appropriate to hold a coronene molecule. In our calculations the mixing of both C1 and C2 phases is very often observed. The

26

ACS Paragon Plus Environment

Page 26 of 34

Page 27 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

mixing is possible, because the zigzag chain of four-molecule pores is a common building block of both C1 and C2 phases, and two short-short dimeric interactions, the only interactions which make the energies of both phases different, have very similar magnitudes. Moreover, here we do not have the two-phase co-existence, which is characteristic to two structures of different molecular density accompanied by the first-order phase transition, but a simple mixing of two phases of the same density. Nevertheless, our calculations for ess1 ≈ ess2 demonstrate a certain domination of the C2 phase over the C1 phase at full occupancy of pores by coronene. The question might arise if the existence of coronene molecules prompts the mixing of phases C1 and C2. The answer is, most likely, negative: the MC calculation for BHTC molecules for ess1 ≈ ess2 without addition of coronene, but with too large trio interaction instead, often results in mixed C1 and C2 structure with the domination of the latter phase. We relate this result to higher entropy of the C2 phase at finite temperature. However, at low coronene concentrations ccor the six-molecule pores characteristic to the C1 structure are more preferred, because larger number of neighboring BHTC molecules are interacting with coronene molecule (more dimeric BHTC bonds are being shared to establish the BHTC-coronene interaction). But with increase of ccor the dominating structure might change, and the share of each phase depends on a competition between values of the short-short dimeric interactions, difference in entropy of phases C1 and C2, and coronene concentration. Since the six-molecule pore in a symmetry-reduced BHTC system is the same as the pore created in a honeycomb phase of symmetric TMA molecules, we studied the effect of corone on ordering in TMA system. The model was reduced to the 3-state Bell-Lavis model used 21,27 for TMA molecules ordering. We found that coronene molecules easily fill the pores in a honeycomb phase. For higher TMA densities the coronene still fills the hexagonal pores, but its insertion affects the subtle balance of existing phases. Addition of a small number of coronene molecules leads to a phase co-existence between the honeycomb, third flower

27

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 34

and superflower phases. For slightly higher ccor , this co-existence transforms into a nicely defined co-existence between coronene-filled honeycomb phase and coronene-free superflower structure. The pores in a honeycomb phase of BTB molecular system are much larger than those found in TMA and any of BHTC phases. Analyzing the sizes of pores which might hold the coronene molecule, we detected the hexagonal three-molecule pore which might be formed in a a structure with trimeric interactions only. These pores of the BTB superflower phase are perfect to hold a coronene, though the host-coronene interaction should be lower than in a TMA-corone six-molecule pore. On the other hand, the BTB system is known as not forming dense phases due to low conformational flexibility of BTB molecules. Thus, our MC calculations was supposed to demonstrate if the coronene-induced superflower phase can occur in a BTB molecular system when BTB-coronene interaction is rather small. Our results predict the formation of the coronene-filled superflower phase. At lower concentrations of coronene they demonstrate the co-existence of a coronene-filled superflower phase with a coronene-free honeycomb structure. At even lower BTB concentration, corresponding to that of the honeycomb phase, and high corone concentration, our calculations demonstrate that three coronene molecules are stable in a hexagonal pore of BTB molecules.

References (1) Bartels, L. Tailoring Molecular Layers at Metal Surfaces. Nature Chemistry 2010, 2, 87-95. (2) Barth, J. V. Molecular Architectonic on Metal Surfaces. Annu. Rev. Phys. Chem. 2007, 58, 375-407. (3) Wang, D.; Wan, L. J.; Bai, C. L. Formation and Structural Transition of Molecular Self-Assembly on Solid Surface Investigated by Scanning Tunneling Microscopy. Materials Science and Engineering R 2010, 70, 169-187.

28

ACS Paragon Plus Environment

Page 29 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(4) K¨ uhnle, A. Self-Assembly of Organic Molecules at Metal Surfaces. Current Opinion in Colloid and Interface Science 2009, 14, 157-168. (5) Rahe, P.; Nimmrich, M.; K¨ uhnle, A. Substrate Templating upon Self-Assembly of Hydrogen Bonded Molecular Networks on an Insulating Surface. Small 2012, 8, 29692977. (6) Ha, N. T. N.; Gopakumar, T. G.; Hietschold, M. Polymorphs of Trimesic Acid Controlled by Solvent Polarity and Concentration of Solute at Solid-Liquid Interface. Surface Science 2013, 607, 68-73. (7) Ye, Y. C.; Sun, W.; Wang, Y. F.; Shao, X.; Xu, X. G.; Cheng, F.; Li, J. L.; Wu, K. A Unified Model: Self-Assembly of Trimesic Acid on Gold. J. Phys. Chem. C 2007, 111, 10138-10141. (8) Steiner, T. The Hydrogen Bond in the Solid State. Angew. Chem. Int. Ed. Engl. 2002, 41, 48-76. (9) Mura, M.; Gulans, A.; Thonhauser, T.; Kantorovich, L. Role of van der Waals Interaction in Forming Molecule-Metal Junctions: Flat Organic Molecules on the Au(111) Surface. Phys. Chem. Chem. Phys. 2010, 12, 4759-4767. (10) Mura, M.; Martsinovich, N.; Kantorovich, L. Theoretical Study of Melamine Superstructures and Their Interaction With the Au(111) Surface. Nanotechnology 2008, 19, 465704. (11) Mura, M.; Silly, F.; Burlakov, V.; Castell, M. R.; Briggs, G. A. D.; Kantorovich, L. N. Formation Mechanism for a Hybrid Supramolecular Network Involving Cooperative Interactions. Phys. Rev. Lett. 2012, 108, 176103. (12) Blunt, M. O.; Adisoejoso, J.; Tahara, K.; Katayama, K.; Van der Auweraer, M.;

29

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Tobe, Y.; De Feyter, S. Temperature-Induced Structural Phase Transitions in a TwoDimensional Self-Assembled Network. J. Am. Chem. Soc. 2013, 135, 12068-12075. (13) Ha, N. T. N.; Gopakumar, T. G.; Gutzler, R.; Lackinger, M.; Tang H.; Hietschold, M. Influence of Solvophobic Effects on Self-Assembly of Trimesic Acid at the Liquid-Solid Interface. J. Phys. Chem. C 2010, 114, 3531-3536. (14) Martsinovich, N.; Troisi, A. Modeling the Self-Assembly of Benzenedicarboxylic Acids Using Monte Carlo and Molecular Dynamics Simulations. J. Phys. Chem. C 2010, 114, 4376-4388. (15) Eder, G.; Kloft, S.; Martsinovich, N.; Mahata, K.; Schmittel, M.; Heckl, W. M.; Lackinger, M. Incorporation Dynamics of Molecular Guests into Two-Dimensional Supramolecular Host Networks at the Liquid-Solid Interface. Langmuir 2011, 27, 13563-13571. (16) Weber, U. K.; Burlakov, V. M.; Perdigao, L. M. A.; Fawcett, R. H. J.; Beton, P. H.; Champness, N. R.; Jefferson, J. H.; Briggs, G. A. D.; Pettifor, D. G. Role of Interaction Anisotropy in the Formation and Stability of Molecular Templates. Phys. Rev. Lett. 2008, 100, 156101. (17) Silly, F.; Weber, U. K.; Shaw, A. Q.; Burlakov, V. M.; Castell, M. R.; Briggs, G. A. D.; Pettifor, D. G. Deriving Molecular Bonding from a Macromolecular Self-Assembly Using Kinetic Monte Carlo Simulations. Phys. Rev. B 2008, 77, 201408. (18) Balb´as Gambra, M.; Rohr, C.; Gruber, K.; Hermann, B. A.; Franosch, T. Simulating Self-Organized Molecular Patterns Using Interaction-Site Models. Eur. Phys. J. E 2012, 35, 25. (19) Szabelski, P.; De Feyter, S. Chiral Occlusion in Two-Dimensional Binary Supramolecular Networks Studied by the Monte Carlo Method. CrystEngComm 2011, 13, 55425550. 30

ACS Paragon Plus Environment

Page 30 of 34

Page 31 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(20) Kasperski, A.; Szabelski, P. Theoretical Modeling of the Formation of Chiral Molecular Patterns in Self-Assembled Overlayers. Surface Science 2014, 629, 57-64. (21) Misi¯ unas, T.; Tornau, E. E. Ordered Assemblies of Triangular-Shaped Molecules with Strongly Interacting Vertices: Phase Diagrams for Honeycomb and Zigzag Structures on Triangular Lattice. J. Phys. Chem. B 2012, 116, 2472-2482. (22) Ibenskas A.; Tornau, E. E. Statistical Model for Self-Assembly of Trimesic Acid Molecules into Homologous Series of Flower Phases. Phys. Rev. E 2012, 86, 051118. ˇ enas, M.; Tornau, E. E. Pin-Wheel Hexagons: A Model for Anthraquinone Order(23) Sim˙ ing on Cu(111). J. Chem. Phys. 2013, 139, 154711. ˇ enas, M.; Tornau, E. E. A Model of Melamine Molecules Ordering on Metal Sur(24) Sim˙ faces. J. Chem. Phys. 2014, 141, 054701. (25) Bell, G. M.; Lavis, D. A. Two-Dimensional Bonded Lattice Fluids. II. Orientable Molecule Model. J. Phys. A: Gen. Phys. 1970, 3, 568-581. (26) Barbosa, M. A. A.; Henriques, V. B. Frustration and Anomalous Behavior in the Bell-Lavis Model of Liquid Water. Phys. Rev. E 2008, 77, 051204. ˇ enas, M.; Ibenskas, A.; Tornau, E. E. Phase Transition Properties of the Bell-Lavis (27) Sim˙ Model. Phys. Rev. E 2014, 90, 042124. (28) MacLeod, J. M.; Chaouch, Z. B.; Perepichka, D. F.; Rosei, F. Two-Dimensional SelfAssembly of a Symmetry-Reduced Tricarboxylic Acid. Langmuir, 2013, 29, 7318-7324. (29) Griessl, S.; Lackinger, M.; Edelwirth, M.; Hietschold, M.; Heckl, W. M. Self-Assembled Two-Dimensional Molecular Host-Guest Architectures From Trimesic Acid. Single Mol. 2002, 3, 25-31.

31

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(30) Lackinger, M.; Griessl, S.; Heckl, W. M.; Hietschold, M.; Flynn, G. W. Self-Assembly of Trimesic Acid at the Liquid-Solid Interfacea Study of Solvent-Induced Polymorphism. Langmuir 2005, 31, 4984-4988. (31) Nath, K. G.; Ivasenko, O.; MacLeod, J. M.; Miwa, J. A.; Wuest, J. D.; Nanci, A.; Perepichka, D. F.; Rosei, F. Crystal Engineering in Two Dimensions: An Approach to Molecular Nanopatterning. J. Phys. Chem. C 2007, 111, 16996-17007. (32) Blunt, M. O.; Russell, J. C.; del Carmen Gim´enez-L´opez, M.; Garrahan, J. P.; Lin, X.; Schr¨oder, M.; Champness, N. R.; Beton, P. H. Random Tiling and Topological Defects in a Two-Dimensional Molecular Network. Science 2008, 322, 1077-1081. (33) Landau, D. P.; Binder, K. A guide to Monte Carlo simulations in statistical physics, Cambridge Uni. Press, 3rd. edition, 2009. (34) Wang, Y.; Lingenfelder, M.; Fabris, S.; Fratesi, G.; Ferrando, R.; Classen, T.; Kern, K.; Costantini, G. Programming Hierarchical Supramolecular Nanostructures by Molecular Design. J. Phys. Chem. C 2013, 117, 3440-3445. (35) Griessl, S.; Lackinger, M.; Jamitzky, F.; Markert, T.; Hietschold, M.; Heckl, W. M. Incorporation and Manipulation of Coronene in an Organic Template Structure. Langmuir 2004, 20, 9403-9407. (36) Blunt, M.; Lin, X.; del Carmen Gim´enez-L´opez, M.; Schr¨oder, M.; Champness, N. R.; Beton, P. H. Directing Two-Dimensional Molecular Crystallization Using Guest Templates. Chem. Commun. 2008, 2304-2306. (37) Li, M.; Deng, K.; Yang, Y. L; Zeng, Q. D; He, M.; Wang, C. Electronically Engineered Interface Molecular Superlattices: STM Study of Aromatic Molecules on Graphite. Phys. Rev. B 2007, 76, 155438.

32

ACS Paragon Plus Environment

Page 32 of 34

Page 33 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(38) Kampschulte, L.; Lackinger, M.; Maier, A. K.; Kishore, R. S. K.; Griessl, S.; Schmittel, M.; Heckl, W. M. Solvent Induced Polymorphism in Supramolecular 1,3,5Benzenetribenzoic Acid Monolayers. J. Phys. Chem. B 2006, 110, 10829-10836. (39) Gutzler, R.; Sirtl, T.; Dienstmaier, J. F.; Mahata, K.; Heckl, W. M.; Schmittel M.; Lackinger, M. Reversible Phase Transitions in Self-Assembled Monolayers at the Liquid-Solid Interface: Temperature-Controlled Opening and Closing of Nanopores. J. Am. Chem. Soc. 2010, 132, 5084-5090. (40) Kampschulte, L.; Werblowsky, T. L.; Kishore, R. S. K.; Schmittel, M.; Heckl, W. M.; Lackinger, M. Thermodynamical Equilibrium of Binary Supramolecular Networks at the LiquidSolid Interface. J. Am. Chem. Soc. 2008, 130, 8502-8507. (41) MacLeod, J. M.; Ivasenko, O.; Perepichka, D. F.; Rosei, F. Stabilization of Exotic Minority Phases in a Multicomponent Self-Assembled Molecular Network. Nanotechnology 2007, 18, 424031.

33

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

34

ACS Paragon Plus Environment

Page 34 of 34