Local self-interaction correction method with a simple scaling factor

In the PZSIC methodPerdew and Zunger (1981), SIE is removed on an orbital by orbital basis from the DFA energy as

$\displaystyle E^{PZSIC-DFA}$

$\displaystyle=E^{DFA}[\rho_{\uparrow},\rho_{\downarrow}]-\sum_{i\sigma}^{occ}\left\{U[\rho_{i\sigma}]+E_{XC}^{DFA}[\rho_{i\sigma},0]\right\},$

(1)

where $i$ is the orbital index, $\sigma$ is the spin index, $\rho$ ($\rho_{i\sigma}$) is the electron density (local orbital density), $U[\rho_{i\sigma}]$ is the exact self-Coulomb energy, and $E_{XC}^{DFA}[\rho_{i\sigma},0]$ is the self-exchange-correlation energy for a given DFA XC functional. Perdew and Zunger applied this scheme to atoms using the Kohn-Sham orbitals. For larger systems the Kohn-Sham orbitals can be delocalized which would result in the violation of size extensivity. Therefore local orbitals are required. This was recognized long ago by Slater and WoodSlater and Wood (1970) in 1971 and was also emphasized by GopinathanGopinathan (1977) in the context of self-interaction-correction of Hartree-Slater method and later by Perdew and Zunger in the context of approximate Kohn-Sham calculations. Subsequent PZSIC calculations by Wisconsin groupHeaton et al. (1983); Harrison et al. (1983); Pederson et al. (1984, 1985); Harrison et al. (1983) used local orbitals in variational implementation. It was shown by Pederson and coworkers that local orbitals used in the Eq. (1) must satisfy the localization equations (LE) for variational minimization of energy. The LE for the orbitals $\phi_{i\sigma}$ is a pairwise condition and is given as

$$\langle\phi_{i\sigma}|V_{i\sigma}^{SIC}-V_{j\sigma}^{SIC}|\phi_{j\sigma}\rangle=0.$$

(2)

In the FLOSIC approach, FLOs are used in stead of directly solving the Eq. (2). First, FOs $\phi^{FO}$ are constructed with the density matrix and spin density at special positions in space called Fermi orbital descriptor (FOD) positions as

$$\phi_{i}^{FO}(\vec{r})=\frac{\sum_{j}^{N_{occ}}\psi_{j}(\vec{a_{i}})\psi_{j}(\vec{r})}{\sqrt{\rho_{i}(\vec{a_{i}})}}.$$

(3)

Here, $i$ and $j$ are the orbital indexes, and $\psi$ is the KS orbital, $\rho_{i}$ is the electron spin density, and $\vec{a_{i}}$ is the FOD position. The FOs are then orthogonalized with the Löwdin’s scheme to form the FLOs. The FLOs are used for the calculation of the SIC energy and potential. In this method, the optimal set of FLOs are found by finding the FODs that minimizes total energy. This optimization process is similar to that for geometry optimization. We note that FLOs can be used in all three SIC (PZSIC, OSIC, and LSIC) methods.

As mentioned in Sec. I, PZSIC tends to overcorrect the DFA energies and several modifications to PZISC were proposed to scale down the PZSIC correction. In the OSIC method of Vydrov et alVydrov et al. (2006) mentioned in Introduction Eq. (1) is modified to

$\displaystyle E^{OSIC-DFA}$

$\displaystyle=E_{XC}^{DFA}[\rho_{\uparrow},\rho_{\downarrow}]-\sum_{i\sigma}^{occ}X_{i\sigma}^{k}\left(U[\rho_{i\sigma}]+E_{XC}^{DFA}[\rho_{i\sigma},0]\right),$

(4)

where each local orbitalwise scaling factor $X_{i\sigma}^{k}$ is defined as

$$X^{k}_{i\sigma}=\int z_{\sigma}^{k}(\vec{r})\rho_{i\sigma}(\vec{r})d^{3}\vec{r}.$$

(5)

Here, $i$ indicates the orbital index, $\sigma$ is the spin index, $z_{\sigma}$ is the iso-orbital indicator, and $k$ is an integer. The quantity $z_{\sigma}$ is used to interpolate the single-electron regions ($z_{\sigma}=1$) and uniform density region ($z_{\sigma}=0$). In their original work, Vydrov et al. used $z_{\sigma}=\tau_{\sigma}^{W}/\tau_{\sigma}$ to study the performance of OSIC with various XC functionals where $\tau_{\sigma}^{W}(\vec{r})=|\vec{\nabla}\rho_{\sigma}(\vec{r})|^{2}/(8\rho_{\sigma}(\vec{r}))$ is the von Weiszäcker kinetic energy density and $\tau_{\sigma}(\vec{r})=\frac{1}{2}\sum_{i}|\vec{\nabla}\psi_{i\sigma}(\vec{r})|^{2}$ is the non-interacting kinetic energy density. Satisfying the gradient expansion in $\rho$ requires $k\geq 1$ for LSDA, $k\geq 2$ for GGAs, and $k\geq 3$ for meta-GGA. Vydrov et al., however, used various values of $k$ to study its effect on the OSIC performance.

In their subsequent work, Vydrov et al.Vydrov and Scuseria (2006) used

$$w_{i\sigma}^{k}(\vec{r})=\left(\frac{\rho_{i\sigma}(\vec{r})}{\rho_{\sigma}(\vec{r})}\right)^{k},$$

(6)

the weight used by Slater in averaging Hartree-Fock potential, as a scaling factor instead of kinetic energy ratio. They repeated the OSIC calculations using $w_{i\sigma}$ in place of $z_{\sigma}$ in Eq. (5). Notice that Eq. (6) contains a local orbital index, this weight is thus an orbital dependent quantity. $w_{i\sigma}$ approaches unity at single orbital regions since $\rho_{\sigma}(\vec{r})=\rho_{i\sigma}(\vec{r})$ at this limit. Similarly, $w_{i\sigma}$ approaches zero at many-electron region since $\rho_{\sigma}(\vec{r})\gg\rho_{i\sigma}(\vec{r})$ at this condition. It was reported that the OSIC with Eq. (6) showed comparable performance as $z_{\sigma}=\tau_{\sigma}^{W}/\tau_{\sigma}$ despite of its simpler form.

Though OSIC had some success in improving the performance with SIC, the approach leads to parameter $k$ dependent performance. Also, it gives $-X_{HO}/r$ asymptotic potential instead of $-1/r$ for finite neutral systems and it results in inaccurate description of dissociation behaviorRuzsinszky et al. (2007). In addition, many-electron SIE and fractional-charge dissociation behavior of positively charged dimers reemerge with the OSICRuzsinszky et al. (2006). The recent LSIC method by Zope et al. applies the SIC in a different way than OSIC and retains desirable beneficial features of PZSIC. In LSIC, the SIC energy density is scaled down locally as follows,

$\displaystyle E_{XC}^{LSIC-DFA}$

$\displaystyle=E_{XC}^{DFA}[\rho_{\uparrow},\rho_{\downarrow}]-\sum_{i\sigma}^{occ}\left(U^{LSIC}[\rho_{i\sigma}]+E_{XC}^{LSIC}[\rho_{i\sigma},0]\right),$

(7)

where

$$U^{LSIC}[\rho_{i\sigma}]=\frac{1}{2}\int d^{3}\vec{r}\,z_{\sigma}(\vec{r})^{k}\,\rho_{i\sigma}(\vec{r})\int d^{3}\vec{r^{\prime}}\,\frac{\rho_{i\sigma}(\vec{r^{\prime}})}{|\vec{r}-\vec{r^{\prime}}|},$$

(8)

$$E_{XC}^{LSIC}[\rho_{i\sigma},0]=\int d^{3}\vec{r}\,z_{\sigma}(\vec{r})^{k}\,\rho_{i\sigma}(\vec{r})\epsilon_{XC}^{DFA}([\rho_{i\sigma},0],\vec{r}).$$

(9)

LSIC uses an iso-orbital indicator to apply SIC pointwise in space. An ideal choice of iso-orbital indicator should be such that LSIC reduces to DFA in uniform gas limit and reduces to PZSIC in the pure one-electron limit. To demonstrate the LSIC concept Zope et al. used $z_{\sigma}=\tau_{\sigma}^{W}/\tau_{\sigma}$ as an iso-orbital indicator. In this study, however, we use $w_{i\sigma}(\vec{r})=\rho_{i\sigma}(\vec{r})/\rho_{\sigma}(\vec{r})$ in place for $z_{\sigma}$ in Eqs. (8) and (9). We refer to the LSIC with $z_{\sigma}(\vec{r})$ as LSIC($z$) and LSIC with $w_{i\sigma}(\vec{r})$ as LSIC($w$) to differentiate the two cases.

All of the calculations were performed using the developmental version of FLOSIC codeZope et al. ; Yamamoto et al. , a software based on the UTEP-NRLMOL code. PZSIC, OSIC, and LSIC methods using FLOs are implemented in this code. FLOSIC/NRLMOL code uses Gaussian type orbitalsPorezag and Pederson (1999) whose default basis sets are in similar quality as quadruple zeta basis sets. We used the NRLMOL default basis sets throughout our calculations. For calculations of atomic anions, long range s, p, and d single Gaussian orbitals are added to give a better description of the extended nature of anions. The exponents $\beta$ of these added single Gaussians were obtained using the relation, $\beta(N+1)=\beta(N)^{2}/\beta(N-1)$, where $N$ is the $N$-th exponent. FLOSIC code uses a variational integration meshPederson and Jackson (1990) that provides accurate numerical integration.

In this work, our focus is on the LSDA functional because LSIC applied to LSDA is free from the gauge problemBhattarai et al. (2020) unlike GGAs and meta-GGAs where a gauge transformation is needed since their XC potentials are not in the Hartree gauge. We used an SCF energy convergence criteria of $10^{-6}$ Ha for the total energy and an FOD force tolerance of $10^{-3}$ Ha/bohr for FOD optimizations in FLOSIC calculations. For OSIC and LSIC calculations, we used respective FLOSIC densities and FODs as a starting point and performed a non-self-consistent calculation of energy on the FLOSIC densities. Several values for the scaling power $k$ are used in the LSIC($w$) and OSIC($w$) calculations. The additional computational cost of the scaling factor in OSIC and LSIC is very small compared to a regular FLO-PZSIC calculation.

In this section, we present our results on total energies, ionization potentials, and electron affinities for atoms.

To study the performance of LSIC($w$) for molecules, first, we calculated the atomization energies (AEs) of 37 selected molecules. Many of these molecules are subset of the G2/97 test setCurtiss et al. (1991). The 37 molecules set includes systems from the AE6 setLynch and Truhlar (2003), small but a good representative of the main group atomization energy (MGAE109) setPeverati and Truhlar (2011). The AEs were calculated by taking the energy difference of fragment atoms and the complex, that is, $AE=\sum_{i}^{N_{atom}}E_{i}-E_{mol}>0.$ $E_{i}$ is the total energy of an atom, $E_{mol}$ is the total energy of the molecule, and $N_{atom}$ is the number of atoms in the molecule. The calculated AEs were compared to the non-spin-orbit coupling reference valuesPeverati and Truhlar (2011) for AE6 set and to the experimental valuesNational Institute of Standards and Technology for the entire set of 37 molecules. The percentage errors obtained through various methods are shown in Fig. 4. The overestimation of AEs with PZSIC-LSDA due to overcorrection is rectified in LSIC($w$). We have summarized MAEs and mean absolute percentage errors (MAPEs) of AE6 and 37 molecules from G2 set in Table 5. For AE6 set, MAEs for PZSIC, LSIC($z$), LSIC($w,k=1$), and OSIC($w,k=1$) are $57.9$, $9.9$, $13.8$, and $33.7$ kcal/mol respectively. The MAE in LSIC($z$) is about 4 kcal/mol larger than LSIC($w,k=1$) but substantially better than the PZSICs or OSIC($w$). For the larger $k$ in LSIC($w$), however, the performance starts to degrade with consistent increase in the MAE of $33.5$ kcal/mol for $k=4$. This is in contrast to OSIC where the performance improves for $k=2$ and $3$ compared to $k=1$. The scaling thus affect differently in the two methods. OSIC($w,k=1$) tends to slightly underestimate total energies. By increasing $k$, total energies shift toward the LSDA total energies and improves performance for moderate increase in $k$. On the contrary, total energies are slightly overestimated for LSIC($w,k=1$), and increasing $k$ makes the energies deviate away from the accurate estimates. OSIC($w,k=3$) and LSIC($w,k=1$) have a similar core orbital SIC energy. In their study of OSIC($w$), Vydrov and ScuseriaVydrov and Scuseria (2006) used values of $k$ up to 5 and found the smallest error of $k=5$ (MAE, $11.5$ kcal/mol). But we expect the OSIC performance to degrade eventually for large $k$ since increase in $k$ results in increase in quenching of the SIC correction thus the results will eventually approach to those of DFA, in this case LSDA. For the full set of 37 molecules, PZSIC, LSIC($z$), LSIC($w,k=1$), and OSIC($w,k=1$) show the MAPEs of $13.4$, $6.9$, $9.5$ and $11.9$%, respectively. OSIC($w$) shows a slight improvement in MAPE for $k=2$ and $3$. For the larger set, LSIC($w$) consistently shows smaller MAPEs than OSIC($w$) for $k=1-3$. All four values of $k$ with the LSIC($w$) in this study showed better performance than PZSIC for the 37 molecules set.

Accurate description of chemical reaction barrier is challenging for DFAs since it involves calculation of energies in non-equilibrium situations. In most of the cases, the saddle point energies are underestimated since DFAs do not perform well for a non-equilibrium state that involves a stretched bond. This shortcoming of DFAs in a stretched bond case arises from SIE; when an electron is shared and stretched out, SIE incorrectly lowers the energy of transition state. SIC handles the stretched bond states accurately and provides a correct picture in chemical reaction paths. We studied the reaction barriers using the BH6Lynch and Truhlar (2003) set of molecules for LSIC($w$) method. BH6 is a representative subset of the larger BH24Zheng et al. (2007) set consisting of three reactions OH + CH₄ $\rightarrow$ CH₃ + H₂O, H + OH $\rightarrow$ H₂ + O, and H + H₂S $\rightarrow$ H₂ + HS. We calculated the total energies of left- and right-hand side and at the saddle point of these chemical reactions. The barrier heights for the forward (f) and reverse (r) reactions were obtained by taking the energy differences of their corresponding reaction states.

The mean errors (MEs) and MAEs of computed barrier heights against the reference valuesLynch and Truhlar (2003) are compared in Table 6. MAEs for PZSIC, LSIC($z$), LSIC($w,k=1$), and OSIC($w,k=1$) are $4.8$, $1.3$, $3.6$, and $3.6$ kcal/mol, respectively. PZSIC significantly improves MAE compared to LSDA (MAE, 17.6 kcal/mol), LSIC($w,k=1$) further reduces the error from PZSIC. Its ME and MAE indicate that there is no systematic underestimation or overestimation. LSIC($w,k=1$) also further improves the PZSIC numbers but not to the same level as LSIC($z$). For $k\geq 2$, MAEs increases systematically for LSIC($w,k\geq 2$) though small MEs are seen for LSIC($w,k=2,3$). The performance deteriorates for $k>2$ beyond that of PZSIC. OSIC($w$) shows marginally better performance than PZSIC. Vydrov and ScuseriaVydrov and Scuseria (2006) showed that the best performance is achieved with $k=1$ (MAE, $3.5$ kcal/mol). The performance improvement with OSIC is not as dramatic as LSICs in terms of MEs and MAEs where the rather large MEs are seen i Overall LSIC($w$) performs better than OSIC($w$) for barrier heights.

A pronounced manifestation of SIE is seen in dissociation of positively charged dimers $X_{2}^{+}$. SIE causes the system to dissociate into two fractionally charges cations instead of $X$ and $X^{+}$. Here we use the SIE4x4Goerigk et al. (2017) and SIE11Goerigk and Grimme (2010) sets to study the performance of LSIC($w$) and OSIC($w$) in correcting the SIEs. The SIE4x4 set consists of dissociation energy calculations of four positively charged dimers at varying bond distances $R$ from their equilibrium distance $R_{e}$ such that $R/R_{e}$ = 1.0, 1.25, 1.5 and 1.75. The dissociation energy $E_{D}$ is calculated as

$$E_{D}=E(X)+E(X^{+})-E(X_{2}^{+}).$$

(11)

The SIE11 set consists of eleven reaction energy calculations: five cationic reactions and six neutral reactions. These two sets are commonly used for studying the SIE related problems. The calculated dissociation and reaction energies are compared against the CCSD(T) reference valuesGoerigk et al. (2017); Goerigk and Grimme (2010), and MAEs are summarized in Table 7. For the SIE4x4 set, PZSIC, LSIC($z$), LSIC($w,k=1$), and OSIC($w,k=1$) show MAEs of $3.0$, $2.6$, $4.7$ and $5.2$ kcal/mol. LSIC($z$) provides small improvement in equilibrium energies while keeping accurate behavior of PZSIC at the dissociation limit resulting in marginally better performance. LSIC($w$) shows errors a few kcal/mol larger than PZSIC. This increase in error arises because LSIC($w$) alters the (NH₃)${}_{2}^{+}$ and (H₂O)${}^{+}_{2}$ dissociation curves. In LSIC($z$) the scaling of SIC occurs mostly for the core orbitals (Cf. Table 2) whereas LSIC($w$) also includes some noticeable scaling down effect from valence orbitals. This different scaling behavior seems to contribute to different dissociation curves. Lastly, OSIC($w$) has a slightly larger error than LSIC($w$).

For the SIE11 set, MAEs are $11.5$, $4.5$, $8.3$, and $11.1$ kcal/mol for PZSIC, LSIC($z$), LSIC($w,k=1$), and OSIC($w,k=1$), respectively. All scaled-down approaches we considered, LSIC($z$) and LSIC($w$), and OSIC($w$) showed performance improvement over PZSIC. LSIC($z$) shows the largest error reduction by 60%, while LSIC($w,k=1$) shows 28% decrease in error with respect to PZSIC. OSIC($w$) with $k=1-3$ has slightly smaller MAEs within 1 kcal/mol of PZSIC. LSIC($z$) method improves cationic reactions more than neutral reactions with respect to PZSIC. Increase in $k$ beyond 2 results in too much suppression of SIC and leads to increase in error for LSIC($w,k\geq 2$). LSIC($w$) yielded consistently smaller MAEs than OSIC($w$) but larger than LSIC($z$) over the whole SIE11 reactions.

Finally, we show the ground-state dissociation curves for H${}^{+}_{2}$ and He${}^{+}_{2}$ in Fig. 5. As previously discussed in literature Ruzsinszky et al. (2005), DFAs at large separation cause the complexes to dissociate into two fragments atoms. PZSIC restores the correct dissociation behavior at the large separation distance. When LSIC is applied, the behavior of PZSIC at the dissociation limit is preserved in both LSIC($z$) and present LSIC($w$). For H${}_{2}^{+}$, a one-electron system, LSIC reproduces the identical behavior as PZSIC [Fig. 5 (a)]. For He${}_{2}^{+}$, a three-electron system, LSIC applies the correction to PZSIC only near equilibrium regime [Fig. 5 (b)]. LSIC brings the equilibrium energy closer to the CCSD energy compared to PZSIC energy. The implication of Fig. 5 is that the present scaling factor $w$ performs well in differentiating the single-orbital like regions and many-electron like regions.

Kamal et al.Sharkas et al. (2020) recently studied binding energies of small water clusters using the PZSIC method in conjunction with FLOs to examine the effect of SIC on binding energies of these systems. Water clusters are bonded by weaker hydrogen bonds and provide a different class of systems to test the performance of the LSIC method. Earlier studies using LSIC($z$) on the polarizabilities and ionization have shown that LSIC($z$) provides an excellent descriptions of these properties when compared to the CCSD(T) resultsVargas et al. (2020); Akter et al. . Here, we study the binding energies of the water clusters. We find that the choice of iso-orbital indicator plays crucial role in water cluster binding energies. The structures considered in this work are (H₂O)_n ($n=1-6$) whose geometries are from the WATER27 setBryantsev et al. (2009) optimized at the B3LYP/6-311++G(2d,2p) theory. The hexamer structure has a few known isomers, and we considered the book (b), cage (c), prism (p), and ring (r) isomers. The results are compared against the CCSD(T)-F12b values from Ref. [Manna et al., 2017] in Table 8. We obtained the MAEs of 118.9, 172.1, and 46.9 meV/H₂O for PZSIC, LSIC($z$), and LSIC($w$), respectively. Thus, LSIC($z$) underestimates binding energies of water cluster by roughly similar magnitude as LSDA (MAE, 183.4 meV/H₂O). This is one case where LSIC($z$) does not improve over PZSIC. A simple explanation for this behavior of LSIC($z$) is that although $z_{\sigma}$ used in LSIC($z$) can detect the weak bond regions, it ($z_{\sigma}$) cannot differentiate the slow-varying density regions from weak bond regions. The $z_{\sigma}\rightarrow 0$ in the both situations causing the weak regions to be improperly treated. Fig. 6 (a) shows $z_{\sigma}$ for water dimer where both slow-varying density and weak interaction regions are detected but not differentiated. As a result, the total energies of the complex shift too much in comparison to the individual water molecules. Thus, the underestimation of water cluster binding energies is due to the choice of $z$ and not the LSIC method. Indeed by choosing the $w$ as a scaling parameter, the binding energies are much improved. Fig. 6 (b) shows there is no discontinuity of $w$ between the two water molecules ($w_{i}$’s of two FLOs along the hydrogen bond are plotted together in the figure). Hence unlike in LSIC($z$), weak interacting region is not improperly scaled down with LSIC($w$). LSIC($w$) shows MAE of 46.9 meV/H₂O comparable to SCAN (MAE, 35.2 meV/H₂O). This result is interesting as SCAN uses a function that can identify weak bond interaction. So LSIC($w$)-LSDA may be behaving qualitatively similar to the detection function in SCAN in weak bond regions. The study of water binding energies is so far a unique case where the original LSIC($z$) performed poorly. But LSIC can be improved by simply using a different iso-orbital indicator. This case serves as a motivation in identifying appropriate iso-orbital indicator that would work for all bonding regions in LSIC.

We now provide a qualitative explanation of why LSIC($w$) gives improved results over PZSIC. This reasoning is along the same line as offered by Zope et al. Zope et al. (2019). As mentioned in Sec. I, when the self-interaction-errors are removed using PZSIC, improved description of barrier heights which involve stretch bonds is obtained but the equilibrium properties like total energies, atomization energies etc. are usually deteriorated compared to the uncorrected functional. This is especially so for the functionals that go beyond the simple LSDA. Typically this is because of over correcting tendency of PZSIC. The non-empirical semilocal DFA functionals are designed to be exact in the uniform electron gas limit, this exact condition is no longer satisfied when PZSIC is applied to the functionalsSantra and Perdew (2019). This can be seen from the exchange energies of noble gas atoms and the extrapolation using the large-$Z$ expansion of $E_{X}$ as shown in Fig. 7. Following Ref. [Santra and Perdew, 2019] we plot atomic exchange energies as a function of $z^{-1/3}$. Thus, the region near origin corresponds to the uniform gas limit. The plot was obtained by fitting the exchange exchange energies ($E_{X}$) of Ne, Ar, Kr, and Xe atoms (within LSIC($w$)-LSDA, LSIC($z$)-LSDA, and LSDA) for Ne, Ar, Kr, and Xe atoms to the curve using the following fitting functionSantra and Perdew (2019).

$$\frac{E_{X}^{approx}-E_{X}^{exact}}{E_{X}^{exact}}\times 100\%=a+bx^{2}+cx^{3},$$

(12)

where $x=Z^{-1/3}$ and a, b, and c are the fitting parameters. The LSDA is exact in the uniform gas limit. So also is LSIC($z$) since the scaling factor $z_{\sigma}$ approaches zero as the gradient of electron density vanishes while the kinetic energy density in the denominator remains finite. The small deviation from zero seen near origin (in Fig. 7) for LSIC($z$) is due to the fitting error (due to limited data point). This error is $-0.62\%$ for LSIC($z$). Thus correcting LSDA using PZSIC introduces large error in the uniform gas limit. The scaling factor $w$ used here identifies single-electron region since the density ratio approaches one in this limit. Fig. 7 shows that present LSIC($w$) approach also recovers the lost uniform limit gas. This partly explains the success of LSIC($w$). Though performance of LSIC ($w$) is substantially better than PZSIC-LSDA it falls short of LSIC($z$). On the other hand, unlike LSIC($z$) it provides good description of weak hydrogen bonds highlighting the need of identifying suitable iso-orbital indicators or scaling factor(s) to apply pointwise SIC using LSIC method. One possible choice may be scaling factor that are functions of $\alpha$ used in construction of SCAN meta-GGA and recently proposedFurness and Sun (2019) $\beta$ parameter. A scaling factor containing $\beta$ recently used by Yamamoto and coworkers with OSIC scheme showed improved resultsYamamoto et al. (2020). Future work would involve designing suitable scaling factors involving $\beta$ for use in LSIC method.