Atomic and electronic structure of Lomer dislocations at CdTe bicrystal interface

Extended defects are of considerable importance in determining the electronic properties of semiconductors, especially in photovoltaics (PVs), due to their effects on electron-hole recombination. We employ model systems to study the effects of dislocations in CdTe by constructing grain boundaries using wafer bonding. Atomic-resolution scanning transmission electron microscopy (STEM) of a [1–10]/(110) 4.8° tilt grain boundary reveals that the interface is composed of three distinct types of Lomer dislocations. Geometrical phase analysis is used to map strain fields, while STEM and density functional theory (DFT) modeling determine the atomic structure at the interface. The electronic structure of the dislocation cores calculated using DFT shows significant mid-gap states and different charge-channeling tendencies. Cl-doping is shown to reduce the midgap states, while maintaining the charge separation effects. This report offers novel avenues for exploring grain boundary effects in CdTe-based solar cells by fabricating controlled bicrystal interfaces and systematic atomic-scale analysis.

To date, theoretical and experimental studies of CdTe GBs have covered only tilt Σ 3 (111), and Σ 9 (112) interfaces, as found within polycrystalline CdTe crystals 8,[17][18][19][20][21][22] . Such an approach of studying GB atomic structures in poly-crystalline thin films is dependent upon luck of being able to find two grains with exactly aligned and viewable zone axis for STEM imaging. A better approach would be to fabricate specific grain boundaries with well-defined misorientation angles using two bulk CdTe wafers. Here, we have developed a method for bonding two single crystal samples together into one bicrystal, the so called wafer bonding method, which is used to create artificial GBs with well-defined misorientation angles, and interfacial planes, thereby allowing us to explore a large portion of the grain boundaries' configuration space [23][24][25][26] . The wafer bonding technique has been used to investigate interfaces in MgO, Bi 2 Te, and Al 2 O 3 11,[27][28][29] .
In this paper, we have combined high-angle annular dark-field (HAADF) imaging in an AC-STEM, first-principles density functional theory (DFT) modeling, and geometrical phase analysis (GPA), to study a [1][2][3][4][5][6][7][8][9][10]/(110) 4.8° tilt GB created by wafer boding. We find that the grain boundary is composed of several different types of Lomer dislocations, and determine the atomic structures of the dislocation cores, which provide evidence of the reaction mechanism of these dislocations. The electronic structure, computed using first-principles DFT on realistic, STEM-matched, models containing up to 1917 atoms, shows significant creation of mid-gap states and distinct charge-channeling tendencies at the dislocation cores. Using the precise atomistic model obtained from STEM imaging on the created artificial grain boundary, the effects of Cl passivants on the mid-gap states and charge-channeling effects are shown with DFT calculations, demonstrating the effectiveness of the current approach towards understanding and improving grain boundaries in semiconductor electronics.

Validation of dislocation array in bicrystal.
A [1][2][3][4][5][6][7][8][9][10]/(110) 4.8° CdTe tilt GB was fabricated by wafer bonding. The symmetrical θ = 4.8° tilt angle was achieved by cutting two single crystal CdTe crystals by 2.4° from the (110) plane as shown in Fig. 1a,b. Selected area diffraction pattern (SAPD) taken from the bonded interface confirm the misorientation angle of the grain boundary using the splitting of the Bragg spots shown in Fig. 1c. Figure 1d shows an atomic-resolution HAADF image of the [1-10]/(110) 4.8° CdTe tilt GB. In order to accommodate a tilt angle θ, the interface is expected to contain a periodic array of evenly spaced dislocations having on average |b| edge components. The dislocation spacing, D, is determined by Frank's formula D = |b|/θ, where |b| is the magnitude of the dislocation's Burgers vector. Using atomic-resolution HAADF images, we find that the average spacing between two adjacent dislocations is 5.6 ± 1.0 nm. This is consistent with D = 5.5 nm, as expected from Frank's formula for a 4.8° tilt grain boundary with a dislocation's Burgers vector of b = (a/2)[110]. Such a Burgers vector indicates the pure edge character of the dislocations at the grain boundary. In addition, we did not find any significant interfacial faceting or interfacial amorphous layers, demonstrating that our bicrystal fabrication has been successful.
Multiple configurations of Lomer dislocations. The great majority of dislocations found along the grain boundary have a Burgers vector b = (a/2)[110]. This type of dislocations can also be found in materials with diamond cubic crystal structure (Si, diamond), and is called a Lomer dislocation. These pure-edge dislocations are immobile since their slip plane would take place on {001} planes, which do not belong to the usual < 110 > {111} slip system of CdTe. Three distinct types of Lomer dislocations can be found in our HAADF-STEM images, which we will call Type I, II and III throughout the remainder of this discussion. Burgers circuits around the different types of dislocations are shown in Figs 2a,b and 3a, respectively. The DFT-relaxed atomic structures of Type I, II and III dislocations are shown overlaid on experimental HAADF-STEM images in Figs 2c,d and 3d, respectively, demonstrating that the DFT-relaxed atomic structures match elegantly the STEM images. The atomistic model construction procedures and DFT calculations are described in the Methods Section and Supplementary Fig. 1. Although all the Lomer dislocation types have the same total Burgers vector, b, the atomic structures are quite different. In the Type I configuration, shown in Fig. 2a,c, the dislocation core structure consists of five-and seven-fold ring units, with two dumbbell pairs terminating the core. This type of coordination has been predicted for the edge dislocations in diamond cubic materials by Hornstra 30 and has been observed in classical simulations by Nandedkar and Narayan 31 . Type II dislocations contain eight-fold ring, as shown in  Fig. 2d-f) reveals the extra column between Cd 1 and Cd 2 , and the distance between the two Cd columns increases to around 0.72 nm. Based on the surrounding columns and DFT structural calculations, discussed later, we identify the extra column as Cd. The structure of Type III is given in Supplementary Fig. 3.
These Lomer dislocations can be seen as two 60° elemental dislocations that have merged from different {111} planes. The intersection angle is acute and the reaction can be written as: One interpretation of the difference between Type I and II cores is that the intersection may not happen at the vertex defined by the two {111} planes on which dislocations glide, but they may slip past each other before being 'locked-in' by strain fields. Future molecular dynamics calculations may help shed further light on this reaction.
An example of a type III dislocation, which consists of two cores separated by an intrinsic stacking fault, is shown in Fig. 3a. From the Burgers circuit surrounding the entire configuration, the Burgers vector is identified as b = (a/2)[110] (Fig. 3a). Images of the upper and lower dislocation cores, along with Burgers circuits, are shown separately in Fig. 3b,c. The projected component of the Burgers vector in Fig. 3b is (a/12) [11− 2]. This can readily be identified as a glide set Shockley 30° partial dislocation, with a corresponding Burgers vector of b 30° = (a/6)[2− 1− 1]. The core structure of the 30° partial dislocation consists of a single Cd atomic column. Earlier reports on the atomic and electronic structures of Shockley dislocations in CdTe have shown that the Cd dislocations act as shallow donors while the Te cores act as shallow acceptors, as indicated by the Fermi levels in the band structure diagrams 5 . Unlike the upper 30° partial dislocation, in the lower dislocation core in Fig. 3c, it is not readily apparent in which crystallographic direction the projected Burgers vector is pointing. In order to find the Burgers vector, we take the total Burgers vector b from Fig. 3a and subtract from it the 30° partial according to equation (2): 30 The Burgers vector of the lower core is therefore b′ = (a/6)[141]. According to previous studies in other materials, the unusual (a/6) < 411> partial dislocation can be seen as the reaction of a perfect 60° (a/2) < 110> and 90° (a/6) < 211> partial dislocation 32,33 .
As shown in equation (1), Lomer dislocations can be interpreted as two merged 60° elemental dislocations, and can be visualized using Thompsons tetrahedron notation ( Supplementary Fig. 5): In order to understand how Type III dislocation formed, we may consider two probable scenarios: i) the dissociation of a previously perfect Lomer dislocation ( Supplementary Fig. 5) or ii) the reaction of a 60° dislocation with an already dissociated 60° dislocation. Frank's b 2 criteria 34 , which compares elastic energy states of dislocations prior and post reaction, tells us that it is energetically more favorable for the latter process to take place. In fact, it can be seen that the undissociated Lomer dislocation has lower elastic energy compared to the final product after dissociation into two partials and a stacking fault. However, the elastic energy is reduced when two isolated dislocations, one of which is already dissociated, merge. In the following we will assume the latter scenario and discuss the reaction in more details.
The dissociation of a 60° dislocation on the (111) plane can be written as:°→°+°-  (4) or using the Thomson tetrahedron notation, as:°→ An intrinsic stacking fault is created in this step and the 30° partial is placed at the upper core of Fig. 3a. The second 60° dislocation reacts with the 90° partial (a/6) [11− 2] or Bδ ' in Supplementary Fig. 4 to form the partial dislocation at the lower core of Fig. 3c, described as:°+°→ The dislocation reactions of the whole Type III dislocation are, therefore, given as: The two partial dislocations can be seen to possess screw components of equal magnitude but opposite direction, such that the whole configuration is still of pure-edge character. The distance between the partial dislocations is on average ~4.0 nm. The dissociation can take place on either one of the two distinct {111} planes, as seen in our STEM imaging analysis. However, we have not observed Lomer dislocations dissociated on both planes simultaneously. We note that, to our knowledge, this type of dislocation has not been reported in CdTe before.
Geometrical phase analysis (GPA) to determine strain field. To quantify the local structural differences of the different types of Lomer dislocations, we perform a GPA of the experimental HAADF-STEM images measuring the strain fields associated with each dislocation. These strain fields will affect the electrical and electronic properties of semiconductors by distorting the local bonding character. Charge transport and recombination may be altered through this effect on the electrostatic potential and via associated Cottrell atmospheres 35 . GPA is a method for obtaining local phase, which can be directly related to displacement field distorting the lattice fringes in the HRTEM or HAADF-STEM images with respect to the reference lattice 36 . Figure 4a-c presents the strain tensor component ε xx of Type I, II, and III dislocations, respectively, and Fig. 4d present the εyy of Type III dislocation.
The strain fields of Lomer dislocations Type I and II extend over 3 nm radius showing the large region of influence around these cores. Type I and II dislocations show two compression-tension strain field pairs, in Fig. 4a,b. We note that GPA is not well suited to infer strain on the scale smaller than a unit cell and thus the fine features are likely an artifact. The long-ranged strain is real and we are particularly interested in its extent. Figure 4c,d shows the strain tensor component ε xx and ε yy based on the HAADF-STEM image of the Type III dislocation obtained from the same area as shown in Fig. 3a. There is no in-plane strain expected due to the stacking fault itself. GPA shows very localized strain region around the stacking fault, which is again due to the small scale and thus can be ignored. The 30° Shockley partial dislocation has a very compact core and does not produce noticeable strain in this projection, although its screw component would need to be considered separately.
The lower (a/6)[141] core strain field is clearly much wider than that of the upper 30° partial dislocation, but appears less extended than the strain field of Type I and II cores. Strain fields from adjacent dislocation cores along the bicrystal interface are very close to overlapping with adjacent fields, but no actual interaction between adjacent strain fields has been observed, due to the average spacing between dislocation cores of ~5.6 nm. Thus, the individual cores can be considered as discrete and isolated from each other but the interface as a whole presents a continuously strained region.
Electronic structure of CdTe dislocations. First principles DFT calculations provide information about the local electronic structures of the dislocation cores. The amount of charge transfer near the Type I and II dislocation cores, calculated using Bader charge analysis, is given in Fig. 5a,b, respectively. In the Type I dislocation (Fig. 5a), the Cd (Te) atoms in the bulk (i.e. far from the core), have an average charge of + 0.55 (− 0.55) |e|, due to the polar-covalent character of the tetrahedral bonds between Cd and Te. In comparison, Te atoms surrounding the dislocation core have ~0.2 |e| smaller absolute Bader charges and the Cd atoms that form a dumbbell with these Te atoms have ~0.1 |e| smaller absolute Bader charges. The most significant effect is seen in the Type II dislocation (Fig. 5b), where the extra Cd atom at the center of the core has ~0.4 |e| lower absolute charge, becoming almost neutral in sharp contrast to the Cd atoms far from the dislocation core. The deviation of the Bader charges at the dislocation core compared to the bulk values is attributable to the change in coordination, which can give rise to the formation of additional states in the band gap.
Accordingly, the electronic structure of the dislocations is analyzed in terms of density of states (DOS). We integrate the projected DOS of atoms in a cylindrical volume around the dislocation core, and compare the result with the DOS of the bulk reference area, which is selected from the same computational cell far from the dislocation cores. Figure 5c, Fermi level, compared to bulk reference in the simulation cell. Figure 5c,d clearly shows that the dislocation cores give rise to additional electronic states within the bulk band gap. While the bulk-like region in the calculation may not reproduce bulk CdTe exactly due to the residual strain field created by the dislocation cores, the bulk-like region exhibits a Kohn-Sham gap similar to that of bulk CdTe (~0.6 eV). The bulk-like region also shows negligible DOS within a gap of about 1.3 eV, similar to bulk CdTe.
The overall negative charge at the Type I dislocation core and positive charge at the Type II dislocation core shown in Fig. 5a,b are expected to lead to electron repulsion and attraction, respectively. Indeed, the radial profiles of the electrostatic potential (Fig. 5e,f) show a repulsive potential for electrons (i.e. attractive for holes) at the Type I core and attractive potential for electrons (repulsive for holes) at the Type II core. We thus see an opposing mechanism to recombination via the mid-gap states provided by the presence of charged Type I and II dislocation cores which can be expected to help separate photo-generated electron-hole pairs.
The amount of charge transfer near the Type III dislocation cores and stacking fault is shown in Fig. 6a, which indicate that the charge distribution around the stacking fault is the same as the bulk region away from the stacking fault. At the Cd-terminated dislocation core, as shown in Fig. 6b, the charge distribution is not significantly different compared to bulk region, with Cd atoms having only ~0.1 |e| smaller charge than the bulk. At the Te terminated core in Fig. 6c, the Te atoms surrounding the dislocation core have a charge of − 0.3 |e|, which is ~0.2 |e| smaller in absolute value compared to the bulk region. Interestingly, the Te column at the dislocation core is charge neutral with almost zero charge (+ 0.04 |e|). The pertinent changes in the local atomic environment at the dislocation core will affect the electronic structure of the CdTe. In Fig. 6d, we calculate the integrated projected DOS of atoms in cylindrical regions centered at the dislocation cores at each end of the stacking fault, and in a rectangular region in the middle of the stacking fault, and compared the result with the DOS of the bulk reference area, which is selected from the same computational cell far from the stacking fault. The actual regions, where the DOS integration is carried out, is given in Supplementary Fig. 6. Figure 6d shows that stacking fault and bulk-like regions have similar DOS around the Fermi level, indicating that stacking fault does not alter the electronic structure of CdTe significantly. The Te-terminated core, on the other hand, shows formation of additional states in the midgap region up to the conduction band minimum. The Cd-terminated core shows additional states at the band gap that are closer to the valence band maximum.
The different charge states at the Te and Cd terminated cores at each end of the stacking fault are expected to act differently on charge carriers. To estimate the response of these defects on charge carriers, we computed the planar-averaged electrostatic potential along a direction perpendicular to the stacking fault. We considered different locations that are labeled in Fig. 6a, but used the same direction for averaging the electrostatic potential, as shown in Supplementary Fig. 7. Figure 6e shows the electrostatic potential profile along line 1, containing the Te-terminated dislocation core, which exhibits a lower electrostatic potential around the Te-core. The electrostatic potential profile along line 2, which passes through the middle of the stacking fault, indicates that the stacking fault exhibits bulk-like electrostatic potential profile (Fig. 6f). Similarly, electrostatic potential along line 3, containing Cd-terminated core, also appears bulk-like. Consequently, stacking fault and the Cd-terminated core is expected to have insignificant effects on charge separation and photovoltaic efficiency, similar to intrinsic twin boundaries 37 .

Figure 5. Electronic structure of Type I and II dislocation cores.
Charge transfer around the dislocation cores calculated using Bader analysis is given in (a,b) for Type I and II dislocations, respectively. Radially integrated density of states (DOS) analysis is given in (c,d) for Type I and II, respectively. The DOS of a region far from the dislocation cores was included for a bulk reference. Radially averaged electrostatic potential profiles as a function of distance from the center of the dislocation cores for (e) Type I and (f) Type II. The proportion of Type I and Type II dislocations is more than 90%, since they are more stable than Type III dislocation according to Frank's b 2 criteria, as we discussed above. Type II dislocation needs a dangling bond, which means that it needs more energy to break the Cd-Te dumbbell structure, so the proportion of Type I is more than that of Type II.
Depending on whether the core has Cd or Te dangling bonds, Type I, II and III dislocation cores significantly changed the electronic structure of CdTe pertaining to photovoltaic properties, as shown in Figs 5 and 6. Here, we explore passivating the dangling bonds at the dislocation core by Cl atoms with the aim of understanding its role on improving photovoltaic efficiency of CdTe that was reported recently 38,39 . As Cl impurities prefer to segregate at dislocation cores and grain boundaries 38 , we substitute all of the Te atoms at the Type I core with Cl to observe the maximum effect of the changes in the electronic structure and fully relax the atomic positions. Figure 7a shows the relaxed structure of the Cl-doped Type I dislocation core and Bader charge distribution. Cl atoms had − 0.6 |e| charge and the nearest Cd atoms had a charge of + 0.7 |e|. The absolute charge values on Cd and Cl at the core are larger than the Type I core in Fig. 5a. The density of states of Cl-doped core compared to original Type I core is given in Fig. 7b, showing substantial reduction in the electronic states in the gap due to Cl doping, which indicates that Cl is effectively passivating the dangling bonds at the dislocation core. However, the amount of Cl considered here was not sufficient to completely eliminate the additional states at the midgap. On the other hand, the radially averaged electrostatic potential of the Cl-doped core does not show a significant change compared to the original Type I core as shown in Fig. 7c, such that Cl-doped core still had a negative charge and would expect to attract holes and facilitate charge separation. Therefore, Cl doping of dislocation cores reduces the midgap states but in the mean time maintains charge separation that would effectively be expected to increase the lifetimes and hence, the photovoltaic efficiency of CdTe.
In a previous study 5 , single Cd and Te columns have been found at the CdTe intragrain Shockley partial dislocation cores, where the stacking faults terminate. The electronic structure of the individual dislocation cores reveals that the Cd cores act as shallow donors while the Te cores act as shallow acceptors 5 . While each individual core modeled separately may appear to have benign characteristics, there is unmistakable charge transfer between those cores when they appear at the ends of the same stacking fault, which is consistent with our results. Moreover, our models, which are not idealized models of separated cores, but realistic models consisting of up to ~1900 atoms and shown to match the STEM images, showed that more generally, Type I, II and III dislocation cores are not electrically benign, as the Cd and Te terminating atoms at the dislocation cores are under-coordinated.
In summary, we create a [1-10]/(110) 4.8° tilt grain boundary via wafer bonding. This study reveals two types of Lomer (a/2)[110] undissociated dislocations with different atom configurations, and one kind of Lomer (a/2) [110] dislocation dissociated to an intrinsic stacking fault with a 30° Shockley partial and an unusual (a/6)[141] partial dislocations. We study the reaction mechanism combining the HAADF-STEM image, GPA strain maps, and Thomson tetrahedron. Structure models derived from the STEM images have been constructed and relaxed by DFT calculations. The relaxed structure is found to be in good agreement with the experiments images. DFT calculations show all three types dislocations present midgap states that are expected to increase recombination, but also electrostatic potential profiles that may aid charge separation. We find that Cl passivation at type I core reduces the midgap states without affecting charge separation. Further work will focus on electrical measurement to verify the solar cell performance on this kind of GB. The combination of STEM, GPA and DFT calculation to explore the atomic structure, dislocation formation, strain map and electronic structure of dislocation cores should be applicable to other complicated GBs in a wide range of II-VI and III-V semiconductors.

Methods
Bicrystal fabrication. High purity, one side polished, 5 × 10 mm 2 p-type CdTe (111) wafers with a thickness of 3 mm were used for the bonding experiments. Out-of-plane (2θ/ω scan) and in-plane (ϕ scan) X-ray diffraction (XRD) (Ultima III, Rigaku) was used to determine the crystal orientations of the planes and flats, respectively. The root mean square surface roughness (RMS) was measured by atomic force microscopy (AFM) (Dimension 3100, Veeco) in tapping mode. For the cleaning process, samples were rinsed in de-ionized water for 5 s, briefly dipped in hydrochloric acid (15% HCl) for 30 s and then rinsed in de-ionized water for 5 s. The samples were dried under a flow of nitrogen gas, and heated on a hot plate at 110 °C for 5 s to evaporate the de-ionized water absorbed on the surface. After that, the wafer pair was loaded into a vacuum bonder (EV 501, EV Group) immediately, to prevent oxidation of the CdTe surface. Bonding was performed at 400 °C for 20 hrs under a pressure of 1 MPa.
TEM imaging conditions and GPA. Local structure and diffraction patterns were acquired using HRTEM (2100F, JEOL) of Focused Ion Beam (FIB) (Nova 200, FEI) milled cross-sections. The atomic arrangement of the interface was characterized using HAADF-STEM (JEM-ARM200F, JEOL). Geometrical phase analysis (GPA) was performed based on the HAADF-STEM images using a GPA plug-in package, which was implemented in the Digital Micrograph software (Gatan).

Construction of dislocation core atomic structures using image analysis and empirical potentials.
Atomistic models of pairs (dipoles) of dislocation cores are constructed from HAADF-STEM images and empirical potentials. First, a black and white filter was applied on the raw HAADF-STEM image to remove the noise and enhance the bright spots located at atomic columns. The positions of bright spots that correspond to atomic columns were identified using image analysis techniques with the Scikit-Image package. To extract the coordinates of atomic columns corresponding to bright spots, we located the maximum peak intensities using a maximum filter with a determined threshold on the image. To remove multiple maximum peaks appearing between atomic columns, we calculated the nearest neighbors of each peak using a radius of 0.9 Å, and combined the peak locations that are closer than this value by taking the center of mass of the nearest neighbors of each peak. After locating the positions of atomic columns, we used the crystallographic information in the CdTe crystal to identify the separate Cd and Te columns and obtained the atomic structure. To build a periodic supercell, we juxtaposed the structure with its mirror image and created a dislocation dipole where necessary. Before first principles calculations, we relaxed the dislocation dipole structure using available Stillinger-Weber 40 and bond order empirical potentials 41 for CdTe. These calculations have lower computational costs compared to DFT calculations and enabled us to determine near-minimum-energy dislocation core structures which match well with the HAADF-STEM images. We considered many variations of atomic configuration at the dislocation cores by adding/subtracting atoms and changing the dimensions of the supercell so that the dislocation cores are not strongly affected by the periodic boundary conditions imposed on the system. In all cases, both interatomic potentials relaxed the atoms to similar positions. The atomic coordinates in the dislocation models are further relaxed using DFT.
The Type I dislocation dipole cell consisted of 1440 atoms with 720 Cd and 720 Te. Type II differs from Type I by only a single Cd/Te atom at the dislocation dipole cores. Type III consists of a stacking fault with two dislocation cores at the two ends. The Cd-terminated upper core structure was apparent from STEM image, but there is ambiguity in the lower core atomic structure. Initially, we created a supercell of 56 × 57 × 15 Å containing 1401 atoms (697 Cd and 704 Te), and attempted different variations of single Te, single Cd, Cd-Te dumbbell atomic columns at the lower core. We compared the relaxed structures with the STEM image and found that this size was insufficient for overcoming the periodic boundary effects, so that we increased the cell size to 70 × 64 × 15 Å containing 1917 atoms (957 Cd and 960 Te), which is sufficient for reducing strain imposed on atoms around the dislocation cores. We tested similar atomic and elemental variations at the lower dislocation core using this large simulation cell as well. Finally, we chose the structure with Te-termination at the lower core that had the best match with STEM around the dislocation core and stacking fault in order to determine the electronic properties.
DFT calculations. The electronic structure of dislocation cores was analyzed in terms of electronic density of states, charge transfer, and electrostatic potential change at and near the dislocation core using planewave-based density functional theory (DFT) method. All DFT calculations were carried out using Vienna Ab initio Simulation Package (VASP) 42,43 . Projector augmented wave (PAW) potentials 44 were used, and the exchange-correlation was treated with the generalized gradient approximation (GGA) parameterized by Perdew, Burke, and Ernzerhof (PBE) 45 . Calculations were carried out at the Γ -point in the Brilluoin zone using planewave kinetic energy cutoff of 343 eV. Atomic positions were relaxed until forces are below 0.05 eV/Å and electronic relaxation was converged to 10 −5 eV to ensure convergence to 1-2 meV/atom in total energy. The electronic charges on atoms in DFT calculations are obtained by partitioning the charge density according to atomic positions in the crystal based on the Bader charge analysis 46 , as implemented in VTST-scripts 47 . Some DFT calculations are performed on Extreme Science and Engineering Discovery Environment (XSEDE) resources 48 .