Combined DFT and Molecular Mechanics Modeling of the Adsorption of Semiconducting Molecules on an Ionic Substrate: PTCDA and CuPc on NaCl.

Experimental results suggest that molecular geometry and energies can be influenced by the presence of thin film substrates as well as surrounding molecules. It is imperative that computational models take this influence into account. The accurate computational modeling of these molecules is an efficient way of carrying out chemistry calculations and reinforcing experimental findings. In our study, density functional theory (DFT) and molecular mechanics (MM) are used to model the configurations of the organic semiconducting materials, 3,4,9,10-perylene tetracarboxylic dianhydride, C24H8O6 (PTCDA), and copper(II) phthalocyanine, C32H16CuN8 (CuPc), as adsorbed on single- and double-layer NaCl substrates of various dimensions and charge settings. After a geometry and charge optimization of the molecules using DFT, the molecular geometries are optimized under different environments using computational calculations with specific force-field settings in HyperChem Professional 8.0(TM) software using MM. Energies and geometries of the molecules are then recorded, and our data are compared to experimental results of similar systems. We find that, with the appropriate choice of substrate properties, the calculated molecular configurations directly reflect those found experimentally. Our results support the idea that this method of simulation can produce reliable models in the field of physical chemistry.

Introduction

The prototypical semiconducting molecules 3,4,9,10-perylene tetracarboxylic dianhydride (PTCDA) and copper(II) phthalocyanine (CuPc) (see Figure ) have garnered great interest in the fields of physical chemistry and materials science. Because of their optical and photovoltaic properties, these molecules have potential in a variety of electronic applications, such as organic field-effect transistors, light-emitting diodes, solar cells, and molecular switches.[1−6] Extensive studies have been completed on the morphology of these molecules on various kind of substrates, including noble metals, graphene, and insulating substrates.[7−13] The molecular arrangement and adsorption to the substrate is of particular interest when stable low-energy configurations are studied. The ability to predict the geometries and the respective energies of these organic semiconducting molecules is critical to the manufacturing of devices and photovoltaics, which emphasize lightweight and low-cost bulk materials.[14,15]

Figure 1

Lewis structures of CuPc and PTCDA.

Because of adverse charging effects and interference with electrons in bulk conducting surfaces, the observation of individual organic semiconducting molecules on inorganic insulator surfaces provides new avenues to study adsorption behavior. Insulating materials are incredibly useful in relation to device manufacturing for the purpose of isolating device regions.[7] Research has shown that insulating ultrathin films such as KBr, KCl, or NaCl can be used as a surface for the adsorption of organic molecules.[16] These configurations limit charging effects; however, when, for instance, imaged with scanning tunneling microscopy (STM) the electrons can tunnel through the the insulating film to an underlying metal substrate.[17−19] The geometric change and subsequent binding energies of organic molecules can be studied as they adsorb to an insulator.

It is important to study how these semiconducting molecules behave in bulk structures. This research can be used in the fields of device processing and materials, where an understanding of interactions between molecules is imperative. Experimental work has been done to study the behavior of these molecules individually as well as in bulk settings; however, computational models of these systems have not been very extensive so far.[20]

Although computational models are common, many methods and software can prove to be expensive and inaccessible.[21] Our goal is to add supplementary data to the experimental findings without losing accuracy of observation or increasing the complexity of the process. Here, we use specific computational methods to add to recent experimental work.[20] Our computational software for molecular mechanics (MM) calculations, HyperChem 8.0, and for density functional theory (DFT) to obtain partial atomic charge, ORCA, can run calculations on a standard desktop computer system.[22−24] One factor for the accuracy of MM is the individual charge configuration of the molecule and the substrates that we have investigated in this work. Molecular mechanics is especially suitable in this context since the ionic nature of the substrate allows a straightforward assignment of partial atomic charges to these templates. The determination of the atomic charges of the semiconducting molecules can be carried out using computational tools employing density functional theory. Comparing results of computational chemistry methods performing geometry optimizations to the experimental work successfully allows us to understand the resilience and legitimacy of our computational methodology. This approach provides an avenue for future combinations of experimental and computational analysis of molecular adsorption on insulating substrates.

Computational Approach

PTCDA and CuPc were modeled using HyperChem 8.0.[24] Once the bonds and atomic locations were optimized, the coordinates of HyperChem-optimized molecules were used as input for the ORCA calculations in order to determine partial atomic charges using DFT.[22,23] We optimized our molecules with the DFT B3LYP/G method, a hybrid generalized gradient approximation (GGA), which is well-established for its robustness and speed. The self-consistent field (SCF) convergence was set to the default for ORCA as SCF = Tight. In order to keep our approach the most versatile, we focused on the Mulliken charges that are generated as one result of the DFT optimization. Experimental studies of both CuPc and PTCDA have suggested various possible charge states when adsorbing on a NaCl/metal template.[25−27] However, in this study we stick with a neutral configuration, as this situation allows for the most general applicability with respect to the underlying substrate. These molecules with newly assigned partial charges were then exported back into HyperChem for another geometric optimization using the Assisted Model Building with Energy Refinement (AMBER 3) force field, which took into account the Mulliken charge assignments from the ORCA software.

Hyperchem files/substrates of NaCl were computationally fabricated using a bulk lattice constant of a = 5.64 Å.[28] The substrates we created function as a local environment for the adsorption of the molecules. As mentioned before, molecular mechanics uses partial atomic charges in order to calculate interaction energies and geometries. Therefore, the configuration of the substrate is important in these adsorption processes. Previous works with different approaches have suggested three possible charge assignments for the ionic NaCl substrate.[29−31] The Pauling charge assignments are based on empirical data according to the formula: Percent Ionic = 100% ·(1 – e – 1/4(χa – χb)2), with χ as the electronegativity of the two ionic components. This formula gave an ionization of ±0.67e for the ions in the NaCl substrate.[29] Other charge settings used for separate substrates are the density derived electrostatic and chemical (DDEC) charge setting and the plane wave charge setting.[30,31] The DDEC setting gave an ionization of ±0.94e, and the plane wave gave ±0.42e, respectively, for the individual ions in the NaCl substrate. The sodium atoms were assigned a positive ionization, and the chlorine atoms were assigned a negative ionization. For each charge type, a single-layer NaCl substrate and a double-layer NaCl substrate were created. To explore possible effects of the edge of the substrate related to the long-range interactions of adsorbate and substrate, we created different lateral substrate sizes: 8 × 8, 16 × 16, 32 × 32, 64 × 64, and 128 × 128 (each dimension is in units of atoms).

To preserve the integrity and consistency of the data, we must avoid our molecules being influenced by the boundaries of the substrate. We ran models for the adsorption of the molecules in vacuo mimicking adsorption starting with the molecule at an initial height of 4 Å above each substrate using an AMBER3 molecular mechanics force field and a convergence limit of 0.001 kcal/mol. AMBER force fields were often used in the study of both synthetic and naturally occurring polymers.[32,33]

The molecules were individually adsorbed in 15° rotation intervals about the z-axis (the z-axis is perpendicular to the substrate), with a single point energy calculation recorded after each adsorption. Different configurations and adsorption energies occur as the molecule was rotated about the z-axis and bound to the surface of the NaCl substrate. This molecular rotation is imperative in finding the global minimum energy of the adsorbed molecule. Often, when a molecule is adsorbed onto a substrate, a local minimum-energy configuration will be found even with a very low energy gradient. When the molecule was rotated on the substrate, energies were compared to find a global minimum-energy configuration. The lowest energy state of each rotation was recorded, and the lowest single point energy was considered the optimal configuration of the molecule. Finding a global minimum in energy as a function of configuration yields results about optimal geometries and locations for the adsorbed molecules.

Results and Discussion

After comparing adsorption energies and geometries for all different sizes of substrate, the 64 × 64 atom substrate (with lateral dimensions of ∼175 Å × 175 Å) was deemed to be the ideal size for the adsorption of the molecules based on the consistency of data, suggesting a lack of influence from the environmental boundaries. Data for the change in adsorption energies when increasing the substrate size for CuPc and PTCDA can be seen in Figure . One can see an exponential change of the differences when moving from smaller to subsequent larger substrates. However, when changing from the 64 × 64 atoms substrate to 128 × 128 atoms substrate, the change in energy is either below or on the order of the convergence limit of 0.001 kcal/mol. Therefore, the 64 × 64 substrate was chosen for all subsequent calculations.

Figure 2

Adsorption energies of various-sized substrates were compared. The energy difference (in log scale) is plotted against the change for the substrate. The energy convergence for the adsorption is set to 0.001 kcal/mol. When the size of the substrate is increased from 64 × 64 to 128 × 128, the change in adsorption energy is on the order (for PTCDA) or smaller (for CuPc) than the convergence limit. Because of the dramatic increase in computational time for a 128 × 128 substrate and no significant change in adsorption energy and especially geometry, the 64 × 64 substrate was deemed suitable without loss of accuracy.

For each type of 64 × 64-atoms substrate, the molecules were individually adsorbed using the same parameters and methods as used before for the determination of the ideal substrate. Energy calculations were recorded each time, and the lowest energy result was considered as the ideal configuration. Using the data we collected for the adsorption energies of each rotation for both molecules, we captured images of the adsorbed molecules at their optimal configurations (see Figures and 4). Additionally, the recorded energies for the adsorbed geometries were subtracted from energies of a molecule/substrate system where the two entities do not interact; for example, the molecule is far away from the substrate (200 Å above the substrate where no discernible interactions are calculated). These energy differences are summarized in Table and represent adsorption energies for the molecules on the various NaCl substrates.

Figure 3

CuPc geometries on single- and double-layer NaCl substrates (white: Cl, blue: Na) with three differing charge assignments, Plane Wave (±0.42e) (top), Pauling (±0.67e) (middle), and DDEC (±0.94e) (bottom). (left column) Single-layer NaCl. (right column) Double-layer NaCl.

Figure 4

PTCDA geometries on single- and double-layer NaCl substrates (white: Cl, blue: Na) with three differing charge assignments, Plane Wave (±0.42e) (top), Pauling (±0.67e) (middle), and DDEC (±0.94e) (bottom). (left column) Single-layer NaCl. (right column) Double-layer NaCl.

Table 1

Adsorption Energies of PTCDA and CuPc on NaCl Substrate, Single Layer (SL), and Double Layer (DL) with Various Charge Schemes

	PTCDA		CuPc
	SL	DL	SL	DL
Plane Wave (0.42e)	1.54 eV	1.72 eV	1.90 eV	2.14 eV
Pauling (0.67e)	1.67 eV	1.85 eV	2.13 eV	2.23 eV
DDEC (0.94e)	1.82 eV	2.01 eV	2.27 eV	2.34 eV

Looking at the adsorption geometries for CuPc on the various NaCl substrates in Figure , two configurations can be observed: (I) one with the central Cu atom of the molecule above a Na+ ion in the substrate, and (II) with the central Cu atom above a Cl– ion of the substrate. Additionally, on the one hand, configuration I leads to an alignment of the symmetry axes of the CuPc with the ⟨100⟩ (diagonal) and the ⟨110⟩ (horizontal and vertical) directions, respectively. On the other hand, configuration II with the Cu atoms above the Cl– ion, the symmetry axis of the molecule through the benzene ring is at an angle of ∼15° c with respect to the horizontal ⟨110⟩ direction of the substrate. Comparing the adsorption energies in Table , one can see that the energies increase with increasing polarity of the substrate and also that the double layer of NaCl produces higher binding to the molecule, both of which are to be expected for the system.

For PTCDA on the various NaCl substrates in Figure , a consistent adsorption geometry could be observed. The central benzene ring of the molecule is centered above a Cl– ion in the substrate. The symmetry axes of the PTCDA molecule are nearly parallel with the ⟨110⟩ directions (horizontal and vertical) of the NaCl substrate. Regarding the adsorption energies (Table ), one can observe a similar trend as in CuPc with increasing value for increasing polarity of the substrate and higher energies for the double layer NaCl.

One of the goals of the present project is to support experimental work without the cost-prohibiting use of a full ab initio DFT calculation and instead to use readily available and relatively fast computational methods such as molecular mechanics. One important input into these MM calculations, though, is an appropriate assignment of partial atomic charges. This can be accomplished for smaller molecules, such as PTCDA and CuPc, using the DFT methods we employed here; however, it is prohibitive for the substrate consisting of ∼4000 atoms. Therefore, a suitable choice for ionic templates must be made. Comparing our results with experimental results, such as the scanning tunneling microscopy study of PTCDA and CuPc on a double layer of NaCl, allows one to make inferences about the choice of substrate and the computational accuracy.[7,20,34]

As mentioned above, the adsorption of the PTCDA molecule on the various NaCl substrates resulted in practically identical configurations. This geometry with the central benzene ring located above the Cl– ion and the symmetry axes of the molecule parallel to the ⟨110⟩ axes of the NaCl substrate agrees perfectly with the experimental results of PTCDA on a NaCl layer.[7,20] Therefore, we can not make any comments on the computational choice of substrate configuration based on the PTCDA molecule alone. However, as noticed for the CuPc molecule, two configurations emerged from our calculations. Experimental data suggest that the central Cu atom of the molecule sits above a Na+ ion,[20,34] as was found in our calculation for the single-layer NaCl substrate using either the Pauling method or the DDEC method to assign charges to the substrate. This suggests that the substrate polarity should be higher than suggested by the Plane Wave method, which only assigns a partial charge of ±0.42e to the ions in the substrate. Another important point to make here is the difference between a single layer and a double layer in our calculations, compared to the double layer of NaCl in the experiments. In the aforementioned STM studies, NaCl is grown as a thin film on metal substrates. That means the polarity in the film’s first layer (adjacent to the metal substrate) is influenced by the electrons in the substrate. In our case of the computational approach, this metal substrate is not present, which means the first layer (in a double-layer setup) has the same polarity as the second layer. Therefore, it might be more insightful for calculations to just consider the top layer of NaCl, if it is computationally not suitable to include the metal substrate below the NaCl layer as in our case. Furthermore, our calculations (SL Pauling or DDEC) and experimental STM results[34] suggest a fourfold rotational symmetry, again corroborating the experimental results.

An adsorption model for both molecules that is self-consistent and supports previous experimental work can be seen in Figure . A single layer of NaCl with a charge scheme according to the Pauling approach was chosen as the substrate. The PTCDA molecule is situated on the substrate with the central ring above a Cl– ion and its symmetry axes parallel to the ⟨110⟩ directions of the substrate. The CuPc molecule adsorbs with its central Cu atom above a Na+ ion of the substrate, and the molecule shows fourfold rotational symmetry with the symmetry axes of molecule and substrate aligned.

Figure 5

Computational adsorption geometry with a single-layer NaCl substrate with partial atomic charges of +0.67e for Na and −0.67e for Cl (according to Pauling[29]). (left) PTCDA. (right) CuPc. (center top inset) Experimental STM image and deduced geometry. Reprinted from Cochrane et al.[20] Copyright 2018 American Chemical Society.

Besides comparing the adsorption geometries with experimental values, we also looked at the validity of our results with respect to previous ab initio calculations on similar systems. The adsorption geometries as depicted in top views in Figure for CuPc and PTCDA determined by our methodology are verified by other studies using first-principles calculations.[35,36] Additionally, the height of the molecule above the substrate is another geometrical factor comparing the results of various studies. Figure shows this configuration for our results. We determined distances of 3.20 Å for PTCDA and 3.22 Å for CuPc between the substrate atoms and the central portion of molecules. For PTCDA, Hochheim and Bredow in their DFT study report 3.37 Å,[37] whereas Burke et al. in an experimental work report a distance of 3.4 ± 0.2 Å.[7] Doležal et al. in their combined experimental and calculation work indicate an adsorption height of ∼3.1 ± 0.1 Å for a CuPc molecule on NaCl.[36] All these values are very similar to the ones found by our approaches, therefore validating our results further. We note here that most first-principles calculations allow the substrate atoms to relax as well, whereas in our study the substrate geometry is fixed. That leads to a slight vertical shift of some individual substrate atoms of ∼0.1 Å up or down, which would not substantially alter the distances reported here.[35,37]

Figure 6

Side view of adsorption geometry with a single-layer Pauling (±0.67e) NaCl substrate. (left) PTCDA. (right) CuPc. The distances between substrate and molecules are dPTCDA = 3.20 Å and dCuPc = 3.22 Å.

Although most experimental studies including the aforementioned STM studies do not provide adsorption energy values, we can compare our energy values to other computational approaches of similar systems. When using DFT calculations, Aldahhak et el. found adsorption energies for PTCDA on a NaCl substrate of ∼2.8 eV, whereas Hocheim and Bredow report ∼2.6 eV.[35,37] Both of these values are higher than the ones we summarized in Table ; the difference between DFT and MM can most often attributed to the different approaches and therefore contributions, in addition to the electrostatic part, of van der Waals interactions.[35] Although the CuPc on the NaCl system has been investigated experimentally and computationally with various methods, no literature values for adsorption energies were reported.

Conclusions

Computationally studying the geometries and interactions of organic semiconducting molecules on inorganic insulating substrates can support and enhance experimental work. A variety of comparisons can be completed to verify computational results correspond to the lab-based work. When one is interested in reinforcing experimental results, this computational method can provide calculations that reflect experimental data. Our approach of using DFT-supported molecular mechanics calculations has two advantages with respect to first-principles ab initio approaches without sacrificing accuracy: MM is easier to access and computationally cheaper for experimental researchers, and it allows one to use a larger substrate with more atoms in order to avoid edge effects.

After a careful consideration of the adsorbed geometries, we conclude that a single-layer NaCl substrate with a partial atomic charge of at least ±0.67e (e.g., according to Pauling or the DDEC approach) is the environment that best reflects the experimental results. Our specific computational process is promising for future work in the realm of molecular mechanics research. From the detailed similarities between the molecules adsorbed onto the SL NaCl substrates with Pauling or DDEC and the experimental data, we gain confidence in this particular method of optimizing these molecules and adsorbing them onto NaCl substrates.

Future work would involve an investigation of other molecules adsorbed on ionic substrates and expanding the choice of partial atomic charges also for the adsorbates.