Reliable Computational Prediction of the Supramolecular Ordering of Complex Molecules under Electrochemical Conditions.

We propose a computationally lean, two-stage approach that reliably predicts self-assembly behavior of complex charged molecules on metallic surfaces under electrochemical conditions. Stage one uses ab initio simulations to provide reference data for the energies (evaluated for archetypical configurations) to fit the parameters of a conceptually much simpler and computationally less expensive force field of the molecules: classical, spherical particles, representing the respective atomic entities; a flat and perfectly conducting wall represents the metallic surface. Stage two feeds the energies that emerge from this force field into highly efficient and reliable optimization techniques to identify via energy minimization the ordered ground-state configurations of the molecules. We demonstrate the power of our approach by successfully reproducing, on a semiquantitative level, the intricate supramolecular ordering observed experimentally for PQP+ and ClO4- molecules at an Au(111)-electrolyte interface, including the formation of open-porous, self-host-guest, and stratified bilayer phases as a function of the electric field at the solid-liquid interface. We also discuss the role of the perchlorate ions in the self-assembly process, whose positions could not be identified in the related experimental investigations.

Introduction

Supramolecular chemistry deals with intermolecular interactions and structure formation beyond individual molecules and as such lies at the base of many nano- and mesoscopic structures found in biology. In recent decades, impressive progress in the experimental branches of this field has resulted in at least two Nobel Prizes in chemistry. By contrast, the theoretical understanding and especially the in silico prediction of supramolecular ordering has lagged behind somewhat. This is easily understood if one considers the sheer size of the systems under study, requiring in many cases consideration of a solid substrate, a sufficiently large number of molecular building blocks or tectons, and a condensed matter medium (i.e., a solvent or electrolyte solution). The interaction of these three components, each with its intrinsic properties and with optional extrinsic steering (e.g., by light, heat, electric field), will determine the observed supramolecular structures and govern the transitions between them.[1,2]

In this paper, we propose a new theoretical framework to predict the supramolecular ordering of complex molecules at an electrochemical solid–liquid interface. The calculations were inspired by recent experimental work[3] in which particularly clear-cut transitions between supramolecular structures were observed as a function of the applied electric field at a metal–electrolyte interface. The target molecules whose supramolecular ordering is considered constitute an organic salt that consists of a large, disc-shaped polyaromatic cation (PQP+) and a much smaller, inorganic anion (perchlorate, ClO4–).[4-600]

The concept of choice to investigate these scenarios would rely on (i) a faithful description of the properties of the system (notably a reliable evaluation of its energy) via ab initio simulations and (ii) in a subsequent step the identification of the optimized (ordered) arrangement of the molecules on the substrate by minimizing this energy via efficient and reliable numerical tools; this optimization has to be performed in a high-dimensional search space, spanning all possible geometries of the unit cell and all possible coordinates and orientations of the molecules within that cell. Both these approaches, considered separately from each other, are conceptually highly complex and from the numerical point of view very expensive, which precludes the application of this combined concept even for a single set of external parameters (such as temperature, density, and external field); it is thus obvious that systematic investigations of the self-assembly scenarios of such systems are definitely out of reach.

In this contribution, we propose an approach to overcome these limitations via the following strategy: in the first step, we map the ab initio-based energies onto the energy of a related classical model (or a classical force field), where the atomistic units of the molecules are featured as spherical, charged units with Lennard-Jones-type interactions and where the electrolyte is treated as a homogeneous, dielectric medium; the interaction between the atomic entities and the metallic surface is modeled by a classical, perfectly conductive, Lennard-Jones-like wall potential. The as-yet open parameters of the resulting force field (energy- and length-scales, charges, etc.) are fixed by matching the ab initio energies of the system with the related energies of this force field: this is achieved by considering archetypical configurations of the system’s building blocks (molecules and surface) and by systematically varying characteristic distances between these units over a representative range. These ab initio energies were then fitted along these “trajectories” by the parameters of the classical force field: the energy- and length-scales of the involved interatomic Lennard-Jones or Mie potentials as well as the atom–wall interaction parameters.

It turns out that this force field is indeed able to reproduce the ab initio-based energies along these “trajectories” faithfully and with high accuracy. Even though the emerging force field is still quite complex (as it features both short-range and long-range Coulomb interactions and involves mirror charges), it is now amenable to the aforementioned optimization techniques, which thus brings systematic investigations of the self-assembly scenarios of these molecules under the variations of external parameters within reach.

As a benchmark test for our approach, we have considered the above-mentioned system, studied in recent experimental investigations: the cation is PQP+ (9-phenylbenzo[1,2] quinolizino[3,4,5,6-fed]phenanthridinylium, a disk-shaped polyaromatic molecule), while the anion is perchlorate, ClO4–; the self-assembly of these ions on a Au(111) surface under the influence of an external electric field was studied. The high accuracy with which the ensuing energies calculated from the force field reproduce the ab initio simulation data make us confident about the applicability of the force field for the subsequent optimization step: using an optimization technique that is based on ideas of evolutionary algorithms, we have then identified the self-assembly scenarios of the ions on the Au surface, for a given set of external parameters (temperature, density, and external field). These first results provide evidence that our approach is quite promising. This concept is furthermore completely flexible as it can easily be extended to other organic molecules of similar (or even higher) complexity. The computational cost of this optimization step is still substantial. Therefore, we postpone a detailed, quantitative, and, in particular, systematic investigation of the self-assembly scenarios of the PQP+ and the ClO4– ions on the Au surface for a broad range of external parameters to a later contribution. Instead, we demonstrate in this contribution for selected sets of parameters that our approach is indeed able to reproduce several of the experimentally identified self-assembly scenarios.

In this context, it has to be emphasized that such a type of optimization problem is highly nontrivial since the huge number of possible local minima in the potential energy surface (embedded in a high-dimensional parameter space) increases exponentially with the number of particles (and their degrees of freedom) of the system;[6] thus, exhaustive search strategies hit the computational limits or even exceed the capacities of present-day supercomputers. Yet, another complication in structure prediction algorithms is caused by the fact that different polymorphs of a system can be kinetically trapped and a vast number of other minima, having values of the internal energy comparable to the global minimum, may also play an important role in structure formation processes.[6,7]

At this point, we owe an explanation to the reader why we have chosen the possibly unconventional approach. Of course, it is obvious that an optimization of the molecular configurations on the basis of full ab initio calculations is from the computational point of view by far out of reach. However, one can argue that suitable force fields (available in the literature) or machine-learning (ML) potentials such as high-dimensional neural network potentials,[8−14] kernel-based ML methods[15] (such as Gaussian approximation potentials[16−20]), or more specialized, effective potentials for selected molecular motives such as the SAMPLE approach[7,21,22] might represent a more conventional approach to tackle this problem (note that the field of ML potentials is rapidly growing and the above list is far from comprehensive). Such arguments represent fully legitimate objections against our approach.

The problem we are addressing in this contribution is, however, a nonstandard problem and thus requires to be treated with a custom force field: the justification for our strategy is that we wanted to endow the atomic units of the molecules with “real” physical properties (such as “size” or “charge”), which will help us to perform the second step in our structural search that we have envisaged (and that we are currently working on): as the computational costs of our approach are still considerable, large-scale investigations are still prohibitively expensive. In an effort to overcome these limitations, we plan to proceed to even more coarse-grained models, which grasp, nevertheless, the essential features of our complex molecules. On the basis of such models, we would then be able to identify with rather low computational costs first trends in structural identification processes. Investigations along this direction will be explored in forthcoming contributions.

Finally, we point out that we are well aware of the limitations and deficiencies of our present model. Features such as the response of the metallic electronic distribution of the gold surface due to the presence of an external bias, variable electrostatic properties of the molecular species (allowing thus for polarization effects), or a space-dependent permittivity cannot be included in our concept. However, at this point, it is fair to say that, to the best of our knowledge, none of the aforementioned alternative approaches (such as the use of conventional force fields or machine-learning frameworks) is able to take all of these effects faithfully into account, either.

The manuscript is organized as follows: In Section , we describe the essential features of the experimental setup, introduce an ab initio and a classical representation thereof, and discuss the mapping procedure between those different instances. In Section , we put forward the memetic optimization procedure based on ideas of evolutionary algorithms to identify ordered ground-state configurations of complex molecules under electrochemical conditions, and in Section , we present selected numerical results, which demonstrate a semiquantitative agreement with the experimentally observed self-assembly scenarios of PQP+ and ClO4– ions on a Au(111)–electrolyte interface under the influence of an external electrostatic field. We conclude our findings in Section .

System and Its Representations

System

Both the density functional theory (DFT) calculations and the related force field are based on a framework that mimics the essential features of the experimental setup, put forward (and discussed) in ref (3); this framework is schematically depicted in Figure : PQP+ and ClO4– ions are immersed in an electrolyte (0.1 M aqueous perchloric acid). From below, the system is confined by a Au(111) surface, which in the experiment serves as the solid substrate for adsorption. An electric field, E, can be applied between a reference electrode located within the electrolyte and the Au surface. The PQP+ and the ClO4– ions are first treated via DFT-based ab initio calculations (see Section ). The calculated energies are then used to fix the force fields of classical particles (notably their sizes, energy parameters, and charges), which represent the atomic entities of the respective ions; the interaction between the atomic entities and the Au(111) substrate is described by means of a classical wall-particle force field (see Section ). Throughout, the electrolyte molecules are not considered explicitly. The electrolyte is rather assumed to be an effective homogeneous medium with a permittivity of water, i.e., ϵr = 78.36, at T = 25 °C,[23−25] corresponding to the temperature at which the experiments by Cui et al.[3] were carried out and assuming that the low concentration of perchloric acid does not change the value of ϵr substantially.[26−29] Hence, in this contribution, we use “electrolyte” as a synonym for a “solvent” unless explicit use is required.

Figure 1

Schematic visualization of the experimental setup to control the pattern formation of PQP+ (and ClO4–) molecules (structure formulas given in top-right insets) close to a Au(111) surface: two Au layers are explicitly shown; the golden, shiny area represents the conductive Au bulk, and the black dashed line marks the surface of the electronic density, which we interpret as a mirror plane. The ions are immersed into an electrolyte (gray, shaded region), which is considered as an effective, homogeneous medium. In the region close to the Au surface (red to blue shaded areas), a homogeneous, electrostatic field E (bold, colored arrow), oriented in the z-direction, features the electrostatic potential drop between the Au surface and the reference electrode inside the electrolyte. The colors of the atoms in the electrolyte correspond to their type, while the color of the mirror atoms (located in the Au bulk) specifies their partial charges, quantified by the color bar (see bottom right) in units of the electron charge, e.

We emphasize at this point that in the experiment, an exact specification of the electric field strength is not possible: as detailed in the supporting information of ref (3), the authors of the related experimental investigations have estimated rather the degree of charge compensation on the Au surface by the adsorbed PQP+ ions as a function of their changing coverage, which does not allow the estimation of the electric field directly. This fact limits the degree of quantitative comparison between experiment and theory.

Ab Initio Simulations

The density functional theory calculations were performed with the software package GPAW,[30,31] and the structures were handled by the atomic simulation environment.[32] The electronic density and the Kohn–Sham orbitals were represented within the projector-augmented wave method,[33] where the smooth parts were represented on real space grids with a grid spacing of 0.2 Å for the orbitals and 0.1 Å for the electron density. The exchange-correlation energy is approximated as proposed by Perdew, Burke, and Ernzerhof (PBE),[34] and weak interactions missing in the PBE functional are described as proposed by Tkatchenko and Scheffler (TS09).[35] The TS09 approximation assumes that long-range dispersive contributions are absent in the PBE functional, such that these can be applied as a correction. The total energy is written aswhere EPBE is the PBE energy and EvdW is the TS09 correction. We have introduced a weight factor wS that will allow the incorporation of electrolyte effects into the dispersive contributions as discussed below. For interactions in vacuum, wS = 1. The presence of the aqueous environment on the electronic and nuclear degrees of freedom included in EPBE is modeled by a continuum solvent model.[36]

Molecular interactions were studied on simulation grids with Dirichlet (zero) boundary conditions. Neumann (periodic) boundary conditions were applied in x- and y-directions in the surface plane for simulations involving the gold surface, while zero boundary conditions were applied in the perpendicular z-direction. The simulation grid was chosen such that at least 4 Å of space around the position of each atom in the nonperiodic directions was ensured. The Au(111) gold substrate was modeled by two layers of 54 atoms, each using the experimental lattice constant of face-centered cubic (fcc) gold of a = 4.08 Å. These settings result in a rectangular unit cell of 26.0 × 15.0 Å2. The Brillouin zone was sampled by 3 × 3 Monkhorst–Pack[37] distributed k-points in the periodic directions.

Potentials were scanned by fixing all gold atoms and a central atom of PQP+ (the nitrogen atom) and/or of ClO4– (the chlorine atom) to given positions, while all other atoms were allowed to relax without any symmetry constraints until all forces were below 0.05 eV/Å.

The interaction of two perchlorate anions in dependence of their distance is shown in Figure for different approximations for the total energy in eq . As expected, the potentials follow the screened electrostatic repulsion [ε∞(0)RCl–Cl]−1 for large distances RCl–Cl, where ε∞(0) = εr is the static relative permittivity of water. There is a slight attractive part in the potential around RCl–Cl ≃ 5.2 Å already in the PBE potential, which leads to a very shallow local minimum. The main reason for this minimum is the decrease in the effective surface ΔA when the solute cavities (the solvent-excluded regions) begin to overlap. This decreases the energetic cost to form the surface due to the effective surface tension γ = 18.4 dyn cm–1 (γ also contains attractive contributions and is therefore much lower than the experimental surface tension of water[36]). Subtracting γΔA nearly removes all of the minima as demonstrated in Figure . By including the full dispersion contribution [ε∞(ωopt) = 1], this local minimum substantially deepens and becomes the total minimum of the potential. An attractive contribution to the potential is not to be expected for the interaction of two anions and needs further discussion. We suspect an overestimation of dispersion interactions if these are treated as in vacuum and no screening through the electrolyte is considered.

Figure 2

Relative energy of two ClO4– anions as a function of the distance between their chlorine atoms, RCl–Cl, where the separated anions define the energy reference; ε(ωopt) = 1 with full van der Waals (vdW) corrections and ε(ωopt) = 1.7 with scaled vdW corrections. The dash-dotted line shows the PBE energy, where the energy contribution of effective surface tension γΔA is subtracted (see the text).

The aqueous environment influences the van der Waals (vdW) interactions as these are of Coulombic origin.[38] To derive an approximate expression for the screened vdW interaction of two molecules A and B at distance R inside the electrolyte, we express the C6 coefficient defining the vdW energy C6/R6 by the Casimir–Polder integral[35,39]where αA,B* is the polarizability of the interacting molecules and ϕ is determined by the propagation of the electric field through the embedding medium[39] with ϕ = 1 in vacuum. Both αA,B* and ϕ are modified relative to vacuum in solution. In the simplest model,[39] we may write ϕ(iξ) = ε∞–2(iξ) with the frequency-dependent relative permittivity of the electrolyte ε∞. The effects of the electrolyte on the polarizabilities αA,B* should, at least partly, already be included in the TS09 description through the effective atomic polarizabilities derived from the self-consistent electron density calculated within the electrolyte. What is left is the effect of the permittivity entering through the function ϕ(iξ) in eq . We assume that the main contribution of ϕ(iξ) is at the resonance frequencies of αA,B*, which are in the optical region for usual molecules. We further assume that ε∞(ωopt) is approximately constant in this frequency region such that we may pull ϕ = [ε∞(ωopt)]−2 out of the integral. This factor scales the C6 coefficient and therefore the vdW contribution. In other words, we apply the weight wS = [ε∞(ωopt)]−2 in eq with the experimental permittivity of water in the optical region of ϵ∞(ωopt) = 1.7; see ref (40). This approach reduces considerably the depth of the suspiciously deep minimum as seen in Figure such that only a shallow local minimum remains similar to the PBE potential. The reduction obtained is quite strong in respect of the small contributions of Axilrod–Teller–Muto interactions commonly assumed.[41,42] The quantitative connection between the screening of dispersive interactions in polarizable media and the many-body effects neglected in TS09[43] is not immediately clear and is certainly worth further investigation. In what follows, we use the same scaling for all of the vdW contributions of the DFT potentials in this work.

Force Field Model

In this section, we describe how we cast our setup into force fields where the atoms in the molecular constituents are described as spherical particles, each of them carrying a charge. The mapping is guided by the energies obtained via the ab initio simulations detailed above.

The Au(111) surface is modeled as a flat and perfectly conducting surface involving mirror charges as detailed below. However, we note that the position of the corresponding surface in the DFT calculations does not coincide with the position of the atoms. Before proceeding, the following comment is in order: in this mapping procedure, the distance of a point charge to a metallic surface is unambiguously defined through the electrons leaking out of the potential defined by the nuclei.[44,45] This feature can explicitly be seen in jellium models[46] but also emerges in implicit calculations[47] where electrons spill out of the surface of metal clusters.[48] From the latter study, we estimate an effective spill out of the surface of 0.5 Å, a value that agrees qualitatively with estimates from the jellium models, extrapolated to large structures.[45] This value will be used below for our problem.

Atomistic Model

In our atomistic model, the molecules are represented as rigid entities composed of atomistic constituents. The molecules are immersed in a microscopic electrolyte, which is treated as a continuous medium of given permittivity. From below, the system is confined by a conducting Au(111) surface (which is assumed to extend in the x- and y-directions), and an external field (with respect to the electrolyte) can be applied in the z-direction, i.e., perpendicular to the surface (or wall). Figure schematically depicts all details of this atomistic model for the PQP+ClO4– system, confined by the Au surface.

To specify the different entities of the system and their force fields, we use the following notation:

Each of a total number of N molecules is uniquely labeled by capital Latin indices I: for each of these units, this index is assigned to its center-of-mass (COM) position vector, R; to a vector P, specifying its orientation within the lab frame in terms of the angle-axis framework[49,50] (see Supporting Information (S.I.), Section S2.2 for more details), and to the set of coordinates, r, of the respective N atomistic constituents of the molecule in its COM frame (to which we also refer as its blueprint). The set of COM positions and orientation vectors of all N molecules are denoted R and P, respectively. The set of all atom positions in the lab frame is given by r, and the position of each atom in the lab frame is uniquely defined by a vector r, labeled with Latin indices (i = 1, ..., n).

Between all atoms, we consider long-range Coulombic interactions (index “C”)with the interatomic distance r = |r – r| and charges q and q of the units i and j, respectively; the dielectric constant ϵ0 and the relative permittivity ϵr specify the implicit electrolyte. Further, we introduce short-range force fields (index “sr”) for which we have considered two options: first, a Lennard-Jones potential (index “LJ”), i.e.for the energy and length parameters, ϵ and σ, we have opted for the standard Lorentz–Berthelot mixing rules,[51] i.e., σ = (1/2)(σ + σ) and , respectively. Alternatively, we have also considered for the short-range interactions the Mie potential[52] (index “Mie”), which can be considered as a generalization of the LJ interaction; its functional form is given byand allows for a variation of the exponents of the repulsive and attractive contributions to the potential, γ( and γ(, respectively. ϵ and σ are again parameters for the energy- and the length-scales. The C is defined as a function of the exponents[52]For the exponents, we apply arithmetic mixing laws, i.e., γ( = (1/2)(γ( + γ() and γ( = (1/2)(γ( +γ().

We assume the Au surface to be perfectly conductive. Consequently, we need to explicitly consider mirror charges in our model; when further assuming z = 0 as the plane of reflection, the Coulombic interaction becomeswith the mirror charges q = −q and their positions r = (x, y, z) = (x, y, −z).

We describe the solid–liquid interface in terms of a slab geometry with a lower confining wall, i.e., we assume periodicity in the x- and y-directions, but a finite extent, c, of the geometry in the z-direction, which is chosen such that no restriction in the orientation of any molecule occurs; thus, c ≈ 1.2–2 nm, given their size and the slab width. We define the (orthorhombic) lattice vectors, a = (a, 0, 0), b = (b, b, 0), and c = (0, 0, c), which, without the loss of generality, define the volume of the unit cell, V = abc, and which we collect within the matrix . Together with the molecular basis, given by R, P, and all N (rigid) molecular blueprints, r, we now define the supramolecular latticewhich gives rise to all atomic coordinates in the lab frame, r, i.e., the molecular crystal structure of the system (see S.I. Section S2.2).

The force field between the atomic entities and the Au surface is described via an LJ-type wall potential[53]In the above relation, z is the height of atom i above the surface and σw and ϵw are the length and energy parameters of the interactions of each atom i with the wall, respectively.

Finally, we express the electrostatic interfacial potential between the electrode and the Au surface by an external, homogeneous electrostatic field, E (i.e., oriented perpendicular to the surface): we account for this potential via U(field)(z) = zqE.[54]

Thus, and eventually, the total potential energy of our model is given by the expressionwith “sr” standing for “LJ” or “Mie”; we recall that r is the set of all n atomic positions r in a lattice with slab geometry (defined by the unit cell ). If not present (and not explicitly addressed), the electric field will be dropped in the argument list of eq , that is, . The notation “∑*” indicates that summation is only carried out over atoms, labeled with Latin indices i and j, which belong to different molecules I and J (with I ≠ J), i.e., molecules being labeled with capital indices. The energy given in eq and the corresponding force fields are efficiently evaluated using the software package LAMMPS.[55]

To evaluate the long-range Coulomb term, , in the given slab geometry, we use numerically reliable and efficient slab-corrected three-dimensional (3D) Ewald-summation techniques.[56−58] The other terms in eq are evaluated via direct lattice summation techniques.

Parametrizing the Classical Model via Ab Initio Calculations

In this work, the blueprint of each molecule r is obtained from electronic structure calculations based on density functional theory (DFT), using dispersion-corrected ab initio structure optimization,[34,35] as described in Section . The partial charges of the atoms, q, are parameterized via a Bader analysis[59] and are collected in Tables S2 and S3 in the S.I. Section S2.3. These charges are directly transferred to the atomic entities. We repeat that throughout the electrolyte molecules have not been considered explicitly: instead, we treat within the force field the electrolyte as an effective, homogeneous medium, introducing the electric permittivity of water ϵr.

To fix the remaining model parameters that specify the interactions in eq , we search for each atomistic entity (labeled i) the set of atomistic model parameters (specified below), which reproduces via eq the ab initio energies as good as possible. On the one side, we consider either the length and the energy parameters of the LJ potential (denoted ) or the length and the energy parameters together with the exponents of the Mie potentials (denoted ), as well as the wall parameters, . To fix these parameters, we proceed as follows:

We first perform ab initio structure optimization for different characteristic molecular configurations specified below. Here, molecules are either positioned next to each other (without considering the wall) or above the Au surface: in the former case, we fix the positions of two selected atoms belonging to different molecules, with the atoms being separated by r; in the latter case, we keep the height, z, of one selected atom above the surface constant. Relaxation of all other degrees of freedom leads in the ab initio simulations to spatially and orientationally optimized molecular structures; they are denoted d(r) and d(z), respectively, with corresponding energies and ; they are, themselves, functions of the interatomic distance, r, and the atom–wall separation, z, of the selected atoms.

For every optimized ab initio structure, d(r) and d(z), obtained in this manner, we define a corresponding molecular configuration r(r) and r(z), which is based on the above-introduced rigid molecular blueprint model r (with the index I running now over all N molecules present in the respective DFT structure). To this end, we synchronize the COM positions of each molecule I in the ab initio simulation with the corresponding COM positions R of its classical counterparts and align their orientation P accordingly.

Finally, we evaluate the corresponding energies with the help of the force field via eq at zero electric field, i.e., and . We search for the best set of parameters (or ) and via simultaneously minimizingOf course, in the model, the same unit cell, , and the same number of particles, n, as in the respective ab initio simulations have to be used. Note that in eq the wall term included in eq is obsolete since the surface atoms are not considered.

These fits are based on five particularly chosen, archetypical configurations, to be discussed below. In the panels of Figure , we display schematic sketches of these configurations of the PQP+ and ClO4– molecules; these panels show the corresponding energy curves obtained from the force field, with parameters based on a fitting procedure to the ab initio energy profiles.

Figure 3

Energies as obtained in ab initio simulations (black crosses) and fitted data, using the force field (involving LJ interactions (open blue circles) or Mie interactions (open orange diamonds)); see Section . Also shown are, with labels (a)–(e), five schematic sketches of the five archetypical configurations of the molecules (along with their relative displacements, schematically indicated via the arrows as the distances vary along the abscissa); the related energy curves are used to fit the parameters of the force field, as outlined in the text; the labels correspond to the itemization (a)–(e) used in Section . (f) PQP+ and the ClO4– molecules, drawn to scale and using the Mie force field for the short-range interactions: atomic entities are shown as transparent spheres with their diameters fixed by their respective optimized σ values and their Bader charges (see the color code).

Tail-to-tail configuration (see the inset in Figure a): We have considered a series of ab initio structure optimizations at constant but successively increasing nitrogen–nitrogen distances, rNN, in the x-direction (while keeping yNN and zNN constant) of an antiparallel oriented pair of PQP+ molecules; both cations are vertically decorated with a ClO4– molecule. The aromatic parts of the PQP+ molecules lie flat in the x- and y-directions such that their tails face each other.

Face-to-face configuration (see the inset in Figure b): In this case, we consider antiparallel oriented but vertically stacked PQP+ molecules (both being horizontally decorated by ClO4– molecules) under the variation of the nitrogen–nitrogen distance, rNN, in the z-direction (while now keeping xNN and yNN constant). Again, the aromatic parts of the PQP+ molecules lie flat in x- and y-directions; however, and in contrast to case (a), these units face each other.

ClO4––ClO4– configuration (see the inset in Figure c): Here, two ClO4– molecules are considered, varying the chlorine–chlorine x-distance, rClCl, while keeping yClCl and zClCl constant.

Face-to-wall topped configuration (see the inset in Figure d): In this case, a single PQP+ molecule, lying flat and parallel to the (x, y)-plane, is located above two layers of Au and is vertically decorated by a ClO4– molecule. The cell geometry is assumed to be periodic in the x- and y-directions and finite along the z-axis; in an effort to scan along the z-direction, we have performed a series of ab initio-based structure optimizations for selected fixed values of zN, i.e., the z-position of the nitrogen in PQP+ above the Au surface. The LJ 10-4-3 potential[53] was used between the Au(111) surface and the molecules.

Face-to-wall beside configuration (see the inset in Figure e): In contrast to case (d), the PQP+ cation is now horizontally decorated by the ClO4– anion such that both molecules are adsorbed on the Au surface. Again, the LJ 10-4-3 potential[53] defined in eq9 was used between the Au(111) surface and the molecules.

In practice, we first optimize , given in eq , involving thereby all interatomic force field parameters; their values are listed in Table for the LJ and the Mie models. These parameters are then kept fixed and are used in the subsequent calculations to optimize the wall force field parameters via optimizing , specified in eq 12; the emerging parameters are listed in Table . In Figure f, we present a visualization of the molecules PQP+ and ClO4–, using these optimized parameters and providing information about the charge of the atomic entities via the color code.

Table 1

Numerical Results for the Optimized LJ and Mie Parameters, and , for Each Element i, for the Results Depicted in Figure (σ in Å and ϵ in meV)a

	σH	σC	σN	σO	σCl	ϵH	ϵC	ϵN	ϵO	ϵCl
	2.243	3.658	3.743	2.865	5.953	3.052	1.204	3.311	7.396	0.172
	2.236	3.703	3.328	2.428	4.956	3.999	0.946	2.021	11.481	5.289

Reference values from the literature are listed in Table S1 in the Supporting Information.

Table 2

Numerical Results for LJ Length and Well-Depth Parameters, σw in Å and ϵw in meV, between the Wall and Each Element i = [H, C, N] and j = [O, Cl], Grouped by the Molecules They Belong to (PQP+ and ClO4–), for Intermolecular Short-Range LJ Parameters Listed in Table and Corresponding σw and ϵw Parameters for Intermolecular Short-Range Mie Parameters Also Listed in Table

	σw[H,C,N]	σw[O,Cl]	ϵw[H,C,N]	ϵw[O,Cl]
	3.197	3.625	3.741	15.781
	3.208	3.630	3.698	20.167

Identifications of Self-Assembly Scenarios

With the classical force field for the PQP+ and ClO4– molecules introduced in Section at hand, we are now ready to identify the ordered ground-state configurations of these molecules as they self-assemble on the Au surface, immersed into an electrolyte and exposed to an electric field. While we leave a more comprehensive and systematic investigation of these self-assembly scenarios to a future publication,[60] we focus in this contribution on the technical details of our approach and on a few selected sets of external parameters (i.e., the electric field strength and the particle density).

Our overall objective is to find for our system the global minimum of the total free energy, F, at T = 0 K as a function of the positions and orientations of all molecules per unit cell for a given value of cell volume and E; at T = 0, this task reduces to the minimization of the internal energy U. The minimum has to be found in a huge-dimensional parameter space, spanning the positions and orientations of the molecules, and by the parameters specifying the unit cell. To be more specific, the dimensionality is set by the number of parameters to be optimized, which read 64, 76, and 88 for five, six, and seven molecules per unit cell, respectively. It is a particular strength of our optimization algorithm (as detailed below) to identify in an efficient and reliable manner minima in such high-dimensional search spaces.

For this purpose, we use a memetic search algorithm that combines evolutionary search strategies (EA)[61−67] and local, steepest gradient descent procedures (LG):[68,69] first, a total number of NEA different lattice configurations,  as defined in eq , is generated. We note that among those configurations we have also intentionally included as “educated guesses” molecular configurations, inspired by the experimental self-assembly scenarios identified in ref (3); however, it should be emphasized that this information is only available for the PQP+ ions, as the experiment does not provide any information about the positions of the perchlorate ions. This population, , is exposed to concepts of natural (or, rather, artificial) selection. At every iteration step of the EA, a new configuration, i.e., an offspring, is created from existing configurations of the most recent population via crossover and mutation operation. This new configuration is then subjected to an LG optimization, an operation that represents by far the most time-consuming task in our algorithm and is performed in parallel using the “mpi4py” framework.[70−72] For an optimal load balance, we additionally spawn a master thread on one of the mpi processes to asynchronously distribute optimization tasks of offspring configurations among all idle mpi processes. The relaxed configurations are gathered by this master thread, which then decides, via a criterion primarily based on the respective internal energy of the configurations, whether the new relaxed particle arrangements are accepted or rejected.

Since the experimental observations[3,5] provide evidence of a structural organization of the molecules into supramolecular lattices, the center-of-mass coordinates of the molecules, R, and their orientations, P, as well as the parameters defining the unit cell, , (see Section and Figure for details), are the variables that have to be optimized for the search of ground-state configurations: we minimize , defined in eq , with respect to R, P, and , keeping the number of molecules N, the unit-cell volume V (with fixed slab width c), and the electrostatic field strength E constant.

In more detail, we proceed as follows:

It is common in evolutionary algorithms to define a genome representation of the entity that is subject to optimization.[73,74] In our case, we represent a supramolecular lattice configuration phenotypically (rather than genotypically[74]) by the set as defined by eq , i.e., the set of all COM coordinates and orientations of all molecules as well as the lattice vectors.

In the first step, two configurations (labeled henceforward by Latin indices), and , are chosen at random or via the “roulette wheel” method (see item (v) below) from the evolutionary population;[61,62,75−77] this strategy favors parents of high quality, hence making them more likely to be used for reproduction than “weaker” configurations, i.e., configurations with higher energy from the evolutionary population. Then, these two configurations are combined via a crossover operation (i.e., a cut-and-splice process; for more details, see below), creating thereby an offspring configuration, , with the subscript “i ⊕ j” emphasizing the executed crossover operation between and . The purpose of this operation is to save high-quality blocks of the genetic material (e.g., the relative positions and orientations of molecules within the unit cell) to efficiently sample the parameter space.[61,62,73−77]

In the second step, the newly generated offspring configuration, , is then exposed to random mutation moves: these are either translations or rotations of single molecules, swaps of center-of-mass positions or orientations of pairs of molecules, or deformations of the unit cell, each of them with a certain probability and within preset numerical boundaries. This step of the algorithm has the purpose of exploring disconnected areas in the parameter space, a feature that is indispensable in global minimization techniques.

After these two steps, and assuming that the offspring configuration, , does not represent a local minimum with respect to the potential energy, a local energy minimization is performed. Here, we mainly rely on the “scipy” implementation of the SLSQP gradient descent algorithm[68,69] (allowing us to define numerical boundaries and constraints on the parameters during the optimization), which minimizes the forces and torques between the molecules as well as the stress of the unit cell. These tools are very helpful to keep the unit-cell volume fixed and to prevent re-orientations of the molecules where some of their atomic constituents would be transferred into positions outside the slab geometry, ensuring thus that z > 0 for all atoms. Subsequently, we perform several “basin-dropping” (BD) steps, where we further try to lower the energy of the configuration by applying several small random “moves” in the parameter space of the LG-optimized offspring; from the emerging configurations, only the ones with low energies are accepted. This specific operation turned out to considerably improve the convergence rate of the local optimization, in particular if multiple and alternating sequences of LG and BD runs are applied.

After the local search procedure, the optimized offspring configuration, , becomes a new candidate to enter the evolutionary population, . The objective of the EA is to retain the best configurations (i.e., the energetically most favorable ones) within the population and to include only the candidates with energy values better or comparable to those of the current population. In an effort to quantify the quality of the candidates, their so-called fitness is evaluated,[61,62,74−77] for which we have used in this contribution the functionF(U) is a monotonic function of the energy U of the candidates, whose value ranges within the interval 0 ≤ F(U) ≤ F(Umin) = 1; Umin and Umax are the minimal and maximal energies appearing in the population, respectively. The selection parameter s quantifies the reproduction rate for configurations within the population in the sense that large values of s tend to exclude configurations with low fitness from reproduction; following ref (77), we commonly use s = 3. The aforementioned “roulette wheel” method for choosing suitable parent configurations also relies on the fitness function (and hence the selection parameter): assuming that the configurations within the population are sorted by their respective fitness values in descending order, F(U) > F(U), the probability, f(U), of a configuration, , to be selected for reproduction is given in terms of the relative fitness[61,75−77]NEA is the total number of configurations within the population. With a certain probability (commonly in 20% of all crossover moves), we allow reproduction between randomly chosen configurations.

Once a new configuration is accepted to enter the population, another configuration has to be eliminated. The probability p(U) for a configuration, , to be replaced is given bya value that is again related to the fitness of the configuration, F(U), and the selection parameter s. Thus, configurations with low fitness are more likely to be eliminated. In any case, a few of the best configurations within the population are retained in an effort to keep the so far best solutions as appropriate parent candidates for the above-mentioned crossover procedures (a strategy referred to in the literature as elitism[75]). It should be emphasized that this strategy does not follow biological selection mechanisms,[78] where populations are replaced entirely once that new generations have been formed; however, our strategy ensures to protect the best genetic material from extinction during the entire search procedure.[75,77]

In an effort to maintain diversity within the population, an additional operation (in the literature referred to as nichening[79,80]), is applied: locally optimized offspring configurations will be discarded if the values of their energy are too close to the energy of any configuration currently in the population, avoiding thereby that the population is overrun by structurally identical configurations. At the same time, the maintenance of genetic diversity is guaranteed.

However, this procedure alone cannot cope with “degenerate” configurations, i.e., if structurally distinct configurations have essentially the same energy values (within the specified nichening tolerance). In our approach, we allow configurations to enter the population only if their structures differ significantly from those of the competing, degenerate configurations. To quantify the structural difference between configurations, we associate a feature vector, f, which collects a set of order parameters pertaining to configuration (see S.I. Section S3.2 for details). The degree of similarity between two configurations, and , is then evaluated by taking the Euclidean distance between the corresponding feature vectors, i.e., Δ = |f – f|; similar configurations will have a small distance, while unlike configurations will have a large distance. If Δ is above a certain threshold value, the offspring configuration, , will not be discarded by the energy-nichening operation.

To offer the reader an insight into the computational complexity of our project, we outline via a few characteristic numbers the computational limitations: the bottlenecks of the identifications of self-assembly scenarios are (i) the huge number of calls of energy evaluations in the optimization step, as underlined via some example: per generation, we have at least 104 calls of the energy kernel; for each state point, we need at least 104 generations, which leads to an absolute minimum of 108 calls of the energy kernel for one (!) set of system parameters and (ii) the optimization of the energy in a high-dimensional search space (as specified above), ranging from ∼60 to 90, depending on the number of molecules.

Summarizing, the complexity of the problem at hand forces us to use all of the above-listed advanced optimization tools, including a basin-hopping memetic approach combining the heuristic nature of evolutionary strategies with deterministic local gradient descent algorithms.[61] The gradient descent method deterministically evaluates every local minimum of the basin with high accuracy (which is additionally sped up by the “basin-dropping” procedure), while the evolutionary search gradually adapts its population to the energetically most favorable solution, exploring the search space for the global optimum.

To round up this section, it should be noted that a variety of techniques have been used in the literature for related optimization problems; among those are Monte Carlo or molecular dynamics-based techniques such as simulated annealing,[81,82] basin-hopping,[83−85] minima-hopping,[86,87] and eventually evolutionary approaches such as genetic algorithms.[61−67,79,80,88−91] The decision on the method of choice relies on the specific problem: for instance, as Hofmann et al. used the SAMPLE technique (see refs (7, 21) and (22)), relying on a discretization of the search space into limited, archetypical, intermolecular motives and elaborate data fitting of emerging force fields to describe intermolecular interactions. To the best of our knowledge, this approach has neither been applied to molecular motives beyond monolayer configurations or to charged molecules, so far, nor has it been used in combination with an external control parameter, such as an electric field or systems composed of multiple components. In general, the fact that the number of archetypical intermolecular motives grows rapidly with the increasing size of the molecules bears the risk of hitting very soon the limits of computational feasibility. However, suitable adaptations of this strategy and/or a combination with evolutionary search strategies or with reinforcement learning, which has, for instance, very successfully been applied to protein folding problems[92] in a similar way as AlphaGo[93] was able to master the infamous board game, might represent a viable route to circumvent the aforementioned limitations; thus, future investigations of such intricate problems as the complex monolayer-to-bilayer transition, addressed in this contribution, might come within reach.

Results

General Remarks and System Parameters

Below, we present selected results for self-assembly scenarios of PQP+ and ClO4– molecules on an Au(111)–electrolyte interface under the influence of an external electrostatic field, as obtained via the algorithm presented in the preceding sections. Our choice of parameters is guided by the experimentally observed molecular configurations.[3] We demonstrate that our proposed strategy is indeed able to reproduce on a semiquantitative level the experimentally observed self-assembly scenarios.[3] As a consequence of the still sizable costs of the numerical calculations, we leave more detailed investigations (where we systematically vary the system parameters) and a quantitative comparison of our results with the related experimental findings[3] to a future contribution.[60]

To be more specific, we have used the following values for the (external) system parameters:

An indication for the number of molecules per unit cell is provided by the experiment:[3] we have considered unit cells containing 10, 12, and 14 pairs of PQP+ and ClO4– molecules. These numbers of molecules include, of course, also the related mirror molecules and correspond to 630, 636, and 742 atomic entities per unit cell, respectively (which interact via short-range and long-range potentials, which are subject to particle wall interaction and are sensitive to an external electrostatic field).

Also, the actual values of the surface area A are motivated by estimates taken from the experiment:[3] we have varied A within the range of 6.5–12.25 nm2, assuming a step size of typically 0.5 nm2; systems will be characterized by the surface density of the PQP+ molecules, defined as σPQP = NPQP/A, with NPQP being the number of PQP+ molecules per unit cell.

The range of the experimentally realized values for the electrostatic field strength E is, however, difficult to be estimated since the major drop in voltage occurs near the negatively charged Au surface and the nearby layers of cations,[54] which is not directly accessible in the experiment. Therefore, we have covered, at least in this first contribution, several orders of magnitude in the value for E within a range that extends (on a logarithmic grid) from E = −10 to −10–3 V/nm; in addition, we have also performed calculations at zero electrostatic field.

It should be mentioned that we have used in all of these calculations the Mie potential within the classical model, since the related LJ model is not able to fit the ab initio data with a comparable and sufficient accuracy (see also the discussion in Section ).

We have covered in total approximately 176 combinations of these parameters (that is the unit-cell volume V with a constant slab width c, the number of molecules N, and the electrostatic field strength E); for each of these, we performed independent evolutionary searches with a population size of typically configurations. Some details about the numerical costs of our calculations can be found in S.I. Section S4.1.

Discussion of the Results

Lateral Particle Arrangements

In this section, we discuss the lateral self-assembly scenarios of the PQP+ and of the ClO4– molecules. Selected results for our numerical investigations are presented in Figure and in, on a more quantitative level, in Table . The actual values have been chosen in an effort to reproduce, at least on a qualitative level, the results obtained in the experimental investigations. Indeed, the sequence of the obtained ordered ground-state configurations (shown in panels a–c) clearly indicates the transition from a stratified bilayer configuration (identified at a rather strong electrostatic field strength of E = −0.3 V/nm) to a self-host–guest monolayer structure (obtained by reducing the field down to E = −0.1 V/nm) and eventually to an open-porous configuration (identified at E = −0.01); similar observations have been reported in the related experimental study.[3]

Figure 4

Results for the ground-state configurations of PQP+ and ClO4– molecules, adsorbed on a Au(111) surface under the influence of an external electrostatic field E, as they are obtained via the numerical procedure, as specified in Section ; calculations are based on the classical model for the molecules, involving the Mie potential (for details, see Section ). In the main panels, configurations are shown in a periodically extended view as projections onto the (x, y)-plane and in the respective insets as projections onto the (y, z)-plane; in the main panels, the respective unit cells are highlighted by thick black lines. Results are shown for different values of the number of PQP+ molecules, the surface density σPQP, and the electrostatic field E: see labels in different panels and Table for details. The red (gray) shaded areas, framed with dashed lines, in (a) and (d) emphasize PQP+ molecules that sit on top of other cations, forming a bilayer structure. The dashed, shaded, magenta (gray) rectangular, and green (gray) square areas in (b) represent tilings formed by perchlorate molecules within the dense PQP+ monolayer configuration. The dashed, shaded, and blue (gray) circles in (c) and (e) emphasize, quantitatively, the porous and auto-host–guest tiles identified in the experiment (see Figure A,C in ref (3), respectively).

Table 3

Results of Evolutionary Ground-State Search for Different Electric Field Strengths, E, for Different Unit Cell Areas, A, and Number of PQP+ Molecules, NPQPa

	Ez (V/nm)	U/NPQP (eV)	NPQP	A (nm2)	NPQP/A (nm–2)
(a)	–0.30	–1.5804	6	8.5	0.705882
(b)	–0.10	–1.7276	6	11.75	0.510638
(d)	–0.10	–1.7274	5	8.0	0.625000
(c)	–0.01	–1.6445	6	11.25	0.533333
(f)	–0.01	–1.6364	6	12.25	0.489796
(e)	–0.01	–1.6332	7	11.75	0.595745

Each line represents an evolutionary search. The respective structures are presented in Figure .

From the results of our investigations (which are shown only selectively), we learn that an electrostatic field strength of E = −0.3 V/nm always leads to bilayer configurations, similar to the one shown in Figure a. This stratified bilayer configuration represents the energetically most favorable one as we vary at fixed E the volume of the unit cell and the number of molecules within the respective ranges, specified in the preceding section; the numerical data of the related internal energy are compiled in Table .

As we proceed to E = −0.1 V/nm, we observe self-assembly scenarios as the ones depicted in Figure b,d, which correspond to self-hosts–guest configurations observed in the experiment;[3] for the data presented in these panels, two different values for NPQP (and hence for σPQP) have been considered: the monolayer configuration shown in panel b has a slightly better value for the internal energy (per molecule) than the rhombohedral bilayer configuration shown in panel d; however, as can be seen in Table , the energy differences are very tiny: differences of the order of 10–4 eV correspond to values where we hit the numerical accuracy of the ab initio-based energy values.

Eventually, we arrive at the so-called open-porous structures, observed in the experiment:[3] the ground-state configurations depicted in Figure c,e,f are evaluated at the same electrostatic field strength of E = −0.01 V/nm, assuming different values for NPQP and σPQP; the open-porous pattern emerging in panel c is the most favorable one in terms of energy per molecule (see Table for the numerical details). There are, however, several serious competing structures with minute energy differences at this value of the electric field strength: another open-porous structure, depicted in panel f, with an energy penalty of less than 8.1 meV per PQP+ molecule compared to case (c) and also a considerably denser configuration, depicted in panel e, with an internal energy value worse by only 11.3 meV compared to case (c) and by 3.2 meV compared to case (f).

From the numerical point of view, the following comments are in order: for a fixed state point, the energy differences of competing structures attain values that hit the limits of the accuracy of the ab initio-based simulations, which can be estimated to be of the order of 0.1–0.01 eV per molecule for dispersive interactions.[35,94−96] These values set the limits of our numerical accuracy. For completeness, we note that for the results for the energies obtained via the classical force field (which are based on LAMMPS calculations), we estimate that our results are numerically reliable down to ∼10–6 eV per atom; within the range of such minute energy differences, no competing structures have been found in our investigations. In general, we observe that the energy differences for the energetically optimal ground-state configurations become smaller as the electrostatic field tends toward zero. Even though the optimization algorithm (as outlined in Section ) has turned out to be very efficient and reliable, we observe (in particular for smaller values of the external field) that new configurations are included in the population of the best individuals even after a large number of optimization steps.

Vertical Particle Arrangements

In Figure , we present in separate panels the height distributions of PQP+ and ClO4– as functions of the electrostatic field, E (which is binned for the six different values of E that were investigated); along the vertical axis, we count (for a given value of E) the occurrence of the respective molecules in bins of 1 Å, and normalize by the total number of all considered configurations identified by the evolutionary algorithm, which are located within an interval of at most 1 kBT (or 43 meV) above the configuration with the best energy, which is of the same order of magnitude as the values presented in Table for different electric field strengths, E. Note in this context that the distance of the first layer of PQP+ molecules can be directly estimated by the vertical equilibrium position of carbon atoms, zC(eq) = 3.166 Å, obtained by minimizing eq for a single carbon atom with σwC(Mie) = 3.208, taken from Table .

Figure 5

Height distribution of PQP+ (left) and ClO4– (right) molecules as functions of the considered values of the electrostatic field, E, normalized to the number of respective molecules (see the color code on the right-hand side of the panels; the value of one means that all respective molecules in all considered configurations are counted in one specific bin; see the text for the energy considered in this analysis). Along the vertical axes, the binning is performed in steps of 1 Å: z = 0 Å marks the position of the gold surface, and the slab width amounts to 12 Å.

For the high values of the field strength (i.e., for E = −10 and −1 V/nm), the PQP+ ions are preferentially adsorbed onto the gold surface as a closely packed monolayer (see the left panel of Figure ), while the perchlorate anions are strongly dissociated and assemble as far from the gold surface as possible (corresponding in our investigation to the numerical value of the slab height, which we fixed to 12 Å); see the right panel in Figure . This situation represents an extreme case in the sense that neither the Coulomb nor the short-range Mie interactions between ions and anions can compensate for the strong negative surface potential.

Decreasing now the magnitude of the electrostatic field to the more moderate value of E ∼ −0.3 V/nm reveals the emergence of a bilayer structure, formed by the PQP+ molecules, with pronounced peaks located at zPQP(1) ∼ 3 Å and zPQP(2) ∼ 7 Å, with relative weights of 72 and 21%, respectively. Such stratified bilayer configurations (as depicted in Figure a) are in competition with structures similar to ones shown in Figure d. Note that in parallel a far more complex height distribution of the perchlorate molecules sets in as soon as the now moderate electrostatic field allows them to proceed toward the interior of the slab: now, more than half of the ClO4– anions are located “in between” the PQP+ “layers”, trying on the one side to compensate the charges of one or several PQP+ “partners” in the slab region and “filling spatial holes” wherever they can, on the other side. A large portion of perchlorate ions is even allowed to adsorb on the surface at a distance of zClO(1) ∼ 4.074 Å; note that these COM positions above the interface are larger than the minimal height of the PQP+ cations due to two reasons: if one face of the oxygen tetrahedron is oriented toward the interface (i.e., parallel to the gold surface), the COM of the ClO4– ion is increased by a value of zCl – zO ≈ 0.492 Å, with respect to the oxygen atoms. These atoms themselves have an equilibrium distance to the surface of zO(eq) ≈ 3.582 (evaluated by minimizing eq for a single oxygen atom with σwO(Mie) = 3.630, cf. Table ), summing up to the presented minimal COM distance of the adsorbed ClO4– molecules from the interface. Note that the height distribution of the ClO4– ions is now rather broad (see Figure ), which is definitely attributed to their relatively smaller size and their considerably higher mobility, as compared to their cationic counterparts (see also the discussion below); these features make a conclusive interpretation of the roles of the ClO4– ions in the structure formation of the entire system rather difficult. Our interpretation is that the perchlorate ions are, due to their small spatial extent and their high mobility, able to compensate for local charge mismatches and to act as spatial spacers between the cations.

Figure 6

Four structurally different configurations of perchlorate ions (framed by gray, dashed lines and differently shaded areas) at optimized, fixed cell geometry (indicated in the bottom-left corner by the black, dashed line) and optimized, fixed positions and orientations of the PQP+ ions; the energies of these four configurations range within an interval of 38 meV (per PQP+–ClO4– pair), as obtained after an evolutionary energy minimization of solely the degrees of freedom of the ClO4– ions, starting from the configuration depicted in Figure a. Changes in the structure as one proceeds from the left to right are highlighted by respective circles (specifying the position of the “moving” perchlorate ion) and arrows.

Decreasing further the magnitude of the electrostatic field down to E ∼ −0.1 and −0.01 V/nm provides unambiguous evidence that the formation of the PQP+ ions into bilayer structures become energetically more and more unfavorable, as the upper peak in the height distribution of the cations vanishes gradually. Concomitantly, an increasing number of perchlorate molecules approach the gold surface and are predominantly located there; possibly, they act as a space filler on the surface itself, while at the same time the small values of the electrostatic field keep the PQP+ molecules near to the surface. In this context, it should be noted that decreasing the electrostatic field is equivalent to decreasing the surface potential; thus, and in combination with the adsorbed perchlorate molecules, a transition from an auto-host–guest to a porous structure is plausible.

Eventually, at zero electric field, the system exclusively gains energy from intramolecular interactions and adsorption on the gold surface. Since the perchlorate molecules are rather spherical in their shape, they can efficiently adsorb onto the gold surface (in an orientation explained above), while the PQP+ molecules are able to efficiently stack, especially without a guiding electrostatic field.

Role of the Perchlorate Anions

We come back to the above-mentioned volatility of the perchlorate ions: in Figure , we present results from yet another evolutionary analysis, i.e., we now fix the positions and orientations of PQP+ molecules as well as the extent and the shape of the unit cell of some optimized configuration (as, for instance, depicted in Figure a) and vary only the degrees of freedom of the perchlorate anions. Figure shows, for fixed cell geometry and fixed PQP+ positions and orientations, four structurally different perchlorate arrangements whose energy ranges within an interval of 38 meV (per PQP+–ClO4– pair): the fact that we obtain completely different configurations of the perchlorates (with essentially comparable energies) undoubtedly indicates the high mobility of the ClO4– ions. Changes in the structure, as one proceeds from the left to right, are highlighted by respective circles (specifying the position of the “moving” perchlorate ion) and related arrows. The ClO4– molecules exhibit a remarkable freedom in their rotation without (or only marginally) changing the energy of a configuration; this fact has rendered the minimization of the energy very difficult. However, it should also be noted that even translations can be performed without a substantial change in energy.

We note that the analysis of these different structures was achieved using a so-called t-SNE[97] analysis on the leading five principal components of a principal component analysis (PCA)[98] of order parameters of all configurations identified by the evolutionary algorithm; for more detailed information on this rather technical issue, we refer to Figure S5 in Section S3.2 of the Supporting Information.

Conclusions and Outlook

The prediction of the supramolecular ordering of complex molecules at a metal–electrolyte interface using DFT-based ab initio calculations is in view of the expected gigantic computational costs, and despite the availability of petascale computers, still an elusive enterprise. In this contribution, we have proposed a two-stage alternative approach: (i) DFT-based ab initio simulations provide reference data for the energies introduced in a classical model for the molecules involved, where each of their atomic entities are represented by a classical, spherical particle (with respective size, energy parameters, and charges). We modeled the interaction between the atomic entities and the metallic surface by a classical, perfectly conductive, Lennard-Jones-like wall potential; the electrolyte is treated as a homogeneous, dielectric medium.

The interparticle and particle–wall parameters were obtained via the following procedure: considering archetypical configurations (involving pairs of ions and/or ions located close to the surface), DFT energies were fitted by the related energy values of the classical model. (ii) The second step identifies the ordered ground-state configurations of the molecules by minimizing the total energy of the now classical system. This optimization is based on evolutionary algorithms, which are known to operate efficiently and reliably even in high-dimensional search spaces and for rugged energy surfaces.

Our new two-stage strategy overcomes the hitherto prohibitive computational cost of modeling the full system while reproducing the key observations of a well-documented experimental system consisting of disc-shaped PQP+ cations and ClO4– anions: as a function of increasing electric field at the metal–electrolyte interface, the molecular building blocks are seen to self-organize into an open-porous structure, a self-host–guest pattern, and a stratified bilayer. Future work will focus on verifying the extent of the predictive power of our method toward molecular self-assembly under electrochemical conditions and on strategies to further streamline and reduce the computational cost of our approach, without sacrificing the reliability of the predicted results.

In view of the high computational costs and the conceptual challenges encountered in our investigations, we have pondered the question whether the complexity of the current model (which, as a classical model, is comprehensive in the sense that it contains all atomistic features) could possibly be further reduced, avoiding thereby conceptual and computational bottlenecks. The idea behind this strategy is to develop, starting from the present model, a hierarchy of ever simpler models where, for instance, larger subunits of the molecule (such as aromatic rings) are replaced by disk-shaped units carrying higher electrostatic moments, as schematically visualized in Figure . Such a model might provide a first, semiquantitative prediction of the self-assembly of the PQP+ and of the ClO4– ions at considerably reduced costs and might help in prescreening possibly promising portions of the huge parameter space for subsequent investigations of the full model. Efforts in this direction are currently pursued.

Figure 7

(a) Atomistic model of the PQP+ molecule as used in this contribution (white: hydrogen atoms, gray: carbon atoms, and blue: nitrogen atom); see also Figure S3 in the S.I. (b) Related coarse-grained model in a hierarchy of ever simpler models, using, e.g., Gay–Berne potentials to account for the van der Waals interaction of all atoms in the specific rings and a multipole expansion to second order (monopole as colored points, dipole moments as small arrows) for the electrostatic interaction.