InterPhon: Ab initio Interface Phonon Calculations within a 3D Electronic Structure Framework

This work provides the community with an easily executable open-source Python package designed to automize the evaluation of Interfacial Phonons (InterPhon). Its strategy of arbitrarily defining the interfacial region and periodicity alleviates the excessive computational cost in applying ab initio phonon calculations to interfaces and enables efficient extraction of interfacial phonons. InterPhon makes it possible to apply all of the phonon-based predictions that have been available for bulk systems, to interfacial systems. The first example, in which this package was applied to InAs surfaces, demonstrates a systematic structure search for unexplored surface reconstructions, navigated by the imaginary mode of surface phonons. It eventually explains the anisotropic surface vibrations of the polar crystal. The second example, involving oxygen adsorption on Cu, reveals adsorption-induced vibrational change and its contribution to energetic stability. The third example, on a Si/GaAs interface, shows distinct vibrational patterns depending on interfacial structures. It leads to a prediction regarding the structural transition of interfaces and unveils the processing conditions for spontaneous growth of GaAs nanowires on Si. High-level automation in InterPhon will be of great help in elucidating interfacial atomic dynamics and in implementing an automated computational workflow for diverse interfacial systems.


Introduction
The computation of collective atomic vibrations, called phonons, is a key to predicting a variety of material properties: transport properties such as diffusivity and conductivity; thermal properties such as heat capacity and entropy; and structural properties such as equilibrium phase and transition pathways [1]. With the aid of tremendous advancements in both hardware (highperformance parallel computers) and software (density functional theory (DFT) codes), a large amount of the energy and force state of arbitrarily given materials, E({RI}) and F({RI}), where {RI} indicates the positions of constituent atoms, can now be concurrently computed using available codes of electronic structure calculations. In addition to the quantum mechanical derivation of energy and force, the advent of several phonon codes, including Phonopy [1], PHON [2], and others [3], has made ab initio phonon calculations routine subjects of research. Nowadays, ab initio phonon calculations are even incorporated as a step in sequences of automatic multistep workflows, eliminating tedious user intervention and manual management during each calculation step [4]. By combining such an automatic platform with a database or crystal structure prediction (CSP) algorithm, massive phonon calculations can be carried out within high-throughput screening approaches to discover unknown materials with optimal properties. In the entire workflow, the phonon calculations usually play a vital role in assessing structural stability [5,6].
However, all of the phonon computations mentioned above are restricted to bulk systems. This is because a great deal of computational cost (approximately one or two orders higher than that of conventional electronic relaxation) is required to execute the ab initio phonon calculations. The required number of computations is proportional to the number of atoms per unit cell, so lower cost is expected for smaller bulk systems. A large number of symmetry operations of bulk crystals also provide an opportunity to substantially reduce the required number of computations [1,2]. In contrast, the symmetrically reduced and thus viable computational cost does not apply to defective systems with broken symmetry, such as surface and interfacial systems. A defective region introduces different kinds of vibrations originating from the different stoichiometry and bonding geometry compared to the bulk environment [7,8]. When performing defect calculations within three-dimensional (3D) DFT codes, however, the supercells of slab geometry or those with confined point defects must be generated [9,10] in order to remove any fictitious interactions between periodic cells [11,12]. Accordingly, a large number of atoms and a low degree of symmetry are intrinsically inevitable, which impose enormous computational cost. This cost is normally not affordable; thus, the effects of localized vibrations are commonly neglected in ab initio calculations for defective systems [13][14][15][16]. The accuracy of the results relies on the assumption that the vibrational contributions to the total free energy are not critical and only the electronic contributions are sufficient, which may not necessarily be the case.
Recently, more and more reports have emphasized the importance of the vibrations in interfacial (including surface) regions, for example, the atomistic mechanisms of a temperature-driven change in metal-insulator surface structure [17][18][19], equilibrium crystal shape [20], and nanowire growth behavior [21,22]. All of these, which were not explicable using only electronic energy contributions, were theoretically elucidated by looking at the subtle interplay between electronic and vibrational energies. Moreover, the emergent phenomena at interfaces [23][24][25][26], such as the interfacial phase-change memory (iPCM) [27] and multiferroics [28,29], are highly related to the intriguing inter-coupling between electronic and atomistic (vibrational) structures within the localized interfacial region. Therefore, the ability to compute interfacial phonon behavior will pave the way for a great deal of further research, and the development of a cost-effective and easily executable phonon code is timely.

In this work, an open-source Python package for interfacial phonon calculations, called
InterPhon, is reported. It characterizes the distinct lattice vibrations with reduced periodicity (onedimensional (1D) or two-dimensional (2D)) by automatically processing the information obtained by 3D-based DFT codes. It fully implements a set of tools to promote efficient phonon calculations and analyses for arbitrary interfacial systems in reduced symmetry and periodicity. Individuals with access to existing DFT programs can easily utilize the InterPhon package with the help of extensible compatibility supported by its adapter modules. The only requirements for package execution are a DFT input file representing the structure of interest and the corresponding output files providing the forces acting on each atom.
To demonstrate the usefulness of InterPhon and its applicability to a variety of material systems, three examples of package usage are presented. The examples investigate changes in interfacial phonon behavior that are associated with interfacial bond formation. They are surface reconstruction of InAs(111)A and (111)B, adsorption of oxygen on Cu(111), and interface of (111)-oriented Si/GaAs. The investigations demonstrate that InterPhon makes it possible to apply all phonon-based investigations that have been available for bulk systems, such as structure search, assessment of structural stability, and temperature-driven phase transition, to interfacial systems.
Such analyses have hardly been achieved due to the scarcity of a cost-effective and easily executable software for interfacial phonons. Since the reconstruction, adsorption, and interface calculations are the general basis of ab initio modeling of material interactions, the InterPhon code will be valuable for designing diverse interactions by enabling predictions of interfacial dynamics at the atomic scale.

Phonon computation
A phonon calculation starts with transforming a displacement vector field defined over lattice domain into a numerical problem of eigenvalues (phonon frequencies) and eigenvectors (phonon modes) of the dynamical matrix [30]. There are essentially two ways to evaluate the dynamical matrix: the direct approach using the finite displacement method (FDM) [31,32] and the linearresponse approach using density functional perturbation theory (DFPT) [33,34]. InterPhon is based on the direct approach within the FDM scheme; the specific formalism used to evaluate the interatomic force constants (Equation S4) and dynamical matrix (Equation S6) is summarized in Section S1 of Supporting Information. Note that the term 'unit cell' is used to refer to the cell of interest for phonon evaluation (e.g., phonons of interface unit cell) throughout this work. When obtaining the force constant elements, atomic forces in response to a set of atomic displacements are required and the forces should be evaluated using supercell consisting of several unit cells. This is because the DFT forces calculated via the Hellmann-Feynman theorem are actually sums over forces between one atom and all periodic images of another [31]. Assuming that atomic interactions for separation larger than supercell size (which eventually determines cutoff distance) are zero, the periodic contribution is neglected and the given DFT forces can be treated as the forces induced by the two-atom interaction. Then, considering the original periodicity, the force constant matrix of supercells is merged back into the unit cell dynamical matrix (i.e., lattice-sum over l' in Equation S6), which encompasses long-range atomic interactions.

InterPhon strategy: selection of the interfacial atoms
The direct approach within the FDM is also the operating mechanism behind widely used phonon codes for 3D bulk systems, such as Phonopy [1] and PHON [2]. These conventional codes, which take advantage of the symmetry in 3D space groups, generate the dynamical matrix of all constituent atoms in a unit cell by requiring 'complete FDM' for all atoms. By contrast, to make the interfacial phonon calculation feasible, InterPhon, which requires 'selective FDM', generates the dynamical matrix of a subset of the constituent atoms in a unit cell and does not assume 3D periodicity. In conventional ab initio interface calculations, the excessive number of atoms is inevitable to represent interfacial region in a 3D periodic cell. Among all of these atoms, however, the bonding environment of most of the atoms located far from the interface (gray atoms in Fig.   1a) is similar to that of the bulk atoms, presenting the same vibrations as the ideal bulk (as intensively tested in Section S4.2 of Supporting Information). Deviation of vibrations from bulk occurs only on the atoms in the vicinity of the interface (green and orange atoms in Fig. 1a).
Therefore, InterPhon focuses on the interfacial atoms by allowing users to easily select atoms to be considered as the interface. In addition, contrary to the compulsory setting limited to the 3D periodic boundary condition (PBC) in existing DFT and phonon codes, user-defined setting of the PBC is allowed in InterPhon. As a result, phonon evaluation proceeds only in the interfacial region with reduced PBC dimensions.
The key strategy behind InterPhon, the user-defined setting of the interfacial region and PBC, reduces the cost of the most time-consuming element of phonon computations (that is, DFT force calculations for displaced supercells) by approximately one order of magnitude. This makes the calculation loads affordable. Eventually, the operation architecture, which is fundamentally different from that of the phonon codes for 3D bulk systems, allows efficient extraction of the localized phonons at arbitrary reciprocal points (called k-points) in 1D or 2D Brillouin zone (BZ).
In addition, focusing on one of two exposed interfaces (top and bottom in slab geometry) makes it  operation. The program execution starts with reading DFT input/output files in the 'inout' subpackage, followed by the numerical calculations in the 'core' sub-package, and ends with printing out the phonon properties by the 'analysis' sub-package. In the middle of the processes, the functionality supported in the 'util' and 'error' sub-packages are internally employed.

Symmetry functionality
The 2D symmetry of interfaces, although its number of operations is usually smaller than that in 3D crystal, enables further reduction of the number of DFT calculations required. In the following explanation on translational and point group symmetry, the notations summarized in Appendices B and C are used to avoid confusion, possibly caused by various representations of vector, matrix, and their product. By translational symmetry, translation by in-plane lattice vectors leaves the interface system unchanged. Therefore, the atomic displacements are only required to the atoms in one unit cell; total number of required displacements corresponds to the number of interfacial atoms in unit cell multiplied by 6, where 6 comes from the forward and backward displacements for each of the three Cartesian directions (see Equation S4). Note that the conventional translational invariance in the force constant matrix [2,31]: which is derived from the fact that the atomic forces must be zero when a crystal shifts rigidly (i.e., = 0 when ′ ′ is constant for all the atoms in Equation B6), is usually not satisfied in interfacial phonon calculations with InterPhon. This is because unlike 2D materials such as graphene, the freely translational motions of interfaces are limited by the presence of bonds with neighboring bulk. Accordingly, acoustic phonon modes with zero frequencies at the -point, corresponding to translational motions, will not be found in interfacial phonons as shown in the calculation results of this study.
In addition, the use of point group symmetry can reduce the number of atomic displacements even further. The displacement vectors do not necessarily have to be the exact Cartesian directions, but just needs to be three linearly independent vectors. If ′ ′ is set and the other two linearly independent displacement vectors, ′ ′ and ′ ′ of the atom at ′ ′ , can be obtained by applying point group symmetry operations { } of the interface plane group, displacement only along the one direction is sufficient for this atom [2,35]. After creating a set of three linearly independent displacement ( ′ ′ , ′ ′ , ′ ′ ) and resulting force vectors ( , ′ ′ , , ′ ′ , , ′ ′ ) for an two-atom pair ( and ′ ′ ), the force constant matrix in orthogonal Cartesian coordinates (Equation B4) can be recovered by: Furthermore, if a point group symmetry operation exists by which another two-atom pair (̃ and ′ ′ ) is mapped onto, the force constant matrix for the image atom pair does not need to be evaluated directly and is simply generated by: which is derived in Appendix D.
In the backend of InterPhon, this translational and point group symmetry mapping is conducted by automatically searching the in-plane symmetry operations of a given interface structure.  [36]. This symmetry functionality is implemented as default in InterPhon and verified in Fig. S1 showing that applying symmetry operations enables a smaller amount of DFT calculations to provide the same phonon dispersion as the one from larger calculations without symmetry usage.

Architecture of InterPhon
InterPhon, which is steered by high-level user interfaces, automatically performs all of the processes required to characterize interfacial phonons from given DFT input and output files. The overall InterPhon package consists of five sub-packages: (i) 'inout', (ii) 'core', (iii) 'analysis', (iv) 'util', and (v) 'error'. Each of the sub-packages takes different roles and the main workflow is as follows: (i) modules in the 'inout' sub-package read and parser DFT input/output files; (ii) this information required to generate the dynamical matrix, such as crystal structure and atomic forces, is provided to modules in the 'core' sub-package that calculate phonon frequencies and modes; (iii) the processed information is provided to modules in the 'analysis' sub-package and the internal processes end with printing of the phonon properties, such as the density of states (DOS), band, vibrational motions (phonon mode), and thermal properties (vibrational entropy and free energy) in the complete forms of both data files and graphics (Fig. 1b). In the graphical representations, inspired by PyProcar [37], atom projection in the phonon band and DOS is effectively visualized in high quality to help analysis of phonon characteristics, using matplotlib [38] as a graphic library. In the middle of the processes, the modules in (iv) 'util' and (v) 'error' sub-packages are internally employed to enrich functionality. The detailed design architecture of the package and interrelation among the sub-modules are shown in Fig. S2, along with a full explanation of the operating procedures, in Section S3 of Supporting Information. The usage of parser modules in the 'inout' sub-package, which communicate with external files, allows a standardized representation of the internal data, irrespective of external DFT codes having different file formats. The parsers for VASP [39][40][41], Quantum ESPRESSO [42], and FHI-aims [43] are currently supported, and the parsers for other codes are under development.

Parameters of InterPhon
The key user-parameters of InterPhon are displacement (default = 0.02), enlargement (default = "2 2 1"), and periodicity (default = "1 1 0"). The default values mean that the displacement length for FDM (i.e., u in Equation S4) is 0.02 Å ; the extension ratios along the a1, a2, and a3 lattice directions are 2, 2, and 1, which determine the size of supercells and accordingly define the cutoff range of the considered interatomic forces; and the periodic (1 or True) boundary conditions are along the a1 and a2 directions and open (0 or False) along the a3 direction. According to the periodicity parameter, the interface plane direction, to which the in-plane symmetry functionality is carried out, is determined and even 0D (like isolated clusters) and 1D PBC systems can be handled. In addition, InterPhon defines the interfacial region according to the automatic reading of the statement of constraints on atom movements, which is supposed to be set by users within the DFT structure file (see Section S4.1 of Supporting Information). By providing users with the opportunity to define the interfacial region, the convergence of the phonon characteristics with respect to the increase in interfacial region thickness (i.e., the number of selected atom layers) can easily be confirmed. Since the explicit boundary between the interfacial region and the bulk cannot be universally defined, this kind of convergence test is essential. An effective methodology to perform the convergence test is suggested in Section S4.2 of Supporting Information, and the convergence of all results in this work was tested. Except for the InAs surface reconstructions, where the bonding states change significantly, fast convergence in vibrations with respect to interfacial thickness was confirmed. This is a very encouraging result because it implies that further reduction of the interfacial region than was done in this work, where three layers were selected as the interfacial region, is possible, which would further decrease the computational cost.

Data and library availability
The following section presents the vibrational change introduced by interfacial bond formation, which was investigated by InterPhon in conjunction with VASP code. All of the data used in these investigations are available in the 'example' folder at the source code repository. The data provenance showing how the data were obtained can be traced back, enabling the reader to reproduce the results. The latest version of the InterPhon package has been released in the Github repository at https://github.com/inwonyeu/interphon, which can be easily installed via pip:

$ pip install interphon
In addition, a complete user manual on the installation, execution, file structure, and operation options is clearly described in the documentation at https://interphon.readthedocs.io. Following this documentation, InterPhon can be executed in two ways: command line and Python script (or interpreter). While the execution by command line is easy to follow, the Python script is more appropriate for dynamic utilization according to the user-specific purpose. Ensuring code reuse and modularity of the Python language, the conjunction of compiled DFT codes with this objectoriented package will be an ideal approach to automated computation projects for heterogeneous interfacial systems [44,45].

Application examples of InterPhon and discussion
The

Surface reconstruction: InAs(111)A and (111)B
The crystal structure of InAs is zinc-blende (ZB, F 4 ̅ 3m ), which is non-centrosymmetric.  (selected atoms), is shown in Fig. 2c. The dispersion was obtained by setting one of the parameters of the 'plot' method in 'Band' class as option = 'projection' (for all of the available plotting options, see the online user manual). An imaginary mode is found at points ̅ and ̅ , presumably due to the short cutoff range of interatomic forces [31] under the assumption that the force between atoms separated by more than one unit cell (whose cell size is a1 = a2 = 8.5 Å ) is zero. Therefore, the interatomic forces were recalculated to incorporate the force interactions outside of the unit    Fig. 3a shows that the bottom In-atoms of the (111)A surface were passivated by pseudo-hydrogen with 1.25 valence electrons. Figure 3b shows the top view of the (111)B surface and the process of creating the Invacancy β(2×2) reconstruction, which was suggested to be stable under In-rich conditions [20].
After setting one vacant site on the topmost As-layer per (2×2) cell and substituting the remaining three As-atoms with In-atoms, a relaxation calculation was performed for the top five layers, resulting in the In-vacancy α(2×2) reconstruction. For the relaxed unit cell, further phonon analysis was carried out for the top three layers (a total of 11 atoms consisting of 3, 4, and 4 atoms in the first, second, and third layer, respectively) by setting the user-parameters of InterPhon as follows: displacement = 0.02, enlargement = "2 2 1", and periodicity = "1 1 0". The corresponding phonon dispersion projected to the first layer (i.e., three In-atoms) for 33 branches is shown in Fig.   3c. Although the cutoff range of the interatomic forces was set to be the twice the length of the reconstruction unit cell (dotted parallelogram in Fig. 3b), which is large enough to incorporate long-range interactions, imaginary modes are found at all of the high-symmetry points, ̅ , ̅ , and ̅ . This implies that the In-vacancy α(2×2) reconstruction is dynamically unstable. The collective atomic motion corresponding to the lowest imaginary frequency at point ̅ (indicated by a red circle and arrow in Fig. 3c) is indicated by the blue arrows on the topmost bilayer in Fig. 3b. Since the imaginary frequency (conventionally represented as negative real numbers) means that there is no restoring force against atomic motion along the corresponding imaginary mode, another relaxation calculation was conducted after slightly deforming the atomic positions along this imaginary mode.
One possible reason why the first relaxation converges to the In-vacancy α(2×2) is the constraint in the DFT energy minimization algorithm, which is caused by the energy landscape and initial symmetry. By escaping from the energy saddle point of the In-vacancy α(2×2) and changing to the In-vacancy β(2×2), the DFT electronic energy was lowered by 0.32 eV, corresponding to a decrease in surface energy by 5.1 meV/Å 2 . This finding demonstrates that such an analysis of a surface-localized phonon provides a great opportunity to systematically search for stable surface structures, which is analogous to a systematic search of the stable bulk phase by 3D ab initio phonon calculations [48].
As the reconstruction changes, the surface phonon dispersion also changes. Figure 3d shows the first-layer-projected band along with the layer-projected DOSs of the (111)B In-vacancy β(2×2).
Although all DFT and InterPhon calculation parameters were the same as those used for the Invacancy α(2×2) in Fig. 3c, the imaginary modes disappear, showing the dynamical stability of the In-vacancy β(2×2). Furthermore, the phonon states contributed by the In-atoms in the first topmost layer are more localized to the high-frequency states in the In-vacancy β(2×2) compared to Invacancy α(2×2), which is clearly shown by the color code presented by "contribution of 1st layer".
On the other hand, the absence of a third layer contribution to the high-frequency vibrations (clearly shown in the layer-projected DOSs) implies that the atomic interactions between the Inatoms in the topmost bilayer and the As-atoms in the third layer become negligible in this reconstruction. This change in bonding states results in weaker bonding and substantially lower frequency than bulk InAs (ZB, dotted line in Fig. 3d).
Finally, the different stoichiometry and bonding geometry between the (111)A In-vacancy(2×2) and (111)B In-vacancy β(2×2) lead to substantial differences in surface phonons, as shown in Fig.   2d and Fig. 3d. Therefore, the corresponding vibrational free energy ( ), which can be evaluated containing the electronic and vibrational contributions, can be predicted as a function of temperature and vapor pressure by a method called ab initio thermodynamics that considers the contact-induced interactions between vapor environment and solid; the detailed method is fully described in method section of the authors' previous papers [21,22]. Comparing the (111)A and (111)B reconstructions, the different slope of is caused by differences in implicit temperature dependence, which comes from the atomic stoichiometry mismatch between surface and bulk [8]. Furthermore, it should be noted that the reduction in γ by the contribution of , called explicit temperature dependence, is much more significant in the (111)B reconstruction because of its substantially lower frequency compared to bulk InAs and the (111)A reconstruction ( Fig. 2d and Fig. 3d). Since a majority of modern functional materials are made of polar crystals (like zinc-blende structures) [3], the ability to elucidate this kind of anisotropic surface vibration will be valuable for other scientifically and technologically important problems such as ferroelectrics and multiferroics.

Adsorption: Oxygen on Cu(111)
The adsorption energy on a metal surface is widely used as an effective descriptor to predict catalytic activity. Based on a database of calculated adsorption energies, new catalysts with improved selectivity for complex surface reactions, such as the oxygen reduction reaction (ORR) and oxygen evolution reaction (OER), have been explored among a variety of metal species, alloy compositions, surface orientation, and so on [50,51]. In this section, the difference in vibration before and after the adsorption of oxygen on Cu(111) and its effects on energetic stability were examined to demonstrate how InterPhon can be applied to promote accurate adsorption-related calculations. This will contribute to various theoretical studies, including studies of deposition and growth from vapor as well as catalytic behavior.
The crystal structure of Cu is face-centered cubic (FCC, Fm3 ̅ m) which is centrosymmetric.
Therefore, contrary to the ZB structure, the two opposite surfaces of (111) are equivalent; hence, a symmetric slab structure without hydrogen passivation was used. Figure 4a shows a side view of the slab along the <111> direction, while the corresponding top view is shown in Fig. 4b. For oxygen adsorption, two sites of Cu(111), FCC and hexagonal close-packed (HCP), were considered using an adsorption density of 0.25 ML, as shown in Fig. 4b. These phonon DOSs were integrated to evaluate vibrational free energy ( ) [1,2,49], which is a constituent term in the change of Gibbs free energy by adsorption: where • and are the Gibbs free energy of the adsorbed surface and pristine surface, respectively; ∆ is the change in electronic energy by adsorption; • and are the vibrational free energy of the adsorbed surface and pristine surface, respectively; and ( 2 ) is the chemical potential of oxygen in vapor phase, which is available in thermochemical tables or can be calculated statistically under the ideal gas assumption [7]. Note that most previous DFT calculations on adsorption considered only the electronic adsorption energy, ∆ .
Recently, there have been attempts to evaluate the Gibbs free energy of adsorption within the following approximation [52]: the vibrations of surface region is not much affected by adsorption, so the free energy contribution from surface vibrations is negligible, Under this approximation, only contributions from vapor species were considered, and the vibrational free energy of an adsorbate (e.g., ( )• ) was calculated using the discrete vibrational frequencies (e.g., total three frequencies in the case of a monoatomic adsorbate) at point ̅ [11,53]. Ignoring the vibrational contribution of the surface may be valid if the density of adsorbates is low. As the coverage and bonds between the adsorbate and surface increase, however, the contribution will no longer be negligible. Figure 4g  increases and becomes comparable to the electronic energy difference. This is consistent with the following tendency observed in several bulk and surface systems: a structure with higher electronic energy (weaker bonding) shows lower vibrational energy, leading to a temperature-and entropydriven phase transition [8,17,18,20,54].

Interface: Si(111)/GaAs(111)
The interfacial region between two different bulk phases is one of the most diverse and interesting material systems where a variety of interactions takes place. However, in numerous cases, structural models proposed by simple electronic energy calculations have failed to explain certain experimental observations, such as intermixing between GeTe/Sb2Te3 superlattices relevant to the resistance change of iPCM [55] and interface stoichiometry at epitaxial antiperovskite/perovskite heterostructures [56]. If an interfacial structure can successfully be described at the atomic scale by considering the contribution of vibrational energy, the interface control will be greatly facilitated, paving the way for unprecedented material design. In this section, using (111)-oriented Si/GaAs as an example, various types of interfacial structures and temperature-driven changes of those structures are investigated to demonstrate how InterPhon can be applied to predict interfacial structural transitions.
The crystal structure of Si is diamond (Fd3 ̅ m), which is centrosymmetric, while that of GaAs is ZB, which is non-centrosymmetric. As explained above with InAs, GaAs shows an intrinsic polarity along the <111> directions, exposing two inequivalent surfaces: Ga-terminated or A-polar (111)A surface; and the opposite As-terminated or B-polar (111)B surface. In contrast, the Si (111) surface is non-polar; thus, vertical growth of GaAs on a Si(111) substrate can be either A-polar or B-polar. Figure 5a shows a side view of the (111)-oriented Si/GaAs interfaces where three bonds between Si and Ga (or As) form. Throughout this section, this kind of interface will be referred to as type-I. On the other hand, at a type-II interface, only one bond between Si and Ga (or As) forms, as shown in Fig. 5b.   The phonon DOSs of the interfaces were integrated to evaluate the phonon contributions to the free energy of interfaces [1,2,49]. Then, the vibrational free energy difference between interface and bulk, called the vibrational interface energy ( ), was calculated: where ,  The calculation method of how to obtain the temperature-pressure dependent interface energy are the same as the method for surface energy introduced with Fig. 3e [21,22]. Note that the relative electronic interface energy within the same polarity (= − / − − − / − ) is the best that could possibly be obtained, since the absolute energy of a strained polar surface or interface cannot be calculated using the current levels of DFT simulation [57,58]. By definition, the positive relative interface energy means that the type-I interface is more stable than the type-II interface.
The results evaluated only by electronic contribution are represented by solid lines, while those reflecting the sum of electronic and vibrational contribution are shown with dotted lines. In particular, Fig. 6c demonstrates that the type-I interface is more stable than type-II only for the Bpolar interface. In addition, the stable temperature ranges of the I-B-polar increase from ~670 to ~740 K by considering the vibrational contribution. Finally, the temperature-pressure boundary between the stability region of the I-B-polar and II-B-polar interfaces are shown in Fig. 6d Followed by a gradual temperature increase to growth temperature (1023 K) and Ga-source injection, the surface pretreatment eventually facilitates the growth of vertically aligned GaAs NWs on a Si(111) substrate, achieving 100% yield [59], since GaAs NWs preferentially grow along the <111>B rather than <111>A direction [47]. According to our prediction, Si-As bonding at a type-I interface leads to vertical B-polar directionality (bottom panel of Fig. 5a), while Si-As bonding at a type-II interface leads to vertical A-polar directionality (top panel of Fig. 5b).
Therefore, taking into account the vibrational contribution, the systematic calculations on the Si/GaAs interface clearly explain why pretreatment of the Si(111) substrate by an As-source had to be carried out under low-temperature conditions to promote vertically aligned growth of B-polar GaAs NWs.

Conclusions
To facilitate theoretical studies on interfacial phonons, a Python package (InterPhon) connected to the 3D framework of available DFT codes was developed, which supports a set of methods to efficiently extract localized lattice vibrations with high-level user interfaces. The InterPhon strategy of user-defined interfacial region and user-defined periodicity was validated showing that the phonons of atom layers far from the interface gradually converge to the bulk phonons.
Furthermore, the instructions of design architecture and package workflow were provided so as to be ready for a potential adaptation to unexpected problems and other computational projects. The

Author contributions
In

Conflicts of interest
There are no conflicts to declare.

Appendix B. Vector and matrix notations
The coordinates of an atom are defined by three Cartesian components: where l and s labels indicate the lattice and basis, respectively.
Displacement vector of the atom located at along an arbitrary direction is represented by a column vector: which describes the proportionality between displacements and forces: The whole force constant matrix becomes 3N×3M matrix where N and M are the number of interfacial atoms in supercell and unit cell, respectively.
A component of the total force vector acting on atom is evaluated by sums over forces contributed by all of the atomic displacements: where the labels of α and β indicate one of the three Cartesian components. This is the same expression as Equation S3 derived from the Newtonian equation of motion.

Appendix C. Symmetry operation
An affine mapping ( , ) transforms into an image point: where is a rotation part represented by 3×3 matrix and is a translation part represented by 3×1 column vector: ).