Ab Initio Modeling of Mixed-Dimensional Heterostructures: A Path Forward

Understanding the electronic structure of mixed-dimensional heterostructures is essential for maximizing their application potential. However, accurately modeling such interfaces is challenging due to the complex interplay between the subsystems. We employ a computational framework integrating first-principles methods, including GW, density functional theory (DFT), and the polarizable continuum model, to elucidate the electronic structure of mixed-dimensional heterojunctions formed by free-base phthalocyanines and monolayer molybdenum disulfide. We assess the impact of dielectric screening across various scenarios, from isolated molecules to organic films on a substrate-supported monolayer. Our findings show that while polarization effects cause significant renormalization of molecular energy levels, band energies and alignments in the most relevant setup can be accurately predicted through DFT simulations of the individual subsystems. Additionally, we analyze orbital hybridization, revealing potential pathways for interfacial charge transfer. This study offers new insights into hybrid inorganic/organic interfaces and provides a practical computational protocol suitable for scaled-up studies.


Theoretical Background and Computational Details
Total-energy calculations for the free-base phthalocyanine (H 2 Pc) molecule are performed with density functional theory (DFT), 1,2 employing a modified copy of version 12.2 of Octopus, 3 using the PBE functional 4 and SG15 pseudopotentials. 5This code represents wave functions and density on a real-space grid, sampling with a spacing of 0.2 Å a simulation box constructed as the union of atom-centered spheres of radius 5 Å.The interface to the polarizable continuum model (PCM) is established by the construction of another sphere-union box, using atom-specific radii (van-der-Waals radii scaled by a factor of 1.05), whose surface is discretized into elements. 6Each of these elements carries a charge representing part of the polarization charges in the environment.The charge on each element is, in turn, computed from the molecular charges, convolving their electrostatic potential at the boundary elements with the PCM response function. 6For layered substrates, this response function contains an image term for substrate polarization as implemented in LayerPCM, 7 an extension of PCM for layered environments.
All simulations employing PBE0 8 and HSE 9,10 functionals, as well as all DFT simulations of MoS 2 , are conducted with the Quantum ESPRESSO suite, 11 using SG15 (PBE) pseudopotentials. 5The code expands wave functions ψ K with crystal momentum K as where ψK is the plane wave representation of the wave function and the sum runs over the reciprocal lattice vectors G.The summation is cut off at values of G corresponding to a kinetic energy of 40 Ry.For MoS 2 , we use a 24 × 24 k-grid and, in case of hybrid functionals, an 80 Ry cutoff for the Fock exchange, as well as a 12 × 12 q-grid.In-plane and out-of-plane lattice constants of 3.18 Å and 15 Å are used, respectively, the former being obtained from geometry optimization, the latter being chosen to include a sufficient amount of vacuum between replicas in the out-of-plane direction.
Eq. (S1) greatly facilitates the calculation of wave packets and unfolded band structures.
The molecule is placed in a simulation cell spanned by lattice vectors A i that are constructed as superpositions of the primitive lattice vectors a j of MoS 2 , where T is a 3 × 3 transformation matrix with integer coefficients.Bloch functions u e iK•r with A i -periodic cell functions u are characterized by a Bloch vector K in the small BZ associated with the A i lattice.If this K is connected to a Bloch vector k within the large BZ related to the a j lattice by a reciprocal lattice vector G of the the square modulus of the corresponding wave packet coefficient in Eq. ( 1) of the main text can be computed as where the sum runs over all reciprocal lattice vectors g of the a j lattice.The unfolded band structure is obtained from the |c nk | 2 of different states ψ nK as where ε nK is the Kohn-Sham eigenvalue of ψ nK .We perform this calculation in postprocessing using an in-house developed code 12 to obtain both the one-dimensional (band structure) and two-dimensional depictions presented in Fig. 3 of the main text.For the two-dimensional representations, the molecule was rotated starting at 15 degrees in steps of 30 degrees, averaging the corresponding spectral functions.For the quasi-degenerate lowest unoccupied orbital, both partner orbitals are averaged over.

S3
To determine the dielectric constants of MoS 2 , we employ the random-phase approximation (RPA) on the basis of a PBE electronic structure, obtained with Quantum Espresso and SG15 pseudopotentials with the same k-mesh and basis set cutoffs as above.The linear-response calculation is conducted with the Yambo code, 13 including 300 unoccupied bands.The dielectric constants associated with the three-dimensional cell containing the MoS 2 monolayer are transformed into intrinsic material values as described in the following section, yielding ε ∥ = 16.52,ε ⊥ = 10.08, and a thickness of 5.46 Å. Substrates are modeled as isotropic polarizable media with dielectric constant ε s , separated from the MoS 2 layer by a vacuum spacer of size 0.69 Å. 14 The G 0 W 0 calculations are performed with Yambo, again using the PBE electronic structure as a starting point, but increasing the out-of-plane lattice constant to 20 Å.A 150-point full-frequency integration is employed to determine the correlation part of the self-energy, using an 80 eV cutoff in the calculation of the RPA-level screening.We repeat the simulation with 100 and 200 unoccupied bands and extrapolate the resulting self-energy correction to infinite bands. 15We account for the anisotropy in the q → 0 component of the screening and use a two-dimensional slab Coulomb cutoff. 16sualization of geometries and orbitals was carried out with the XCrysDen program. 17

Dielectric Models
The dielectric models for the different environments considered in this work as well as the values of the parameters are shown in Fig. 2c).The parameters are determined from first principles as described in the following.

Determination of the Dielectric Constant of a Molecular Film
The dielectric constant ε of a molecular film is determined from DFT in conjunction with the isotropic PCM, using a formula that connects the microscopic polarizability with the macroscopic dielectric constant.Taking a starting guess for ε, we calculate with Octopus (PBE functional) the dipole moment ∆p(ε) = (∆p x (ε), ∆p y (ε), ∆p z (ε)) induced by a static electric field E(1, 1, 1) and the corresponding ε-dependent local cavity field, 18 then update where v is the volume of the PCM cavity, computed with a Monte-Carlo integration.We repeat this procedure until self-consistency, obtaining a value of ε = 5.32 for an H 2 Pc film.

Determination of the Dielectric Constants of a Two-Dimensional Material
The dielectric constants ε ∥ and ε ⊥ of a polarizable slab replacing a MoS 2 monolayer are determined using the approach outlined in Ref. 19.In brief, we perform linear-response calculations (RPA@PBE) with a simulation cell containing MoS 2 layers separated by a large amount of vacuum.The resulting dielectric constants ε∥ and ε⊥ represent an average polarization of the slab and the vacuum spacer.To recover the MoS 2 -intrinsic values, we transform ε∥ and ε⊥ according to which can be derived from electrostatic boundary conditions. 20In these equations, c is the height of the simulation cell and d is the thickness of the slab, which is adjusted until the smoothly transitions into the plane-averaged PBE exchange-correlation potential.This approach corresponds to the explicit construction of a smooth and consistent exchangecorrelation self-energy. 21,22We end up with the values ε ∥ = 16.52,ε ⊥ = 10.08, and d = 5.46 Å.

Bandgap Renormalization For MoS 2 : The ∆W model
The renormalization of the band gap of MoS 2 is estimated with an electrostatic model.
We consider the setup in Fig. S1a), where a charge q is added at the position r within the substrate-sustained MoS 2 monolayer, parametrized as specified in Fig. 2c) in the main text.As shown in Ref. 14, the corresponding Poisson's equation can be solved without further approximations to determine the full electrostatic potential φ, including polarization contributions.The potential φ(r ′ ) felt by a test charge q ′ at r ′ is proportional to the statically screened interaction W between q and q ′ , W (r, r ′ ) = φ(r ′ )/q.Combining this result with the screened interaction W 0 for a reference setup without substrate [Fig.S1b)], we obtain a self-energy correction to the GW band structure of an isolated monolayer due to substrate where "+" and "-" refer to occupied and unoccupied bands, respectively, and we set r to the center of the monolayer.W 0 has to be subtracted as a corresponding contribution is already accounted for in the GW calculation of the isolated monolayer.

Figure S1 :
Figure S1: Polarization charges emerging as a consequence of a charge q added in the MoS 2 layer

Figure
Figure S2: a) Wave packet representation of the highest occupied orbital of H 2 Pc in the MoS 2 Brillouin zone upon rotation around the principal axis.The arrows indicate the folding processes that map features back into the first BZ, based on comparison to the unfolded momentum map in Fig. 8 of Ref. 34. b) Evolution of the azimuthally averaged wave packet upon out-of-plane rotation.