Zero point energy (ZPE) and entropic (TS) correction for heteroatom doped-graphenes (G, GB, GN, GP, and GSi) at 298 K.
In this chapter, we introduce the density functional theory (DFT)-based computational approaches to the study of various electrochemical reactions (hydrogen evolution reaction (HER), oxygen evolution reaction (OER), oxygen reduction reaction (ORR)) occurring on heterogeneous catalysis surfaces. A detailed computational approach to the theoretical interpretation of electrochemical reactions and structure-catalytic activity relationships for graphene-based catalysts will be discussed. The electrocatalytic activity of catalysis can be theoretically evaluated by overpotential value determined from free energy diagram (FED) of electrochemical reactions. By comparing electrocatalytic activity of systematically designed graphene-based catalysts, we will discuss the structure-catalytic activity relationships, especially the electronic and geometrical effects of heteroatom dopants.
With the climate change, fast consumption of fossil fuels, and environment situations due to carbon release, the research and development of clean energy is of vital importance in the coming decades. Promising applications of electrocatalysis for clean energy conversion, for example fuel cells, water electrolysis, metal-air batteries, and CO2 to fuel conversion, are the subjects of both extensive fundamental and utilitarian studies. These technologies play a crucial role in the future of sustainable energy utilization infrastructure, and thus huge research efforts have been dedicated to improving the electrocatalytic activity of these reactions, which include electrocatalytic oxygen reduction reaction (ORR), and hydrogen oxidation reaction (HOR) that occur on the cathode and anode of a hydrogen-oxygen fuel cell, respectively, and hydrogen evolution reaction (HER) and oxygen evolution reaction (OER) at the cathode and the anode of an electrolytic cell producing gaseous molecular hydrogen and oxygen, respectively. These reactions play an important role in regenerative fuel cells and dominate their overall performance [1, 2, 3, 4, 5]. Understanding the HER/OER/ORR mechanisms of various catalysts could provide design guidelines for material and process development, as well as facilitating the discovery of new catalysts. Above all, the detailed OER/ORR mechanisms in acid/alkaline environment are still being studied. Generally, OER/ORR can proceed in Langmuir-Hinshelwood (LH) or Eley-Rideal (ER) mechanisms . The LH mechanism comprises all reactive intermediates on the surface while the ER mechanism includes species from the electrolyte that reacts with the surface intermediate. Despite the controversy over the mechanism, ER mechanism is generally accepted with lower reaction energy barrier than that of LH mechanism , and many researchers have conducted theoretical studies on OER/ORR based on the ER mechanism. However, there are two feasible reaction pathways in ER mechanism, two-step pathway and four-step pathway depending on the relative stability of O* and OOH* intermediates generated after the adsorption of O2 on the catalyst . Thus, we sought to describe the detailed reaction pathway of the OER/ORR as well as proposing solutions for the determination of the preferred reaction pathway on the ER mechanism.
Precious metals such as platinum (Pt), iridium (Ir), and ruthenium (Ru)-based catalysts [9, 10, 11] are generally needed to promote the HER for the generation of hydrogen fuel from the electrochemical splitting of water, the ORR in fuel cells for energy conversion, and the OER in metal-air batteries for energy storage. Besides the requirement for high catalytic activity, other issues related to these catalysts are their limited reserves and comparatively high cost, which have precluded these renewable energy technologies from large-scale commercial applications. In this regard, huge amount of efforts has been devoted to develop novel electrocatalysts to completely or partially replace precious metal catalysts in energy technologies. Along with the intensive research efforts in developing nonprecious electrocatalysts to reduce or to replace precious metal catalysts, various carbon-based, metal-free catalysts have been extensively studied because they have unique advantages for designated catalysis due to their tunable molecular structures, abundance and strong tolerance to acid/alkaline environments when used as alternative HER/OER/ORR catalysts. A rapidly growing field of metal-free catalysis based on carbon-based materials has developed, and a substantial amount of literature in both on the theoretical and experimental fields has been generated. Recent studies have revealed that graphene,  graphite,  vertically aligned nitrogen-doped carbon nanotubes (VA-NCNTs),  heteroatom-doped CNTs,  and nitrogen-doped graphene sheets  have excellent catalytic performance. The presence of N in N-doped graphene leads to more chemically active sites, a high density of defects and high electrochemical activity. Due to these enhanced electronic properties, N-doped catalysts in the C network are attractive for a wide range of applications, including as metal-free catalysts for HER/OER/ORR in fuel cell systems. Recently, carbon nitride-based catalysts (C3N4 and C2N) with N-rich including both graphitic and pyridinic N moieties is a promising catalyst due to its competitiveness over a wide range of electrocatalyst processes, despite pure C3N4 and C2N itself being inert with regard to HER/OER/ORR activity. Here, we attempted to enhance the catalytic activity of graphene, C3N4 and C2N by introducing heteroatoms, which is an effective way to manipulate its electronic structure and electrochemical properties.
In this chapter, we will introduce metal-free bifunctional electrocatalysts of the heteroatom-doped graphenes (GXs, where G and X represent graphene and the heteroatom dopant) for HER  and the heteroatom-doped C3N4 (XY-C3N4s, where X and Y indicate the dopant and doping site on C3N4, respectively) for OER/ORR [18, 19, 20]. From the doping effect which shows better performance for HER/OER/ORR, we first present evidence that structural deformation and periodic lattice defects play the fundamental role in the HER activity of GXs by adjusting the electronic properties of graphene. We found that graphene doped with third row elements has higher HER activity with out-of-plane structural deformation compared to graphene doped with second row elements, in which graphene tends to maintain its planar structure. We systematically describe a structure-activity relationship in GXs for HER based on a thorough understanding of the effects of dopants, respectively. In addition, the third row elements-doped graphenes (GSi, GP and GS) show an interesting regularity described by a simple 3 N rule: GXs give outstanding HER activity with sustained metallic property when its primitive cell size has 3 × 3 N (N is integral) supercell size of pure graphene. Secondly, we describe not only the detailed OER/ORR mechanisms but also improved OER/ORR activity of C3N4 by introducing dopants such as P or S into the C3N4 matrix. Especially, we explore the causes of variation in HER/OER/ORR performance with respect to the type of dopant by comparing geometric and electronic structures of GXs and XY-C3N4s.
From these geometric and electronic structures, we demonstrated that GXs  and XY-C3N4 [18, 19, 20] show outstanding HER/OER/ORR activity with synergistic geometric and electronic effects, which coordinatively increase unsaturated sp3-C via structural deformation and improve electrical conductance by modulating the electronic structure with extra electrons from dopants. Our theoretical investigations suggest that the synergistic effect between geometric and electronic factors plays an important role in HER/OER/ORR catalytic activities. It can be emphasized that there is a close correlation between the geometric/electronic structure and HER/OER/ORR catalytic activities. This understanding of the structure-activity relationship will give an insight into the development of new highly efficient electrocatalytic materials.
2. Theoretical background
2.1. Hydrogen evolution reaction
HER is a multistep process that takes place on the surface of catalyst, and there are two proposed mechanisms: Volmer-Heyrovsky and Volmer-Tafel. Both Volmer-Heyrovsky and Volmer-Tafel mechanisms describe the hydrogen atom adsorption and hydrogen molecule desorption reactions among (1) an initial state , (2) an intermediate adsorbed state or , and (3) a final product state , where the and denote the active site and adsorbed hydrogen atom on the surface of the catalyst, respectively. Because the initial and final states are equivalent at equilibrium reduction potential, U = 0, the Gibbs free energy of the intermediate state, , has been considered as a major descriptor of the HER activity for a wide variety of catalysts. Therefore, the optimum value of should be zero for a spontaneous reaction without activation energy barrier. The Pt catalyst facilitates HER with a low activation energy, ∼0.09 eV . In this work, we have considered the Volmer-Tafel mechanism to calculate the Gibbs free energy of the intermediate state, , , , and with different hydrogen coverage at the active sites, , and , as shown in Figure 1. The Gibbs free energy of the adsorbed hydrogen is calculated as:
The gas phase values were from reference 17, while the values for the adsorbed species were taken from DFT calculations. The same values for the adsorbed species for all the N X N models were used, as vibrational frequencies have been found to depend much less on the surface than the bond strength.
where , and refer to the number of hydrogen atoms, adsorbed hydrogen on the surface, and hydrogen coverage, respectively.
In contrast to the single hydrogen reaction of the Volmer step, two hydrogen atoms mediate the Tafel step. Therefore, we obtain the Gibbs free energy of the intermediate state during the Volmer and Tafel steps with the following equations to determine the different hydrogen coverages of active sites.
2.2. Oxygen evolution reaction and oxygen reduction reaction
2.2.1. Reaction pathways in alkaline media
The generally acceptable OER mechanism is the four-electron associative mechanism in alkaline media. The four elementary steps of OER mechanism are described as follows:
where * represents the active site on the surface, (l) and (g) refer to liquid and gas phases, respectively, and O*, OH* and OOH* are adsorbed intermediates.
In contrast to OER, ORR can proceed either by a two-step or four-step pathways depending on the relative stability of O* and OOH* intermediates generated after the adsorption of O2 on the catalyst. The two-step pathway of ORR in alkaline environment is summarized using the following elementary steps,
whereas the four-step pathway has following elementary steps:
Looking at the elementary reaction steps of ORR, both reaction pathways lead to the same final products as 4OH− and the different intermediate states after the adsorption of O2 on the catalyst, O* in (Eq. (12)) and OOH* in (Eq. (14)). Therefore, it is worth mentioning that the transformation step of
2.2.2. Derivation of the free energy relations
The free energy change of each elementary reaction of OER in alkaline media can be expressed as follows:
Therefore, the free energy change of each elementary reaction can be calculated using (1) the chemical potentials of hydroxide, electron, liquid water and oxygen molecule () and (2) the free energies of each intermediate (, , and ) on the surface of catalyst (*). In the following paragraph, we describe the theoretical method to calculate the free energy change of each elementary reaction of the proposed reaction path based on the relation between the standard electrode potential and the chemical potential of hydroxides and electrons in alkaline environment.
2.2.3. The chemical potentials of , , , and
In alkaline environment, the standard oxygen reduction reaction is described as follows:
and the standard reduction potential (
The left side in (Eq. (24)), the chemical potentials of electron and hydroxide, could be derived further as follows:
at standard and equilibrium conditions . In the case of the right side in Eq. (24), the chemical potentials of liquid water and oxygen molecule (
Moreover, the chemical potential of oxygen molecule (
Because the experimental standard free energy change is −2.46 eV, the equation can be written as:
Therefore, the chemical potential of oxygen molecule (
Therefore, Eq. (24) can be written as:
Finally, we can obtain
|Species||E (eV)||ZPE (eV)||TS (eV)|
|H2O (0.035 bar)||−14.22||0.56||0.67|
Gas phase H2O at 0.035 bar was used as the reference state because at this pressure gas phase H2O is in equilibrium with liquid water at 300 K. The same values for the adsorbed species for all the models were used, as vibrational frequencies have been found to depend much less on the surface than the bond strength. We took the standard entropies from thermodynamic tables for gas phase molecules.
2.2.4. The free energies of each intermediate (, , and ) on the surface of catalyst (*)
The first step is the adsorption step of active site with a release of an electron:
where respectively and could be expressed by DFT energies:
Replacing (Eq. (33)), (Eq. (35)) and (Eq. (36)) in (Eq. (34)) we get:
The second step is oxidation of the OH* species to O* with release of water and an electron:
The relation for
The third step is represented by formation of the OOH* on top of oxygen with a release of an electron:
The relation for and in terms of DFT energies is similar to the relation for (Eq. (35)). The same equations are replaced gradually in the (Eq. (40)) as follows:
The last step is the evolution of oxygen molecule:
Therefore, we leave from the following equation:
Finally, the summation of , , , and should be 1.608 eV.
The reaction free energies of O*, OH* and OOH* species on the surface of catalyst are corrected by ZPE and TS. Moreover, the free energy relations of the ORR in alkaline media can also be explained in the same vein as the OER above mentioned .
2.2.5. Reaction pathways in acidic media
The generally acceptable OER mechanism is the four-electron associative mechanism in acidic media. The four elementary steps of OER mechanism are described as follows:
where * represents the active site on the surface, (l) and refer to liquid and gas phases, respectively, and O*, OH* and OOH* are adsorbed intermediates. The ORR can proceed completely via a direct four-electron process in which O2(g) is reduced directly to water H2O(g), without involvement of hydrogen peroxide. The ORR mechanism is summarized using the following elementary steps,
Here, we took the OER reactions ((44)–(47)) to derive the thermochemistry of both OER/ORR, because the ORR reactions (Eqs. ((48)–(51)) are inversed from the OER reactions (Eqs. ((44)–(47)). The catalytic activity of the OER/ORR processes can be determined by examining the reaction free energies of the different elementary steps.
2.2.6. Derivation of the free energy relations
The free energy change of each elementary reaction of OER can be expressed as follows:
Therefore, the free energy change of each elementary reaction can be calculated using (1) the chemical potentials of proton, electron, liquid water and oxygen molecule (
2.2.7. The chemical potentials of H+, and
In acidic environment, the standard hydrogen electrode is based on the redox half-cell,
and the standard reduction potential () is 0.00 V at T = 298.15 K. In terms of chemical potentials, the equation is expressed as follows:
In (Eq. (57)), the chemical potentials of proton, electron and hydrogen could be derived further as follows:
at standard and equilibrium conditions (
Another approximation is that for liquid water and oxygen molecule,
2.2.8. The free energies of each intermediate (, , and ) on the surface of catalyst (*)
The first step is the adsorption step of OH* on active site with a release of a proton and an electron:
where and are expressed by DFT energies:
Replacing (Eq. (29)), (Eq. (62)), (Eq. (64)) and (Eq. (65)) in (Eq. (63)) we get:
The second step is oxidation of the OH* species to O* with release of a proton and an electron:
The relation for
The third step is represented by formation of the O* on top of oxygen with a release of a proton and an electron:
The relation for
The last step is the evolution of oxygen molecule:
Therefore, we leave from the following equation:
Finally, the summation of , , , and should be 4.92 eV.
The reaction free energies of OH*, O* and OOH* species on the surface of catalyst considered by ZPE and TS. Moreover, the free energy relations of the ORR in acidic media can also be explained in the same vein as the OER above mentioned.
2.3. Free energy diagram (FED) and Overpotential()
We can deduce an important parameter of electrocatalytic activity from the calculated , the magnitude of the potential-determining step (GHER and GOER/ORR) in consecutive reaction steps. This is the specific reaction step with the largest in the HER and OER/ORR elementary reaction steps, that is, the concluding step to achieve a downhill reaction in the free energy diagram (FED) with increasing potential N, as shown in Figure 2 for HER and Figure 3 for OER/ORR:
The theoretical overpotential at standard conditions is then given by (Eq. (76)) in acidic and alkaline conditions, respectively:
In the case of HER pathways, the Volmer-Heyrovsky reaction is an electrochemical reaction involving electrons at all steps, as shown in Figure 2(a). Therefore, when the external potential () is applied, the chemical potential of steps involving electrons shifts by the external potential (), and the Volmer-Heyrovsky reaction changes from uphill to downhill reaction. In contrast to the Volmer-Heyrovsky reaction, the Volmer-Tafel reaction is an electrochemical reaction involving electrons only in the Volmer step, as shown in Figure 2(b). Therefore, when the adsorption of hydrogen is unstable, the downhill reaction is possible by applying the external potential (), but when the adsorption of hydrogen is stable, the reduction of H+ is spontaneous reaction, and the desorption of H2 in the Tafel step is the rate determining step. However, the desorption of H2 is a thermodynamic reaction, which is not involved in electrons, not an electrochemical reaction. In the case of OER/ORR pathways, the theoretical overpotential (ηΟΕΡ/ΟΡΡ) represents the relative stability of the intermediates between (4OH− + *)alkaline, (2H2O(l) + *)acidic and (O2(g) + *), and vice versa, as shown in Figure 2(c)-(e). In addition, it can be calculated by applying standard density functional theory (DFT) in combination with the computational standard hydrogen electrode (SHE) model. Because the equilibrium reduction potential is and , the chemical potential difference between (4OH−)alkaline, (2H2O(l))acidic and O2(g) should be 1.608 and 4.920 V, respectively. However, the actual catalytic behaviors deviate from the ideal case due to correlation with binding energies of the intermediates. Therefore, most catalysts require overpotential (ηHER and ηOER/ORR) in order to achieve an overall downhill reaction. Consequently, the ηHER and ηOER/ORR the important indicators of catalytic activities of a catalyst, and the lower ηHER and ηOER/ORR indicate a thermodynamically superior catalyst.
3. Structure-catalytic activity relationships
3.1. Heteroatom doped graphene (GX) for HER catalyst
Figure 3 shows a schematic of the heteroatom doped-graphene (GX, where G means graphene and X is B, N, Si, P, and S dopant) structures. We have used the second row elements (B and N) and the third row elements (Si, P, and S) in the periodic table in order to investigate the structural and electronic doping effects on HER activity because the third row elements are relatively larger than the second row elements and p- and n-type doping effects can be expected from the electron deficient B, and electron rich N, P, and S elements give in-plane and out-of-plane structures, respectively, due to the size of the dopants.
3.1.1. Structural and electronic doping effects on HER activity
On the atomic orbital hybridization characters of adjacent carbon atoms of dopant, the natural bond orbital (NBO) analysis show an increased p orbital contribution from sp2 to sp3 hybridization of carbons adjacent to the dopant due to structural deformation from in-plane to out-of-plane. Compared to sp2-hybridized carbon, sp3-hybridized carbons more readily form an extra a hydrogen atom without additional structural change. Therefore, in out-of-plane structures, the subsequent two hydrogen atoms prefer to bind to only sp3 hybridized carbons adjacent to the dopant. However, in the case of in-plane structures having only sp2-hybridized carbons, the first hydrogen atom should result in structural deformation to form sp3-hybridized carbons. Therefore, the first hydrogen adsorption on the in-plane structure is less favorable than the reaction on out-of-plane structures. The second hydrogen atom can favorably bind to sp3-hybridized second neighboring carbons of dopant. Consequently, structural deformation with dopants that are third row elements is associated with improved HER activity due to atomic orbital hybridization, as shown in Figure 3(b). Looking at the electronic structures of GXs, p- and n-type doping effects are also expected from electron deficient and rich elements. In the case of in-plane GXs, electron deficient boron shifts its band structure up by withdrawing an electron from graphene, and electron rich nitrogen shifts its band structure down by donating an electron to graphene. Interestingly, in the case of out-of-plane GXs, electron rich phosphorous and sulfur dopant have no associated band shift. The origin of the flat band can be understood based on the localization of an extra electron onto the dopant site. In order to verify the relationship between geometric and electronic structures of GXs, we have systematically changed the structures of GXs from in-plane to out-of-plane structures and vice versa. These calculations clearly show that dopant can produce n- and p-type doping states as well as a localized state depending on the structure of GXs induced by the type of dopants. When the out-of-plane (in-plane) deformation is applied in the in-plane (out-of-plane) GXs, band structures change from the p-type doping state (localized state) to the localized state (p-type doping state). Therefore, it is worth mentioning that the localized electronic states can be associated with physical regularity of HER activities on the out-of-plane GXs.
3.2. Heteroatom doped C3N4 (XY-C3N4) for OER/ORR bi-functional catalyst
In this work, we have systematically investigated a metal-free bifunctional electrocatalyst of heteroatom-doped carbon nitride (XY-C3N4, where X and Y indicate the dopant and doping site on C3N4, respectively) for oxygen evolution and oxygen reduction reactions (OER and ORR) in alkaline media, considering the possible reaction pathways based on the Eley-Rideal (ER) mechanism as well as the doping effects on electrocatalytic activity. Moreover, we determined that the relative stability of O* and OOH* intermediates was a key factor determining the ORR pathway; accordingly, ORR follows a two-step reaction pathway governed by O* rather than a four-step reaction pathway governed by OOH*. Firstly, we calculated relative Ef and the most stable structures of heteroatom-doped C3N4 structures, XY-C3N4, where X = P, S, or PS and Y = C or N [18, 20]. By considering the Ef (formation energy) results as well as the various doping sites, we determined the most stable XY-C3N4 structures, PCA-C3N4, PNB-C3N4, SCB-C3N4, SNB-C3N4, PCASCB-C3N4, PNASNB-C3N4, and PCASNB-C3N4, for the investigation of OER/ORR bifunctional electrocatalytic activities and their reaction mechanisms in alkaline media.
3.2.1. OER and ORR catalytic activity on XY-C3N4
Figure 4 (a) and (c) show the volcano plot of OER and ORR in alkaline media at all possible active sites on XY-C3N4, which represents the apparent catalytic activity, respectively. This theoretical analysis reveals that the PCSC-C3N4 structure has minimum ηOER (0.42 V) and ηORR (0.27 V). The ηOER/ORR value of PCSC-C3N4 is comparable to those of the best conventional catalysts (∼0.42 V for OER on RuO2 and ∼0.45 V for ORR on Pt). Most elementary steps in OER have an uphill reaction at 0.00 V and at equilibrium potential of 0.402 V. Analyzing the FED at the equilibrium potential of 0.402 V, we can define the rate determination step as the formation of OOH* from O* on PCSC-C3N4, as shown in Figure 4 (b). Therefore, the at the rate determination step, 0.42 V, become the ηOER for facilitating the OER as a spontaneously downhill reaction. By applying the increased electrode potential of 0.82 V, all elementary reactions can be downhill reactions for spontaneous OER. In the case of the ORR, we considered the O2 adsorption free energy (ΔGO2*) and ηORR are highly correlated to determine ORR activity. If the O2 adsorption reaction is very difficult to achieve as an initiation step in ORR, that material cannot have ORR activity, even though it can have a possibility to have a very small ηORR. Therefore, we defined a new indicator of ORR activity as figure of merit (FOM), which is defined by FOM = −(ηORR + ΔGO2*). Using FOM, we constructed a volcano plot of ORR at all possible active sites on XY-C3N4 in order to compare the ORR catalytic activity, as shown in Figure 4 (c), where the PCSC-C3N4 structure have minimum ηORR of 0.27 V.
3.2.2. Structural and electronic doping effects on OER/ORR activity
The synergistic effect of P, S co-doping can be explained based on the geometric and electronic effects of heteroatoms [17, 18, 19, 20]. Considering the atomic size of heteroatoms, the relatively larger S and P dopants can cause structural deformation of XY-C3N4, which improves the OER and ORR activity by enhancing the stability of intermediates with increasing p orbital character of active sites adjacent to the dopant from sp2 to sp3 hybridization. Compared to sp2 hybridized orbitals on pure C3N4, sp3 character of active sites on XY-C3N4 is more suitable for forming chemical bonds with intermediate species. Therefore, in out-of-plane structures, the intermediates prefer to bind to sp3-hybridized active sites. Moreover, to verify the relationship between geometric and electronic structures of XY-C3N4, we intentionally changed the structures of XY-C3N4 from in-plane to out-of-plane to increase the activity for binding intermediate species by increasing the sp3 character of active sites on XY-C3N4. We investigated density of state (DOS) of pure C3N4 and XY-C3N4 (X = P, S, or PS and Y = N) to determine at what condition the OER/ORR exhibit outstanding performance [18, 20]. It can be expected that electron-rich heteroatom doping will induce electronic structure changes from a non-metallic doping effect to a metallic doping effect due to a downward band shift. As a result, PCSC-C3N4 shows the best OER/ORR activity by maintaining a metallic property despite the presence of out-of-plane deformation. Consequently, it can be emphasized that there is a close correlation between the electronic/geometric structure and OER/ORR catalytic activities, and the best bifunctional OER/ORR catalytic activity of P,S co-doped PCSC-C3N4 is attributed to a synergistic effect between the electronic and geometric effects.
We have systematically investigated the detailed mechanisms of HER/OER/ORR and the synergistic effect between geometric and electronic factors plays an important role in HER/OER/ORR catalytic activities of GXs and XY-C3N4. In this work, we demonstrated that the HER/OER/ORR activity of GXs and XY-C3N4 can be modulated by structural and electronic factors, including structural deformation with dopants. These structural and electronic factors enhance adsorbent binding strength during the reactions in HER/OER/ORR by generating sp3-hybridized atoms and facilitating charge transfer between adsorbents and GXs/XY-C3N4 as metallic properties. Additionally, we re-evaluated the generally accepted ER mechanism of OER/ORR by comparing the stability of intermediates governing the reactions, where we found that the OER/ORR respectively follows four-step and two-step reaction pathway. We also elucidated the importance of the adsorption free energy (ΔGO2*) in ORR activity. Considering the ΔGO2* with a FOM, FOM = −(ηORR + ΔGO2*), we successfully represented the ORR activity of XY-C3N4. We believe that the understanding of the detailed mechanism as well as the relationship of structure-electrocatalytic activity of the HER/OER/ORR can facilitate development of new electrocatalytic materials.
This research was supported by grants from the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT, and Future Planning (NRF-2018R1A2B6006320). This work was also supported by the Supercomputing Center/Korea Institute of Science and Technology Information with supercomputing resources including technical support (KSC-2017-C3-0032).