Investigation on the stability of the Enol Tautomer of Favipiravir and its derivatives by DFT, QTAIM, NBO, NLO and ¹H-NMR

Abstract

Relative stability energies of enol and keto forms of Favipiravir drug (6-fluoro-3-hydroxy-2-pyrazine carboxamide) and its twelve derivatives were theoretically envisaged using the B3LYP/6-311 + G(3df,2p)//B3LYP/6-311 + G(d,p) level of theory. Influence of the nature of the substituted groups on the geometric structures, IMHB, vibrational frequencies, dipole moments, and global reactivity descriptors was systematically studied. The enol form was found to be more stable than the keto form, regardless of the nature of the substituted group, by an average of 12.5 and 5.4 kcal/mol in the gas phase and in the aqueous solution, respectively. The role of the IMHB strength in the stability of the enol form compared with the keto form was investigated in terms of QTAIM and NBO analysis and ¹H NMR chemical shifts. The NLO properties at CAM-B3LYP with 6-311 + G(d,p) basis set of the two forms and their derivatives were also used to explain the relative stability.

KEYWORDS:

• Favipiravir
• quantum theory of atoms in molecules (QTAIM)
• intramolecular hydrogen bond (IMHB)
• substitution effect
• natural bond orbital (NBO)
• nonlinear optics (NLO)

1. Introduction

Fujifilm group in Japan has developed and demonstrated the activity of FAV T-705 drug(6-fluoro-3-hydroxy-2-pyrazine carboxamide) against a wide range of RNA viruses [1–6]. From a chemistry point of view, as FAV belongs to the 2-hydroxy pyrazine carboxamide family, it is expected to undergo tautomerization such 1,3-proton transfer process [7–9], and therefore, several tautomeric forms have been envisaged for this family. Additionally, the presence of the salicylamide structural moiety plays an important role in the excited-state proton-transfer properties [10–14]. In the last few years, tautomerization of the FAV T-705 and its analogue T-1105 (3-hydroxy-2-pyrazine carboxamide), as well as some of their derivatives, have been extensively studied in the literature in both gas phase and solution [15–20]. For example, Antonov et al. [15,16] studied the relative stability of the keto and enol forms of T-705 and T-1105 using quantum code in both the gas phase and solution, including water, toluene, and acetonitrile. Their results showed that the enol form is more stable than the keto form, regardless of the solvent. Wazzan et al. [17] used the DFT B3LYP method to envisage five different tautomers and conformers in both the gas phase and aqueous solution and they confirmed that the enol tautomer is the most stable one. Furthermore, it was also found that the stability of the keto form is increased due to the dipole–dipole interaction of the keto form with the solvent molecules [15–17].

Umar [18] investigated the tautomeric states and rotational isomers of FAV (T-705) and some of its derivatives using the DFT method, and the results showed that the enol forms are more stable than the keto forms in the range of 7.86–10.72 kcal/mol in the gas phase and 1.07–3.46 kcal/mol in the solution phase. Kavitha and Alivelu [21] used the B3LYP/6-311++G (d,p) level to calculate the resonance energies between keto and enol form by using NBO analysis. Also, they used the quantum theory of atoms in molecules (QTAIM) to calculate the IMHB energies in both forms (Keto: N … .H and enol: O … H). They found that the O..H IMHB is stronger than the N … H IMHB.

Systematic computational investigation was also used to investigate the electronic structure, spectroscopic properties, and tautomerism of some halogenated FAV derivatives (fluorine, chlorine, and bromine) by Almeida La Porta et al. [22]. They suggested that spectroscopic properties are a useful tool to elucidate such tautomeric forms.

Inter- and intra-molecular hydrogen bonds play a vital role in many chemical, physical, and biological processes [23,24]. Of particular, special attention is often paid to IMHB. According to Jeffrey's monograph [2], IMHBs appear to have been discovered very early by Sidgwick and Callow in 1924 [3]. They found that these interactions explained the differences in physical properties between ortho-substituted phenols (where they were present) and meta- and para-substituted phenols (where they were not observed). Different properties were also observed between some 1,2-disubstituted benzenes and their 1,3-benzenes as the 1,4 disubstituted counterparts. This can be explained by the formation of intramolecular rings (quasi-rings) in the former species [3].

Development of Bader’s QTAIM [25–27] is a preferable approach to studying the topological characteristics of electronic distribution within the molecule, which can be used to identify and characterize the type of interactions within the molecules; covalent, coordinate, van der Waal, … etc. [28–33]. In particular, QTAIM is used to study the IMHB, in which a bond critical point (BCP) is located between the bridge hydrogen and the hydrogen acceptor. The topological parameters of the IMHB such as electronic density ${\rho }_{\mathrm{bcp}}$ and its Laplacian ${\mathrm{\nabla }}_{\rho }^2$ at bcp of the HB and the potential energy density (V) [34,35]. E_HB follows the exponential dependence on geometrical parameters of the IMHB [36] and the E_HB is given by the following equation as follows: (1) $E_{\mathrm{HB}}=0.5\times |V|.$(1) In the above equation, 0.5 corresponds to the angular coefficient and the values of E_HB and V are in kcal/mol.

Frontier molecular orbitals (FMO) theory is a potent practical model that is used to describe the optical and electronic properties of a molecule along with its reactivity and stability [37–39]. Quantum chemical parameters such as the HOMO and LUMO play a remarkable role during molecular interaction, especially the IMHB [40,41]. HOMO is the electron-rich orbital having high energy and can donate electrons. On the other hand, LUMO is the lowest-lying empty orbital containing the capability to accept electrons.

To the best of our knowledge, the prototropic mechanism of the keto ↔ enol tautomerism process of FAV T-705 and its analogue derivatives has been reported in the literature, neither experimentally nor theoretically. Therefore, the main aim of this study focuses on the investigation of the main factors that enhance the stability of the enol form of FAV and some of its derivatives compared to the analogue keto forms using different computational tools. To achieve this goal, several issues were considered. The first issue includes a systematic theoretical investigation of the relative stability of the different tautomers of the title compound and its derivatives. Of particular interest is the influence of the nature of the substituted groups (Scheme 1) on the relative stabilities of the different forms, and the change of structural geometries such as bond lengths, bond angles, and torsion angles. The second issue provides a full QTAIM, NBO analysis, proton nuclear magnetic resonance (¹H NMR) chemical shift, and non-linear optical (NLO) properties studies.

Scheme 1. Schematic representation of keto and enol forms proposed for this study.

2. Computational details

The ground state geometries of the investigated species were optimized using the hybrid density functional B3LYP method [42–44]. This approach has been shown to yield reliable geometries for a wide variety of systems. All calculations were performed using the 6-311 + G(d,p) basis set using the Gaussian-09 series of programmes [45]. To ensure that all the stationary points of the potential energy surface correspond to local minima, the harmonic vibrational frequencies were calculated using the level of theory. In addition, the corresponding zero-point energy corrections (ZPE), which were scaled by a factor of 0.9806 [46], and the thermal correction for energies were also estimated. Final energies were calculated using the same functional (B3LYP) in conjunction with the 6-311 + G(3df,2p) basis set [47,48]. The polarizable continuum model (PCM) [49] was used to study the solvent effect for one solvent (water).

The binding characteristics and the IMHB were analysed by means of the QTAIM of Bader [25,27] using the Multiwfn code [50]. In addition, NBO analysis (NBO) [51] as implemented in G09 was used to evaluate the interactions between the donor–acceptor parts of the IMHB by performing a single point energy calculation with B3LYP/6-311 + G(d,p) level of theory at the optimized geometries. The nature of the noncovalent interactions (NCI) between the proton donor and proton acceptor pair in the E2 and K2 forms of FAV drug derivatives was analysed using the Multiwfn code [50] and Visual molecular dynamics (VMD) 1.9.4 [52,53]. Visualization of the 3D-isosurfaces maps of the FAV drug were visualized with the aid of Gunplot version 5.5-git [54] programme. ¹H NMR chemical shielding was evaluated using the GIAO (gauge-including atomic orbitals) approach [55,56] using NMR = GIAO with the B3LYP/6-311 + G(d,p) level of theory at the optimized geometries obtained at the same level of theory in gas phase and in aqueous solution (water). In order to compare with experiment, the calculated absolute shielding was transformed to chemical shifts using the tetramethylsilane (TMS) references of 31.9843 (gas) and 31.9757 (water) ppm, which were calculated at the same level of theory and δ = δ_{calc(reference)} - δ_calc.

Finally, the nonlinear optical (NLO) parameters were computed at the long-range CAM-B3LYP functional since traditional hybrid functionals such as B3LYP usually overestimate hyperpolarizabilities [57,58] combined with 6-311 + G(d,p) basis set in the gas phase and water (by applying the PCM of solvation).

3. Results and discussion

All sets of energetic data obtained in this study are given in Table SD1 and SD2 of the supplementary materials. For comparison, the relative energies, enthalpies, and Gibbs free energies of the tautomeric forms of FAV-F (T-705) using B3LYP/6-311 + G(d,p) in both gas phase and aqueous solution (between brackets) using PCM model are given in Table SD 3 of the supplementary materials. It is important to mention that the PCM is a commonly used method in computational chemistry to model solvation effects. In this model, the solvent is moulded as a polarizable continuum, rather than individual molecules, to reduce the computation costs [49,59].

3.1. Relative stability and tautomerization process

The relative stabilities of envisaged forms of FAV and its derivatives (Scheme 1) are summarized in Table 1. Graphically the relative stabilities of K1 and the transition state (TS), which connects K2 with E1 tautomer with respect to the most stable tautomer (E2) in the gas phase and aqueous solution are shown in Figure 1. Actually, as inferred by using the B3LYP/6-311 + G(3df,2p)//B3LYP/6-311 + G(d,p) level of theory and bearing in mind the effect of solvent for all substituents [60–62], our first results to be mentioned here is that E2 form is the most stable isomer, which in agreement with the previous studies [15–17]. It is also found that for all derivatives, tautomer E2 is largely more stable than K1, K2, and E1 isomers by 7.0–12.3, 8.3–15.4, 6.3–10.0 kcal/mol in the gas phase, respectively, and by 1.6–10.0, 3.0–8.1, and 6.0–7.2 kcal/mol in aqueous solution. Importantly, isomer E1, which is a rotamer for E2, is the second most stable in both media. Gas phase results show that, for all substituents except FAV-CHO and FAV-COOH, the K1 form is more stable than the K2 one by 2.1–4.9 and kcal/mol. In an aqueous solution, the stability of K2 is increased and the relative stability difference lies in the range of 0.9–2.7 kcal/mol, which might be ascribed to the dipole–dipole and hydrogen bond interactions between the solvent (water) and the different tautomers [62,63].

Figure 1. Graphical representation of the relative stabilities energies (kcal/mol) of tautomer K1 (a) and the transition state (b), which connects K2 with E1 with respect to the most stable tautomer (E2) calculated using B3LYP/6-311 + G(3df,2p)//B3LYP/6-311 + G(d,p) level in gas phase (Blue squares) and water (red circles).

Table 1. Gas phase and aqueous solution relative energies (ΔE in kcal/mol) of the investigated tautomers (K1, K2, E1, and E2) and the transition state (TS), which connects K2 and E1, using the B3LYP/6-311 + G(3df,2p)//B3LYP/6-311 + G(d,p) level of theory. The ZPE corrections are included.

	Gas phase					Aqueous solution
Species	K1	K2	TS	E1	E2	K1	K2	TS	E1	E2
FAV-H	8.5	10.7	40.5	8.4	0.0	2.0	3.4	36.8	6.2	0.0
FAV-F	11.2	14.6	42.7	8.4	0.0	4.6	7.1	38.7	6.8	0.0
FAV-Cl	10.8	13.7	42.1	8.5	0.0	4.1	6.4	38.4	6.1	0.0
FAV-Br	10.7	14.0	42.0	8.5	0.0	4.0	6.3	38.2	6.8	0.0
FAV-CN	9.9	11.5	41.1	8.8	0.0	2.5	3.5	37.6	6.5	0.0
FAV-NO2	7.4	8.3	38.1	6.3	0.0	2.0	3.0	37.4	7.1	0.0
FAV-CHO	12.1	9.8	39.4	10.0	0.0	2.3	3.2	36.6	7.2	0.0
FAV-COOH	10.1	9.9	39.6	8.6	0.0	1.6	3.2	36.6	7.1	0.0
FAV-CH3	8.6	12.1	41.0	8.3	0.0	2.7	4.9	37.2	6.0	0.0
FAV-NH2	10.5	15.4	42.9	8.2	0.0	10.0	8.1	38.4	6.6	0.0
FAV-OH	12.3	15.4	43.3	8.5	0.0	4.8	7.5	38.5	6.6	0.0
FAV-OCH3	10.8	14.1	42.5	8.2	0.0	4.2	6.8	38.4	6.6	0.0
FAV-SH	9.4	13.2	42.2	8.4	0.0	4.1	6.6	38.2	6.1	0.0

For K1 form, the parent molecule FAV-H is more stable than other derivatives except for FAV-NO₂ derivatives. In the case of K2 form, the parent FAV-H molecules is more stable than the other derivatives with the exception of FAV-NO₂, FAV-CHO, and FAV-COOH. Whereas, in the case of rotamer E1, the results reveal that FAV-H is more stable than other derivatives except for FAV-NO₂, FAV-CHO, OCH₃, and FAV-COOH. Monitoring of the results in an aqueous solution leads to the same results for both K1 and K2, while in the case of E1, the picture, to a large extent, differs. In this case, the results of E1 illustrate that FAV-H is more stable than other derivatives (except FAV-CH₃ and FAV-SH). According to these findings, it is clearly observed that there is no systematic dependence of the relative stability on the nature of the substituent group [33,64].

One of the most important questions to be addressed here is whether of K1 and E2 forms could be observed in the gas phase and/or in aqueous solution. To answer this question, the energy profiles of the corresponding tautomerization processes for FAV-H and FAV-F derivatives, calculated at B3LYP/6-311 + G(3df,2p)//B3LYP/6-311 + G(d,p), are shown in Figure 2. For FAV- derivative (Figure 5(a)), to get E2 from K1 (Figure 2(a)), it is obligatory, as a first step, to convert K1 to K2 by an internal rotation of the acetamide group, with a rotation barrier energy of 2.1 kcal/mol. The second step is to get E2 from K2 by a mandatory step of 1,3-hydrogen transfer between the N and O atoms to get E1 through the transition state (TS) followed by an internal rotation of the hydrogen atom of the enol group around the oxygen atom, with an energy barrier of 8.4 kcal/mol. The activation energy barrier of the transition state, which connects K2 and E1 is higher is located above the local minimum (E1) by ∼ 30 kcal/mol higher and above the global minimum E2 by ∼41 kcal/mol.

Figure 2. Gas phase energy profiles for keto-enol tautomerization processes in isolated neutrals of substituted FAV drug derivatives: (a) FAV-H and (b) FAV-F. Relative energies are in kcal/mol.

For FAV-F derivative (Figure 2(b)), similar to FAV-H, to get E2 from K1, it is obligatory, as a first step, to make an internal rotation of the acetamide group to convert K1 to K2 with an energy of 3.3 kcal/mol (endothermic process). The second step is a 1,3-proton transfer process between the adjacent N and O atoms followed by an internal rotation of the hydrogen atom around the oxygen atom, with a rotation barrier energy of 8.4 kcal/mol. The activation energy barrier that connects K2 and E1 is ∼ 28 kcal/mol higher than the local minimum (E1) and ∼43 kcal/mol higher than the global minimum (E2).

For all derivatives, the lowest activation energy barrier of the transition states of the K1 ↔ E2 tautomerization process, with respect to the global minimum (E2), corresponds to FAV-NO₂ derivative, while the highest one corresponds to FAV-OH derivative, see Figure 1 and Table 1. In fact, by considering the K2↔E1 tautomerization process, in which E1 is local minimum, it is found that the activation barriers for these transition states lie in the range 29.4–34.8 kcal/mol. It is also found that the lowest activation energy barrier for K1 ↔ E1 process corresponds FAV-CHO derivative, while the one belongs to FAV-OH derivatives. Considering the results K1 ↔ E1 process in an aqueous solution, FAV-CHO derivative has the lowest activation energy of 29.4 kcal/mol, while FAV-Cl derivative has the highest activation barrier with a relative energy of 32.2 kcal/mol (Table 1).

Actually, the analysis of our results led us to conclude, as shown previously [17], that tautomer E2 is predominant in the gas phase, while the possibility of the existence of K1 tautomer is increased in aqueous solution. Moreover, the proposed K1↔E1→E2 process is thermodynamically possible and the existence of these forms in both the gas phase and the aqueous solution is highly acceptable.

3.2. Molecular geometry, vibrational frequency, and dipole moment

One of the most important questions that must be addressed here is why the enol form (E2) is high stable than the keto form (K1). In order to answer this question, the two K1 and E2 forms have been systematically investigated in terms of their molecular geometry, topology, and NBO analysis.

The most relevant geometrical parameters (bond lengths, bond angles, and torsion angles) of the isomers K1 and E2 of FAV derivatives in gas phase are summarized in Table 2. The results in aqueous solution are summarized in Table SD 4 of the supplementary materials. As an example, the optimized geometrical structures of some selected derivatives (isomers E2 and K1 of the FAV-H, FAV-F, FAV-NH₂, and FAV-NO₂) are shown in Figure 3. The optimized structures for all the investigated species are available from the author upon request.

Figure 3. Optimized geometries (Left: Columns 1 and 2) and molecular graphs (Right: columns 3 and 4) of the E2 and K1 forms. Bond lengths are in A°, bond angles are in degree, and electron density at the bond critical point (red circles in au) and ring critical points (yellow circles).

Table 2. The most relevant geometrical properties (bond length/Å, bond angle (OHO (θ_E2) and NHO (θ_K1)) and dihedral angle (OHOC (ψ_E2) and COHN (ψ_K1))/°), vibrational frequencies of O-H and N-H bonds (ν_O-H and ν_N-H) in cm⁻¹, ¹H NMR chemical shifts (δ in ppm) and dipole moments (μ in Debye) for isomers K1 and E2 calculated using B3LYP/6-311 + G(d,p) level in gas phase and in aqueous solution.

	E2							K1
	d5	d6	θE2	ψE2	νO-H	δ	D	d5	d6	θK1	ψK1	νN-H	δ	D
	Gas
H	0.994	1.687	148.7	−0.02	3233	12.599	7.1	1.015	1.907	134.5	0.00	3480	9.281	8.8
F	0.993	1.692	148.3	−0.01	3248	12.531	5.2	1.015	1.916	134.5	0.02	3480	9.194	8.7
Cl	0.995	1.688	148.4	−0.10	3226	12.600	5.3	1.015	1.920	134.3	0.03	3485	9.146	8.6
Br	0.995	1.689	148.3	0.05	3229	12.621	5.4	1.014	1.922	134.2	0.01	3486	9.129	8.6
CN	0.999	1.664	149.1	0.13	3154	13.454	3.2	1.013	1.935	133.6	0.07	3503	8.918	10.4
NO2	1.000	1.661	149.1	−0.01	3135	13.775	3.8	1.013	1.951	133.2	−0.01	3510	8.830	10.5
CHO	0.999	1.663	149.2	−0.06	3151	13.512	6.0	1.013	1.932	133.6	−0.06	3502	9.022	11.7
COOH	0.998	1.669	149.1	0.03	3170	13.385	6.4	1.014	1.932	133.7	0.03	3499	9.006	10.6
CH3	0.993	1.697	148.4	−0.07	3260	12.292	7.3	1.015	1.900	134.9	−0.03	3471	9.374	8.7
NH2	0.990	1.712	147.8	0.19	3304	11.522	7.0	1.016	1.887	135.5	−0.41	3454	9.567	9.1
OH	0.992	1.704	148.0	0.08	3281	11.758	7.51	1.016	1.898	135.1	0.01	3465	9.428	10.4
OCH3	0.991	1.707	148.0	0.05	3288	11.988	8.0	1.016	1.896	135.2	−0.02	3465	9.497	10.3
SH	0.993	1.698	148.2	0.02	3255	12.189	5.8	1.015	1.909	134.7	0.01	3476	9.268	9.7

First inspection of the results indicates that a quasi-ring is created due to the formation of the O-H·××O and N-H·××O IMHB in E2 and K1 tautomers, respectively. One of the important features that we are looking for is the strength of the hydrogen bonds in both isomers. Thus, the N-H, O-H, and H·××O distances and the ∠ O-H·××O and ∠N-H·××O angles can be roughly used as good indicators to estimate the hydrogen bond strength. As known, the shorter the IMHB and the more opening one (closer to 180°) are used as good indications for the strength of the IMHB. The results in Table 2 clearly indicate that the O·××H distance in E2 and K1 forms lies in the range 1.694–1.737 Å and 1.910–1.962 Å, respectively, and in the range 1.661–1.712 Å and 1.896–1.951 Å in the aqueous solution, respectively. On average, the IMHB in E2 is shorter than that in K1 by 0.218 and 0.230 Å in both media, respectively, indicating that the IMHB in E2 is stronger than that in K1, which may be used as the first criterion for the stability of E2.

For all FAV-derivatives, the ∠N-H·××O angle in E2 lies in the range 146.8–148.0° and 147.8–149.2° in gas phase and aqueous solution, respectively. Whereas, the ∠N-H·××O angle lies in the range 132.8–135.5° and 133.2–135.5° in the gas phase and the aqueous solution, respectively. On average, the ∠O-H·××O (147.4° (gas) and 148.5 (water)) is much closer to 180° than the ∠N-H·××O (133.8 (gas) and 134.4 (water)), implying that the IMHB in E2 is stronger than that in the K1, which may be used as another justification for the stability of E2.

Let us now discuss the substitution effect on the geometries of the investigated species. For E2 tautomer, with respect to the parent molecule (FAV-H), the substitution of electron donating groups (EDGs) shortens the O-H, enlarges the O·××H IMHB, and decreases the OHO angle, while the substitution of electron-withdrawing groups (EWGs) enlarges the O-H, shortens the O·××H IMHB and increases the OHO angle. For K1 tautomer, the substitution of EDG shortens the O·××H bond, enlarges the N-H bond, and increases the N·××HO angle, while the reverse is true in the case of EWGs. Considering the E2 tautomer, it is found that the longest O·××H IMHB belongs to the FAV-CHO derivative, while the shortest one corresponds to the FAV-OCH₃ derivative. Whereas, the widest ∠OHO angle belongs to the FAV-CHO derivative, and the closest one corresponds to the FAV-NH₂ one. It is also found, in comparison with the parent FAV- (H) derivative, the hydrogen bonds in the FAV -CN, FAV -NO2, FAV -CHO, and FAV -COOH derivatives are shorter by 0.018, 0.017, 0.021, and 0.012 Å. For the halogen substitutions (F, Cl and Br), the hydrogen bonds are slightly longer by 0.009, 0.004, and 0.006 Å. In the case of the electron donating groups such as CH₃, NH₂, OH, OCH₃, and SH the hydrogen bond distance becomes shorter by 0.009, 0.021, 0.018, 0.022, and 0.010 Å, respectively (see Table 2). On the other hand, considering the K1 tautomer, the N·××H IMHB decreases by 0.021Å in FAV -NH₂ and increases by 0.031 Å in FAV-NO₂, with respect to the parent molecule FAV-H. Furthermore, the changes in OHO and NHO angles in E2 and K1 tautomer are very small and don’t exceed 0.1%.

The stretching vibration frequencies of the O-H and N-H bonds (ν_O-H and ν_N-H) are also depicted in Table 2. As can be deduced from the results in Table 2, the geometric changes caused by the hydrogen-bonding formation are well correlated with ν_O-H and ν_N-H. The results reveal that there is a red-shift of the stretching mode due to the formation of the IMHB. For E2 and K1 tautomers, upon IMHB formation, the O-H and N-H bonds (proton donating bonds) are lengthening, which is escorted by the red-shift of the frequency modes. For E2 tautomer, the greatest shifts were for EDGs (NH₂, OCH₃, OH, and CH₃, as well as F), while the smallest red-shifts were for EWGs (Cl, Br, NO₂, CHO, CN, and COOH). For K2 tautomer, the situation is different, it is found that the greatest shifts for EWGs (NO₂, CHO, CN, and COOH), while the smallest red-shifts for EDGs (NH₂, OCH₃, OH, and CH₃). Whereas the picture is vice-versed in the case of isomer K2. Our theoretical results show that the ν_O-H nicely correlates with the O-H and O·××H distances (R² = 0.999 and 0.974, respectively), similarly, ν_N-H adequately correlates with N-H and O·××H distances (R² = 0.963 and 0.956, respectively). These findings imply that vibrational frequencies can be easily evaluated from O-H and N-H and O·××H distances using the following linear equation: (2) $\begin{array}{rl}\text{E}2:{\nu }_{(\text{O}-\text{H})}& =-17701{\text{d}}_{\text{O}-\text{H}}+20849\text{and}\\ \text{K}2:{\nu }_{(\text{N}-\text{H})}& =-14517{\text{d}}_{\text{N}-\text{H}}+18215\end{array}$(2) Furthermore, it is also found that the δ_H NMR chemical shift of the donated proton (H6) in E2 tautomer adequately correlates with the ν_O–H and O–H and O·××H distances and ∠OHN angle (R² = 0.985, 0.984, 0.968, and 0.881, respectively). For K2 tautomer, similar correlations are also obtained between the δ_H NMR chemical shift of the donated proton (H6) and ν_(N-H) and N-H and O·××H distances and ∠OHN angle (R² = 0.976, 0.929, 0.988, and 0.976, respectively). These permit us to estimate the proton NMR chemical shift from the geometrical properties of the molecule using the appropriate mathematical models.

3.3. Dipole moments

The dipole moment values of K1 and E2 forms in both media are collected in Table 2. For the sake of comparison, the graphical representation of the dipole moments of both forms in the gas phase (G) and the aqueous solution (W) are shown in Figure 4. First inspection that should be noticed is that the dipole moments of FAV-derivatives values in aqueous solution are higher than that in the gas phase, which can be explained in terms of the high dielectric constant of water. Our theoretical results show, for all substitutions and in both media, that the dipole moment values of K1 tautomer are higher than those of E2 one, implying that K1 is more polar than E2 tautomer. The gas phase dipole moment values of E2 and K1 tautomers lie in the range 2.8–6.1 and 4.6–8.1 Debye, respectively, while their values in aqueous solution lie in the range 3.5–8.5 Debye (E2) and 8.6–11.7 Debye, respectively. For E2 form, the electron-withdrawing substitution FAV-CN has the smallest dipole moment (gas: 2.8 and water: 3.5 Debye), while the electron-donating substitution FAV-NH₂ has the largest dipole moment (gas: 7.7 and water: 8.5 Debye). The picture is reversed in K1 form, the largest dipole moment corresponds to the electron-withdrawing substitution FAV-CHO (gas: 8.1 and water: 11.7 Debye), while the smallest dipole moment belongs to the electron-donating substitution FAV-OCH3 (4.6 and 10.3). These results may be explained in terms of consideration of charge values on the most polar atoms (O and N). For E2, it can be stated that in the case of EDGs (FAV-NH₂), N (−0.0889, −0.1382, −0.4273 e) and O (−0.4315, −0.2648) atoms carry the most negative charge but in case of FAV-CN they carry the least negative charge (N: −0.0642, −0.140, −0.4084 and O: −0.4255, −0.2408). The reverse observations are also noticed in the case of K1 form. On the other hand, as expected, the polarity of both forms is increased in aqueous solution and the dipole moment values of both forms become higher than those in the gas phase, see Table 2, which might be attributed to the decrease in the angle between the dipole vectors of the lone pairs on oxygen atoms as the number of hydrogen bonds to that oxygen increases, which perturbs the structure and induced a dipole moment in the molecule.

Figure 4. Graphical representation of the dipole moment (D in Debye) of (a) E2 and (b) K2 of FAV-derivatives calculated using B3LYP/6-311 + G(d,p), blue column (gas) and orange column (aqueous solution). For sake of comparison, the scales of y-axis in the two panels are unified in both panels.

3.4. Frontier molecular orbitals and chemical reactivity

The distribution of the electron density of the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) of some Fav-derivatives, including FAV-H, FAVAV-F. FAV-NH2, FAV-CH3, FAV.OH and FAV-OCH₃, are shown in Figure 5. The electron density of HOMO and LUMO of the other FAV derivatives are shown in Figure SD 1 of the supplementary materials. The results obtained show that the electron distribution of the HOMO and LUMO occupies more than 90% for both forms of FAV drug moiety and its derivatives are concentrated over the whole moiety of the drug. The HOMO regions are characterized by electron-rich centres, while the LUMO regions are characterized by electron-poor centres. Therefore, our results can indicate that the HOMO orbital own to the π-electron delocalization and the lone pairs of electrons over the nitrogen and oxygen atoms, and the LUMO orbital is concentrated over the regions which have the π* character. These results signify that the first excited state might be attributed to the $\pi \to {\pi }^{\ast }$ transition [37–41,65]. It is important to notice that increases in the electron density on the oxygen atom (O1) by substitution of the EDGs increase the strength of the IMHB, while the reverse is true in the case of substitution of the EWGs [41]. On the other hand, it is also noticed that the electronic densities of the hydroxyl oxygen atom in the case of E2 form decrease from HOMO to LUMO orbital, while the electronic densities of the amino nitrogen atom in the case of K1 form increase from HOMO to LUMO. This means that the interaction between O and H atoms is stronger than that between N and H atoms.

Figure 5. The HOMO and LUMO of some FAV derivatives, including FAV-H, FAV-F, F-NH2, FAV-CH₃ and FAV-OCH₃.

Quantum chemical descriptors, including the energies of HOMO and LUMO (E_HOMO and E_LUMO), the energy gap (ΔE_gap= E_LUMO-E_HOMO), global hardness (η = (E_LUMO-E_HOMO)/2), chemical potential (μ =− (E_LUMO+ E_HOMO)/2) and electrophilicity (ω = μ²/2η) of the FAV drug and its derivatives in the gas phase are summarized in Table 3, while the those in aqueous solution are listed in Table SD 5. As a principal consequence, strong electron-acceptor propensity leads to great stabilization of the LUMO, which lowers the energy gap of the chemical species. Previous studies pointed out that a molecule with a low energy gap is more polarizable and is generally associated with high chemical reactivity and low kinetic stability and is termed a soft molecule [66–69]. Lowering of the energy gap guides to increase the reactivity of the molecule. First inspection of our results shows that the parent molecule FAV-H in its two forms (E2 and K1) has the highest energy gap (except FAV-COOH in the case of E2), indicating that the stability of these derivatives increases with respect to the parent molecule, regardless of the nature of the substituents.

Table 3. Frontier molecular orbital energy (E_HOMO and E_LUMO in eV), energy gap (ΔE in eV), chemical hardness (η in eV), softness (σ in ev⁻¹) electrophilicity (ω in eV) and chemical potential (μ in eV) for FAV and its derivatives calculated using the B3LYP/6-311 + G(d,p) level of theory in the gas phase.

	E2							K1
	EHOMO	ELUMO	ΔE	μ	η	σ	ω	EHOMO	ELUMO	ΔE	μ	η	σ	ω
FAV-H	−7.307	−2.515	4.793	−4.911	1.198	0.835	10.066	−7.131	−2.764	4.367	−4.948	1.092	0.916	11.210
FAV-F	−7.367	−2.849	4.518	−5.108	1.130	0.885	11.550	−7.203	−3.178	4.025	−5.191	1.006	0.994	13.387
FAV-Cl	−7.303	−2.823	4.480	−5.063	1.120	0.893	11.444	−7.146	−3.124	4.022	−5.135	1.006	0.995	13.112
FAV-Br	−7.239	−2.824	4.416	−5.032	1.104	0.906	11.468	−7.083	−3.117	3.966	−5.100	0.992	1.009	13.116
FAV-CN	−7.807	−3.135	4.671	−5.471	1.168	0.856	12.813	−7.577	−3.367	4.210	−5.472	1.053	0.950	14.225
FAV-NO2	−8.130	−3.437	4.693	−5.784	1.173	0.852	14.255	−7.644	−3.506	4.138	−5.575	1.035	0.967	15.022
FAV-CHO	−7.622	−2.934	4.688	−5.278	1.172	0.853	11.885	−7.369	−3.123	4.246	−5.246	1.062	0.942	12.963
FAV-COOH	−7.649	−2.816	4.834	−5.233	1.208	0.828	11.330	−7.291	−3.007	4.283	−5.149	1.071	0.934	12.377
FAV-CH3	−7.027	−2.428	4.598	−4.728	1.150	0.870	9.719	−6.886	−2.668	4.218	−4.777	1.055	0.948	10.820
FAV-NH2	−6.281	−2.339	3.942	−4.310	0.986	1.015	9.425	−6.173	−2.691	3.481	−4.432	0.871	1.149	11.282
FAV-OH	−6.281	−2.339	3.942	−4.310	0.986	1.015	9.425	−6.697	−2.878	3.819	−4.788	0.955	1.047	12.003
FAV-OCH3	−6.665	−2.464	4.201	−4.565	1.050	0.952	9.919	−6.586	−2.751	3.835	−4.669	0.959	1.043	11.366
FAV-SH	−6.660	−2.628	4.032	−4.644	1.008	0.992	10.698	−6.501	−2.939	3.562	−4.720	0.891	1.123	12.509

Let us now compare the stability and the chemical reactivity of the FAV derivatives depending on the nature of the substituted group. It is clearly obvious that the EDGs (NH₂, OH, OCH₃, SH) decrease the stability and increase the chemical reactivity of the E2 form of the FAV-derivatives by lowering their energy gaps, which might be ascribed to the increase of the E_LUMO. On the other hand, the EWGs (NO₂, CHO, and CN) decrease the LUMO energies and consequently increase the energy gaps, which destabilize these derivatives. It is also found that the energy gap values lie in the range 3.942–4.834 eV in gas phase and 3.481–4.367 eV in the aqueous solution. For both forms, the lowest energy gap is found for FAV-NH₂, while the highest one corresponds to E2/FAV-COOH (E2) and K1/FAV-H (K1). These results might be explained in terms of the participation of the EDGs or EWGs in enhancing or deducing the conjugation effect in the composition of the LUMO, which decreases or increases the energy level. Similar results are also found, with slight changes in the stability order, for K1 form of FAV-derivatives in an aqueous solution.

It is worth mentioning, for all substitutions, that the energy gaps of E2 tautomer in aqueous solution are lower than those in gas phase (except FAV-NO₂, FAV-NH₂, and FAV-SH), reflecting the less kinetic stability and the higher chemical reactivity of E2 in gas phase compared to those in aqueous solution. On the other hand, the energy gaps of K2 form are lowered upon going from the gas phase to the aqueous solution (except FAV-CN, FAV-NO2, FAV-CHO, and FAV-COOH). These results indicate that the substitution of EWGs decreases the stability and increases the activity of K2 form in an aqueous solution. The exception here is the halogen which has a contrasting behaviour because it acts as a withdrawing group by inductive effect and an electron donating group by resonance.

3.5. Quantum theory of atoms in molecules (QTAIM) analysis

In this work, topological parameters as calculated using the quantum theory of atoms in molecules (QTAIM) are used to get a deeper understanding of the nature of bonds such as IMHBs [33]. The most relevant calculated topological parameters of E2 and K1 forms of the FAV and its derivatives are listed in Table 4. The other topological parameters are given in Table SD 6 of the supplementary materials.

Table 4. The most relevant topological parameters of E2 and K1 tautomers and their derivatives (/au) in gas phase, Density of all electrons at the bond critical point and its Laplacian of O–H and N–H bonds (${\rho }_{O-H}^{\mathrm{BCP}}$ and ${\mathrm{\nabla }}^2{\rho }_{O-H}^{\mathrm{BCP}}$), Density of all electrons at the bond critical point and its Laplacian of O·H hydrogen bonds (${\rho }_{O\cdots H}^{\mathrm{BCP}}$ and ${\mathrm{\nabla }}^2{\rho }_{O\cdots H}^{\mathrm{BCP}}$), hydrogen bond strength ($E_{\mathrm{HB}}^{\ast }$ in kcal/mol), and bond order of O···H hydrogen bond ($B{O}_{O\cdots H}$ (unitless)) in gas Phase. The other topological parameters are given in Table SD 6 of the supplementary materials.

	${\rho }_{\mathrm{quasi}}^{\mathrm{RCP}}$	${\mathrm{\nabla }}_{{\rho }_{\mathrm{quasi}}^{\mathrm{RCP}}}^2$	${\rho }_{O-H}^{\mathrm{BCP}}$	${\mathrm{\nabla }}_{{\rho }_{O-H}^{\mathrm{BCP}}}^2$	${\rho }_{O\cdot H}^{\mathrm{BCP}}$	${\mathrm{\nabla }}_{{\rho }_{\mathrm{OH}}^{\mathrm{BCP}}}^2$	$B{O}_{O\cdot H}$	$E_{\mathrm{HB}}^{\ast }$
FAV-H	0.0258	0.1857	0.3323	−2.3319	0.0465	0.1352	0.89	13.7
FAV-F	0.0262	0.1859	0.3328	−2.3384	0.0455	0.1337	0.89	13.3
FAV-Cl	0.0261	0.1860	0.3318	−2.3304	0.0460	0.1341	0.89	13.5
FAV-Br	0.0261	0.1864	0.3318	−2.3300	0.0459	0.1337	0.89	13.4
FAV-CN	0.0257	0.1834	0.3282	−2.2988	0.0485	0.1371	0.87	14.5
FAV-NO2	0.0264	0.1863	0.3274	−2.2924	0.0484	0.1361	0.87	14.4
FAV-CHO	0.0258	0.1845	0.3284	−2.2984	0.0488	0.1377	0.87	14.6
FAV-COOH	0.0259	0.1851	0.3292	−2.3056	0.0478	0.1360	0.88	14.2
FAV-CH3	0.0258	0.1852	0.3338	−2.3446	0.0455	0.1339	0.90	13.3
FAV-NH2	0.0256	0.1831	0.3357	−2.3607	0.0442	0.1325	0.91	12.8
FAV-OH	0.0260	0.1845	0.3344	−2.3513	0.0445	0.1323	0.90	12.9
FAV-OCH3	0.0259	0.1837	0.3353	−2.3581	0.0441	0.1318	0.91	12.8
FAV-SH	0.0258	0.1851	0.3332	−2.3410	0.0454	0.1335	0.90	13.2
	${\rho }_{\mathrm{quasi}}^{\mathrm{RCP}}$	${\mathrm{\nabla }}_{{\rho }_{\mathrm{quasi}}^{\mathrm{RCP}}}^2$	${\rho }_{N-H}^{\mathrm{BCP}}$	${\mathrm{\nabla }}_{{\rho }_{N-H}^{\mathrm{BCP}}}^2$	${\rho }_{N\cdot H}^{\mathrm{BCP}}$	${\mathrm{\nabla }}_{{\rho }_{O\cdot H}^{\mathrm{BCP}}}^2$	$B{O}_{N\cdot H}$	$E_{\mathrm{HB}}^{\ast }$
FAV-H	0.0228	0.1691	0.3334	−1.7324	0.0284	0.1069	1.11	6.93
FAV-F	0.0232	0.1703	0.3330	−1.7310	0.0283	0.1063	1.10	6.92
FAV-Cl	0.0232	0.1704	0.3334	−1.7328	0.0280	0.1053	1.11	6.78
FAV-Br	0.0232	0.1710	0.3334	−1.7326	0.0278	0.1047	1.11	6.73
FAV-CN	0.0227	0.1668	0.3345	−1.7348	0.0270	0.1027	1.12	6.47
FAV-NO2	0.0233	0.1702	0.3347	−1.7338	0.0264	0.1008	1.13	6.3
FAV-CHO	0.0229	0.1677	0.3346	−1.7337	0.0272	0.1038	1.12	6.57
FAV-COOH	0.0230	0.1687	0.3344	−1.7323	0.0273	0.1038	1.12	6.59
FAV-CH3	0.0229	0.1692	0.3327	−1.7312	0.0290	0.1085	1.10	7.12
FAV-NH2	0.0229	0.1688	0.3316	−1.7301	0.0298	0.1104	1.09	7.38
FAV-OH	0.0230	0.1689	0.3324	−1.7319	0.0293	0.1093	1.09	7.21
FAV-OCH3	0.0231	0.1700	0.3318	−1.7298	0.0292	0.1083	1.09	7.16
FAV-SH	0.0230	0.1698	0.3327	−1.7322	0.0286	0.1069	1.10	6.97

One of the most important parameters that can be calculated from the topological parameters is the IMHB strength in both forms ($E_{\mathrm{HB}}^{\ast }(E2)\text{and}{E}_{\mathrm{HB}}^{\ast }$(K1)), which has been calculated by using the Espinosa method [36]. As can be clearly obvious in Table 4, for all substituents, the $E_{\mathrm{HB}}^{\ast }(E2)$ is much larger than the $E_{\mathrm{HB}}^{\ast }$(K1). Indeed, it is found that the $E_{\mathrm{HB}}^{\ast }(E2)$ lies in the range of 12.76–14.63 kcal/mol, while the $E_{\mathrm{HB}}^{\ast }$(K1) lie in the range of 6.30–7.38 kcal/mol.

For E2 tautomer, considering the parent molecules (FAV-H) as a reference, our results indicate that the greatest $E_{\mathrm{HB}}^{\ast }(\text{E}2)$ is observed for the strongest EWGs in the molecules FAV-CHO, FAV-CN, FAV-COOH and FAV-NO₂, while the smallest ones are observed for the EDGs in the molecules FAV-OCH₃, FAV-NH₂, FAV-OH and FAV-SH. These results indicate that the substitution of the EDG weakens the hydrogen bond strength of the E2, while the EWGs strengthen these bonds. The situation is reversed in the case of K1 tautomer in which EDGs strengthen the IMHB, while the EWGs weaken the strength of the IMHB, which is consistent with the changes in the geometrical properties. This contradiction in the effect of the EDGs and EWGs can be well understood by monitoring the chemical structures of both forms. The resonance effects in E2, due to the presence of the pyrazine moiety is higher than that in the K1 tautomer, which contains dihydopyrazine moiety, instead. This is also can be observed by looking at the values of electron densities at the ring critical point (ρ_RCP) of the quasi-rings of pyrazine and dihydopyrazine moieties in both forms, respectively. The results show, for all substituents, that the ρ_RCP values of the pyrazine ring in E2 tautomer (∼ 0.026 au) is larger than that of the dihydopyrazine ring in K1 tautomer (0.023 au). The R group is substituted at the para position with respect to the carbonyl group (C = O) and the C-OH in K1 and E2 forms, respectively. Consequently, in case of E2 form, substitution of an EWG reduces the partial negative charge of the O5 atom and increases the partial positive charge of the H6 atom with respect to the parent molecule, whereas substitution of an EDG is expected to show the opposite effect. On the other hand, the truth is revered in the case of K1 form [70,71].

The halogenated derivative due to the inductive and resonance effect behaves as EDGs. Actually, our results show that the $E_{\mathrm{HB}}^{\ast }(K1)$ values of FAV-F, FAV-Cl and FAV-Br derivatives are 6.92, 6.78 and 6.73 kcal/mol, respectively, whereas, the $E_{\mathrm{HB}}^{\ast }(E2)$ values for the same substitutions are 13.31, 13.50, and 13.43 kcal/mol. These values are lower than the hydrogen bond strength of the parent molecule (FAV-H) in both forms (K1: 6.93 kcal/mol and E2: 13.69 kcal/mol).

In summary, our results deduce that the strongest hydrogen bonds in the case of E2 form occur for EWGs while the weakest ones for EDGs. On the contrary, the weakest hydrogen bonds in the case of K1 tautomer occur for EWGs while the strongest ones for EDGs. Our theoretical calculation shows that the changes in the geometries of the hydrogen bonds are well correlated with the IMHB strength. In fact, we have found that the $E_{\mathrm{HB}}^{\ast }(E2)$ are linearly correlated with the O·××H in both tautomers (determination coefficient $R^2$ for $E2=0.9995$ and for $K1=0.998$), implying that the IMHB length can be used as a good estimator for the evaluation of the hydrogen bond strength.

For E2 tautomer, Table 3 shows that the highest electron density of the O·××H hydrogen bond at BCP corresponds to CHO substitution (EWG), while the least one belongs to the OCH₃ group (EDG), whereas the maximum electron density of O-H at BCP corresponds to FAV-NH₂ molecule and the minimum one belongs to FAV-CN derivative. For K2 form, the maximum electron density of O·××H contact at BCP is found for FAV-NH₂ derivative, while the least one belongs to NO₂ substitution (EWG), whereas the reverse is true in the case of the N-H bond. For all substitutions, the electron density O·××H contact at BCP in E2 tautomer is much larger than that in K1 tautomer. The electron density of O-H is slightly larger than that of the N-H bond, with the exception of FAV-CH₃, FAV-NH₂, FAV-OH, FAV-OCH₃, and FAV-SH substitutions.

As reported in the literature [70,72–74], the nature of the bond can be identified by the bond order (BO = -G/V) ratio, for BO >1 the bond is noncovalent, while for 0.5 < -G/V < 1, it is partly covalent. Our theoretical calculation shows, for E2 tautomer, that the O·××H (IMHB) has low electron density (ranging from 0.0441 - 0.0488 au) and positive ${\mathrm{\nabla }}_{{\rho }_{\text{OH}}^{\mathrm{BCP}}}^2$values (ranging from 0.1318 - 0.1377 au) with BO (ranging from 0.87 - 0.91) and negative H_(BCP) values (ranging from −0.0039–0.0061 au). On the other hand, for K1 tautomer, it is found that the O·××H (IMHB) has low electron density (ranging from 0.0264 - 0.0298) and positive ${\mathrm{\nabla }}_{{\rho }_{\text{OH}}^{\mathrm{BCP}}}^2$values (ranging from 0.1008 - 0.1104 au) with BO (ranging from 1.09–1.13) and positive H_(BCP) values (ranging from 0.0021 - 0.0025 au). These values strongly confirm that the IMHB in E2 tautomer is partly covalent, while that in K2 tautomer is noncovalent interaction.

Furthermore, our results show that the pseudo-ring containing O·××H IMHB is created in both tautomers and thus the RCP exists. The greater the electron density at RCP belongs to the stronger IMHB and vice versa [75,76]. Inspection of the results in Table 3 indicates that the electron density at RCP for E2 tautomer (ranging from 0.0183–0.0190 au) is larger than that in K2 tautomer (ranging from 0.0132–0.0139 au), implying that, for all substitutions, the IMHB in E2 form is stronger than that in K1 one. These values can be used to verify that E2 tautomer is energetically more stable than K1 tautomer.

Interestingly, a series of linear correlations are found between the geometrical and topological parameters. For example, the ρ and distance relationship for O·××H contact is found; R² is very close to unity (0.9997). In addition, excellent correlations between ${\mathrm{\nabla }}_{{\rho }_{\text{peudo}}^{\mathrm{RCP}}}^2$ versus $E_{\mathrm{HB}}^{\ast }$, and ${\rho }_{\text{O}-\text{H}}^{\mathrm{BCP}}$; R²=  0.993 and 0.994, respectively. These results suggest that the properties of the ring critical point values might be very useful in estimating the strength of the IMHB, permitting us to estimate and predict other properties of the systems.

3.6. NBO analysis

The NBO occupation number for ${\sigma }_{\text{O}-\text{H}}^{\ast }$ antibond, oxygen lone pair of proton acceptor and second-order perturbation stabilization energies, $E^{(2)}$, which belongs to charge transfer between oxygen lone pair and ${\sigma }_{\text{O}-\text{H}}^{\ast }(\text{E2})$ and ${\sigma }_{\text{N}-\text{H}}^{\ast }$ (K2) antibonds $E_{\mathrm{LP}(\text{O})\to {\sigma }_{\text{O}-\text{H}}^{\ast }}^{(2)}$(E2) and $E_{\mathrm{LP}(\text{O})\to {\sigma }_{\text{N}-\text{H}}^{\ast }}^{(2)}$ (K2) in kcal/mol ($E_{\mathrm{HB}}^{\mathrm{CT}}$) are presented in Table 5. In the current study, the most important feature in the NBO analysis of the hydrogen-bond system is the charge transfer between the lone pairs of the proton acceptor (LP(O) atom) and the proton donor (${\sigma }_{\text{O}-\text{H}}^{\ast }$ and ${\sigma }_{\text{N}-\text{H}}^{\ast }$ antibonds). Upon the formation of the hydrogen bond, in both forms, the occupancies of the ${\sigma }_{\text{O}-\text{H}}^{\ast }$ and ${\sigma }_{\text{N}-\text{H}}^{\ast }$ antibonds increase and further enlarge and weaken the O-H bond in E2 tautomer and N-H bond in K1 tautomer. Our results show that the gas phase calculated values of the $E_{\mathrm{HB}}^{\mathrm{CT}}$ values for E2 and K2 tautomers lie in the range of 19.66–23.73 and 6.55–8.72 kcal/mol, respectively (Table 5). It is also found that, for all substituents, the $E_{\mathrm{HB}}^{\mathrm{CT}}$ values in E2 tautomer are much stronger than the corresponding $E_{\mathrm{HB}}^{\mathrm{CT}}$ values in K1 form, indicating that the O·××H in E2 is stronger than that in K1 form, which benefits the stabilization of E2 with respect to the K2 tautomer, confirming the results obtained by the geometrical structures (bond lengths and bond angles), QTAIM and the relative energies.

Table 5. Charge transfer energy of the IMHB $(E_{\mathrm{HB}}^{\mathrm{CT}})$, occupation numbers (O.N) and p-character in ${\sigma }_{O-H}\text{and}{\sigma }_{N-H}$ bonds calculated using B3LYP/6-311 + G(d,p) level in gas phase at 298.15 K.

	E2				K1
				p-character				p-character
R	^a$E_{\mathrm{HB}}^{\mathrm{CT}}$	$O.N_{(O)}$	$O.N_{{\sigma }_{O-H}^{\ast }}$	in ${\sigma }_{O-H}$	$E_{\mathrm{LP}(O)\to {\sigma }_{N-H}^{\ast }}^{(2)}$	$O.N_{(O)}$	$O.N_{{\sigma }_{N-H}^{\ast }}$	in ${\sigma }_{N-H}$
FAV-H	21.69	1.857	0.054	3.21	7.72	1.858	0.025	2.06
FAV-F	20.96	1.856	0.052	3.22	7.68	1.855	0.025	2.04
FAV-Cl	21.40	1.856	0.053	3.22	7.51	1.856	0.025	2.07
FAV-Br	21.31	1.856	0.053	3.22	7.43	1.856	0.025	2.07
FAV-CN	23.58	1.853	0.057	3.21	6.84	1.856	0.023	2.07
FAV-NO2	21.42	1.853	0.058	3.21	6.55	1.856	0.023	2.08
FAV-CHO	23.73	1.853	0.058	3.21	6.93	1.858	0.024	2.07
FAV-COOH	22.99	1.854	0.056	3.21	7.02	1.858	0.024	2.07
FAV-CH3	20.80	1.858	0.052	3.22	8.11	1.859	0.026	2.06
FAV-NH2	19.66	1.859	0.049	3.23	8.72	1.858	0.028	2.06
FAV-OH	20.07	1.858	0.050	3.22	8.31	1.857	0.027	2.06
FAV-OCH3	19.66	1.859	0.049	3.22	8.36	1.857	0.027	2.06
FAV-SH	20.82	1.973	0.052	3.22	7.93	1.857	0.026	2.06

^a Charge transfer energy between Oxygen lone pair and ${\sigma }_{O-H}^{\ast }$ and ${\sigma }_{N-H}^{\ast }$ antibonds ($E_{\mathrm{HB}}^{\mathrm{CT}}$), $E_{\mathrm{LP}(O)\to {\sigma }_{O-H}^{\ast }}^{(2)}$ and $E_{\mathrm{LP}(O)\to {\sigma }_{N-H}^{\ast }}^{(2)}$.

^b Occupation numbers (O.N) for oxygen lone pair and each of ${\sigma }_{O-H}^{\ast }$ and ${\sigma }_{N-H}^{\ast }$ antibonds.

In fact, it is found that the smallest $E_{\mathrm{HB}}^{\mathrm{CT}}$(E2) values correspond to FAV-OCH₃, FAV-OH, FAV-CH₃, and FAV-SH derivatives, respectively (19.66, 20.07, 20.80 and 20.82 kcal/mol, respectively), while the highest values correspond to FAV-CHO, FAV-CN and FAV-COOH derivatives, respectively (23.73, 23.58 and 22.99 kcal/mol, respectively). On the other hand, the smallest $E_{\mathrm{HB}}^{\mathrm{CT}}$(K1) values correspond to FAV-NO₂, FAV-CN, and FAV-CHO derivatives, respectively (6.55, 6.84, and 3.93 kcal/mol, respectively), while the highest ones correspond to FAV-NH₂, FAV-OCH₃ and FAV-OH derivatives, respectively (8.72, 8.36 and 8.31 kcal/mol, respectively). These results are also confirmed which is in agreement with the EHB as calculated by the QTAIM. Because of the induction effect the highest one belongs to FAV -CHO (23.73 kcal/mol). On the other hand, the least $E_{\mathrm{HB}}^{\mathrm{CT}}$ (K1) value is found for FAV -NO₂ (6.55 kcal/mol), while the highest one is found for FAV -NH₂ (8.72 kcal/mol). For both tautomers, the summarized data in Table 5 show that the $E_{\mathrm{HB}}^{\mathrm{CT}}$ values for F substitution are higher than those of Cl and Br substitutions, respectively, which might be attributed to the low inductive and field effects of Cl and Br compared to F atom. These results might be attributed to the resonance effects, which reduce the electron-withdrawing strength of the inductive and field effects of the F- group more than in the case of the Cl- and Br- groups. Therefore, the Cl and Br atoms become, as overall, more electron-withdrawing than F atom [76,77]. On the other hand, it is also found that as the electron-withdrawing power increases (CN, NO₂, and CHO), the $E_{\mathrm{HB}}^{\mathrm{CT}}$ interaction increases. Similar to the electron withdrawing group effect, the substitution of the electron-donating groups is expected to reduce the interaction between the proton donor and proton acceptor pair, making the $E_{\mathrm{HB}}^{\mathrm{CT}}$ smaller (see Table 5).

As can be seen in Table 5, the occupancy of the ${\sigma }_{O-H}^{\ast }$ (E2) is much higher than the occupancy of the ${\sigma }_{N-H}^{\ast }$ (K1), and they lie in the range of 0.049–0.058 e and 0.023–0.028 e, respectively. It is also found that the changes of oxidation numbers of oxygen atoms in E2 and K1 tautomers and ${\sigma }_{O-H}^{\ast }$ (E2 tautomer) and ${\sigma }_{N-H}^{\ast }$ (K1 tautomer) are generally in agreement with the $E_{\mathrm{LP}(O)\to {\sigma }_{O-H}^{\ast }}^{(2)}$(E2) and $E_{\mathrm{LP}(O)\to {\sigma }_{N-H}^{\ast }}^{(2)}$ (K2) charge transfer from LP(S) to σ*SH and hydrogen-bond-formation energy. The maximum and minimum occupancy of the ${\sigma }_{O-H}^{\ast }$ (E2) are for R = NH₂ and R = CHO, respectively, while those for the ${\sigma }_{N-H}^{\ast }$ (K1) are for R = NO₂ and NH₂.

Furthermore, the natural hybrid orbitals (NHOs) can be used to describe the bonding as implemented in the NBO analysis. The atomic charge distribution and percentage of the s-character of O1 in E2 and K1 tautomers and the p-character of O5 (E2 tautomer) and N5 (K1 tautomer) atoms are investigated. The p-character of O5 (E2 tautomer) and N5 (K1 tautomer) atoms natural hybrid orbitals of ${\sigma }_{O-H}$ (E2) and ${\sigma }_{N-H}$ (K1) are also presented in Table 5. The first inspection shows that the p-character of O5 atoms natural hybrid orbitals of ${\sigma }_{O-H}$ (E2) is much larger than the p-character of N5 atoms natural hybrid orbitals of ${\sigma }_{N-H}$ (K1), indicating that the O-H bond is shorter than the N-H bond [78]. For E2 tautomer, the results obtained show that the p-character of O5 natural hybrid orbitals of ${\sigma }_{O-H}$ lie in the range of 3.11–3.23. Whereas, the p-character of N5 atom natural hybrid orbitals of ${\sigma }_{N-H}$ are much smaller than those of E2 tautomer and they lie in the range of 2.04–2.08. For E2 tautomer, the p-character of O5 natural hybrid orbitals of ${\sigma }_{O-H}$ in F substitution (sp^3.22) is slightly larger than that in the parent molecule (H-substitution, sp^3.21), similarly, the O-­H bond length in F substitution (0.990 Å) is slightly shorter than the usual O-H bond length (d₅) in the parent molecule of 0.991 Å. Whereas, it is much smaller in NO₂ substitution (sp^3.11). These results might be explained in terms of the nature of the substitution group. A similar conclusion can be also obtained when the K2 tautomer is considered (see Table 5). These findings enhance the conclusion that is raised by the other results, which prove that E2 tautomer is more stable than K1 tautomer.

Of interest, it is shown that the $E_{\mathrm{HB}}^{\mathrm{CT}}$ energy adequately correlates with other geometric and topological parameters. For example, it is found that there is a good correlation between $E_{\mathrm{HB}}^{\mathrm{CT}}$ (E2) versus O-H (R² = 0.9821), O·××H (R² = 0.9928), ${\rho }_{\text{OH}}$ (R² = 0.9933) and $E_{\mathrm{HB}}^V$ (R² = 0.9842). Similarly, it is also found nicely correlations between $E_{\mathrm{HB}}^{\mathrm{CT}}$ (K1) versus N-H, O·××H, ${\rho }_{\text{O}\cdot \text{H}}$ and $E_{\mathrm{HB}}^V$; correlation coefficients are in the range 0.998–0.994. These findings suggest that the properties of the charge transfer between the lone pairs of proton acceptor and antibonds of proton donor could be very useful in estimating the strength of the IMHB.

3.7. Non-covalent interactions (NCI) analysis

The NCI analysis could be very useful in understanding the constructive interactions between the lone pairs of proton acceptors and antibonds of proton donors within a molecule (IMHB). NCI plot is useful in illustrating the correlation of the strong directional attractions with the localized atom-atom contacts and the molecular division having weak interactions. The reduced density gradient (RDG) is an incentive non-dimensional quantity, which involves the density and first derivative, and it is given as: (3) $\text{RDG}(r)=\frac{1}{2{(3{\pi }^2)}^{1/3}}\frac{|{\mathrm{\nabla }}_{\rho }(r)|}{\rho {(r)}^{4/3}}$(3) where ρ(r) represents the electron density, $|{\mathrm{\nabla }}_{\rho }(r)|$ is it electron density gradient strength and sign λ₂ presents the second eigen value of the electron density Hessian matrix and represents the type of bond. In the current study, four different FAV-derivatives for both isomers E2 and K1 were presented, namely: FAV-H, FAV-F, FAV-OCH₃, and FAV-NO₂, representing the different substitution effects (Figure 5). Scatter graphs of RDG in Figure 6 (left) indicate the weak forces of interaction, which occur between the atoms in the molecule. More information on NCI can be obtained from the NCI index (Top of Figure 6), which is built upon the RDG. Statistical information about the nature and power of interactions of molecules can be obtained from the ρ quantity of the RDG against sign (λ2)ρ peaks. Of particular importance in the nature of the interaction is the value of sign (λ2)ρ; (λ2)ρ > 0 indicates a repulsive interaction (non-bonded) and the sign (λ2)ρ < 0 shows an attractive interaction (bonded) and sign (λ2)ρ nearly zero for van der Waals interaction (weak interaction). Multiwfn and VMD software were used to study the interaction of the power in the molecular system, which designates the stronger advisability of the blue colour and the push of red.

Figure 6. Maps of reduced density gradient (RDG) versus the electron density ρ multiplied by the sign of λ₂ (left) and 3D-Isosurface of colour scaling of weak interactions (NIC) of FAV derivatives (right) (FAV-H, FAV-F, FAV-OCH3 and FAV-NO₂).

Inspection of the 3D-Isosurface of colour scaling of weak interactions diagram (NIC, Figure 6 (right)) indicates that the interactions between the proton donor and proton acceptor (O atom) in tautomer E2 is stronger than that in K1, reflecting that the amount of the strength of the IMHB. Among the explored Fav-derivatives, the FAV-NO₂ shows the weakest NCI, while the FAV-OCH₃ shows the strongest NCI.

3.8. NLO properties

Previous studies showed that non-linear optical (NLO) compounds exhibited extraordinary progressions in many areas such as electro-optics [79], fibre optics [80,81], data transformation, data storage in the field of wireless communication and photonic laser [79]. These materials are also stimulating and the most important subject for researchers to design high-performance NLO compounds using organic and inorganic systems. Two important NLO parameters were calculated. The anisotropic polarizability ($\mathrm{\Delta }{\alpha }_o$) which is defined as [82]: (1) $\begin{array}{rl}\mathrm{\Delta }{\alpha }_o& =\frac{1}{\sqrt{2}}[{({\alpha }_{\mathrm{xx}}-{\alpha }_{\mathrm{yy}})}^2+{({\alpha }_{\mathrm{yy}}-{\alpha }_{\mathrm{zz}})}^2+{({\alpha }_{\mathrm{zz}}-{\alpha }_{\mathrm{xx}})}^2\\ & +6{\alpha }_{\mathrm{xz}}^2+6{\alpha }_{\mathrm{xy}}^2+6{\alpha }_{\mathrm{yz}}^2{]}^{\frac{1}{2}}\end{array}$(1)

The static first hyperpolarizability (${\beta }_0$) which is defined as [82]: (2) $\begin{array}{r}\{\begin{array}{c}{\beta }_x={\beta }_{\mathrm{xxx}}+{\beta }_{\mathrm{xyy}}+{\beta }_{\mathrm{xzz}}\\ {\beta }_z={\beta }_{\mathrm{zzz}}+{\beta }_{\mathrm{xxz}}+{\beta }_{\mathrm{yyz}}\\ {\beta }_y={\beta }_{\mathrm{yyy}}+{\beta }_{\mathrm{xxy}}+{\beta }_{\mathrm{yzz}}\end{array}\{{\beta }_0={({\beta }_x^2+{\beta }_y^2+{\beta }_z^2)}^{1/2}\end{array}$(2) As noticed from Table 6 and Figure 7, and as reported in most literature [57,83], the values $\mathrm{\Delta }{\alpha }_o$are much smaller than those of ${\beta }_0$ for all molecules in both media. Since$\mathrm{\Delta }{\alpha }_o$ values are ranging from 86.02–108.60 au (E2 in gas), from 87.78–108.83 au (K1 in gas), from 114.93–146.91 au (E2 in water), and 118.07–147.56 au (K1 in water). While ${\beta }_0$ are ranging from 121.13–604.44 au (E2 in gas), from 189.15–475.58 au (K1 in gas), from 144.05–1812.64 au (E2 in water), and 586.56–1490.09 au (K1 in water). The aqueous medium resulted in a significant improvement in these values especially for E2 form and its derivatives [84]. Another important observation is the lower values of NLO parameters especially ${\beta }_0$ values for E2 and its derivatives with respect to those of K1 form and its derivatives. The enhanced NLO activity of a compound is always related to its lower energy gap. Since literature has highlighted an inverse relationship between energy gaps and hyperpolarizabilities [57,82–86]. Lower gaps and higher hyperpolarizabilities was correlated with higher reactivity and thus less stability. Therefore, the enhanced NLO parameters of K1 and its derivatives confirm their reactivity and less stability compared to E2 and its derivatives.

Figure 7. Anisotropic polarizability ($\mathrm{\Delta }{\alpha }_o$) and static first hyperpolarizability (${\beta }_0$) (in au unit) of (a) E2 and (b) K1 forms and their derivatives in the gas and water media calculated at CAM-B3LYP and PCM/CAM-B3LYP, respectively, with 6-311 + G(d,p) basis set.

Table 6. NLO parameters (in au) of the investigated FAV molecules calculated with CAM-B3LYP/6311 + G(d,p) in the gas phase.

	E2				K1
FAV	Δαo	^1ΔΔαo	βo	^2Δβo	Δαo	^1ΔΔαo	βo	^2Δβo
Gas Phase
FAV-H	86.02	12.51	197.08	1247.59	87.78	10.75	186.90	1257.77
FAV-F	86.40	12.13	121.13	1323.53	88.52	10.01	272.67	1172.00
FAV-Cl	99.63	−1.09	135.67	1309.00	100.85	−2.32	280.12	1164.55
FAV-Br	107.49	−8.96	172.99	1271.68	108.44	−9.91	189.15	1255.51
FAV-CN	104.23	−5.70	301.44	1143.22	104.95	−6.42	218.63	1226.04
FAV-NO2	104.48	−5.94	604.44	840.23	108.44	−9.91	189.15	1255.51
FAV-CHO	103.71	−5.18	589.61	855.05	105.08	−6.55	475.58	969.09
FAV-COOH	106.59	−8.06	485.17	959.50	108.14	−9.60	376.71	1067.96
FAV-CH3	99.36	−0.83	186.55	1258.12	100.93	−2.40	244.11	1200.56
FAV-NH2	98.47	0.06	240.54	1204.13	99.78	−1.25	444.77	999.90
FAV-OH	91.85	6.69	162.41	1282.26	93.50	5.03	358.04	1086.63
FAV-OCH3	105.20	−6.67	152.01	1292.66	106.52	−7.99	324.68	1119.99
FAV-SH	108.60	−10.07	207.32	1237.35	108.83	−10.30	192.44	1252.23
Water medium
FAV-H	114.93	−16.40	304.28	1140.39	118.07	−19.54	586.56	858.11
FAV-F	114.64	−16.11	144.05	1300.62	118.38	−19.85	800.44	644.23
FAV-Cl	132.98	−34.44	190.75	1253.91	135.37	−36.84	918.16	526.51
FAV-Br	144.13	−45.60	184.24	1260.42	146.14	−47.61	794.08	650.59
FAV-CN	139.25	−40.72	863.24	581.42	140.87	−42.33	681.61	763.06
FAV-NO2	141.51	−42.97	1812.64	−367.97	146.14	−47.61	794.08	650.59
FAV-CHO	139.19	−40.66	1396.07	48.60	142.17	−43.64	1202.05	242.62
FAV-COOH	141.74	−43.21	956.79	487.88	144.81	−46.27	812.20	632.47
FAV-CH3	132.51	−33.98	224.10	1220.57	135.62	−37.09	813.26	631.41
FAV-NH2	132.87	−34.34	566.38	878.28	136.03	−37.50	1490.09	−45.42
FAV-OH	122.64	−24.10	340.28	1104.39	126.31	−27.77	1193.91	250.75
FAV-OCH3	139.16	−40.63	357.56	1087.10	142.23	−43.70	1191.17	253.50
FAV-SH	146.91	−48.38	522.29	922.38	147.56	−49.03	664.98	779.69
	${}^1\mathrm{\Delta \Delta }{\alpha }_{\text{o}}=(\mathrm{\Delta }{\alpha }_{\text{o}}{)}_{\text{pNA}}-(\mathrm{\Delta }{\alpha }_{\text{o}}{)}_{\text{FAV}}$
	${}^2\mathrm{\Delta }{\beta }_o=({\beta }_o{)}_{\text{pNA}}-({\beta }_o{)}_{\text{FAV}}$

The NLO properties is mainly quantified by the values of static first hyperpolarizability $({\beta }_0)$, and as illustrated in Table 6 and Figure 7, among the E2 form, in the gas phase, the derivatives with CN, NO₂, CHO, COOH, NH₂, and SH groups recorded higher ${\beta }_0$ values compared to that of the parent form. While in water, additional derivatives included OH, OCH₃ groups recorded better NLO properties with higher ${\beta }_0$ values [87,88]. In contrast, all the derivatives of K1 form showed better ${\beta }_0$ values in both media. In addition, for the sake of assessment, the NLO parameters of the prototypical donor–acceptor molecule of p-nitroaniline (pNA) have been calculated at the same level of theory [89]. We find that the ${\beta }_0$ values of the two forms E2 and K1 and their derivatives are smaller than that of pNA apart from two derivatives E2-NO₂ and K1-NH₂ in water medium. The incorporation of NO₂ and NH₂ groups within E2 and K1 forms, respectively, is expected to enhance the charge separation and thus, hyperpolarizability due to the strong electron-withdrawing and electron-donating power of these two groups, respectively. These two derivatives are hence likely to be good candidates for NLO materials.

4. Conclusions

In this study, theoretical study of the relative stabilities of some FAV derivatives and their tautomers and rotamers have been systematically studied by using DFT calculation using B3lYP/6-311 + G(3df,2p)//B3LYP/6-311 + G(d,p) level. The results revealed that the enol form (E2) is predominant in the gas phase and aqueous solution. The relative stability of the other tautomers and rotamers with respect to E2 are increased in aqueous solution. The enol form of the parent FAV drug is more stable than the keto form by 6.77 kcal/mol. The geometrical properties, topological parameters, NBO analysis, and non-covalent interaction analysis have been used to investigate the stability of the enol form in comparison with the keto form K1, in terms of the strength of the IMHB. The IMHBs in E2 isomer are shorter than those in K1 form by an average of 0.229 Å and the bond angles of the IMHBs in E2 are wider than those in K2 form by an average of 14.1°, proving that IMHBs in E2 are stronger than those in K2 tautomer. The topological parameters indicate that the interaction between the proton donor and proton acceptor in E2 is partly covalent, while it is non-covalent in K2 tautomer. The average electron density at the bond critical point of the IMHB in E2 is 0.0051 au larger than that in K2 tautomer. Furthermore, the electron density at ring critical point of the pseudo ring in the case of E2 is larger than that in K1 tautomer by 0.0157 au. The IMHBs between the O and H atoms (E2 form) are stronger than those between the N and H atoms (K1 form) by an average of 6.73 kcal/mol. NBO analysis indicates that the charge transfer interactions between the lone pairs of electrons of the proton acceptor atoms and the proton donor antibonds are also larger in E2 formed compared to those in K2 tautomer. Non-covalent interactions clearly indicate that the interactions between the proton donor and proton acceptor in E2 are stronger than those in K1 tautomer. These results clearly show why E2 is energetically more stable than K1 tautomer. The inductive and resonance effects are important in the discussion of the effect of substitution on the geometrical, topological, and NBO parameters. Additionally, enhanced NLO parameters of K1 form and its derivatives confirm a lower stability of these compounds as compared to E2 form and its derivatives. In a water medium, E2-NO₂ and K1-NH₂ are likely to function well as NLO materials.

Authors’ contributions

Zaki Safi: Data curation, Formal analysis, Investigation, Methodology, Visualization, Software, Writing-original draft, Writing – review & editing, and Nuha Wazzan: Formal analysis, Investigation, Software, Visualization, Writing – original draft, Writing – review & editing. The two authors agreed to the final version of manuscript.

Acknowledgements

The authors gratefully acknowledge King Abdulaziz University’s High-Performance Computing Centre (Aziz Supercomputer) (http://hpc.kau.edu.sa) for assisting the calculations for the work of this paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The authors confirm that the data supporting this study's findings are available within the article and its supplementary data. Additionally, the datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Additional information

Funding

This project was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah [grant number KEP-PhD: 93-247-1443]; the authors, therefore, gratefully acknowledge DSR technical and financial support.