{'fraction_non_alphanumeric': 0.030385819950689678, 'fraction_numerical': 0.015326181115479443, 'mean_word_length': 4.724775891665077, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Inclusion of Levodopa into β-Cyclodextrin: A Comprehensive Computational Study

This study focused on the inclusion of levodopa (LVDP) into β-cyclodextrin (BCD) using various computational methods such as quantum mechanics (QM), molecular dynamics/steered molecular dynamics (MD/SMD), and QM/molecular mechanics/Poison–Boltzmann surface area (QM/MM/PBSA). The QM results assigned the most significant charge-transfer atoms and the higher stability of LVDP in the aqueous phase. The MD results indicate the formation of a 1:1 complex with a reasonable estimation of the effective radius of the complex, the significant contribution of hydrogen bonding in the binding energy, and the enhancement of the water solubility of LVDP. By accounting for the water hydrogen bonds and their dipolar effects, QM/MM calculations lead to the more accurate IR spectrum and binding energy of the BCD–LVDP complex. By considering carboxylic and amine functional groups’ more precise arrangement, QM/MM assigns stronger hydrogen bonds between LVDP and BCD. While all the methods provide a reasonable estimation of the binding energy, the most accurate value (−4.14 kcal/mol) is obtained from QM/MM/PBSA.


INTRODUCTION
Today, a chemical drug that is a type of precursor and is known as Levodopa (LVDP) and L-Dopa (Figure 1a) is used and prescribed as a medicine for Parkinson's disease. Parkinson's disease is one of the most common diseases in different parts of the world that mostly appear in old age. 1,2 LVDP was able to cause immediate improvement in patients with Parkinson disease due to its conversion to dopamine in both the central nervous system and the peripheral nervous system. 3−5 However, excessive amounts of this drug can cause side effects such as tremors and dyskinesia. The poor aqueous solubility (3.3 mg/mL) and the degree of dissolution of LVDP are two critical factors that affect its formulation and development process and limit its therapeutic index. Therefore, we need to use a suitable drug delivery system to overcome these shortcomings. 6−8 Drug delivery is a process for the proper transfer of medicinal compounds into the target, to increase therapeutic efficacy, to reduce toxicity, to reduce side effects, and to increase solubility. 9,10 Cyclodextrins (CDs) as highly biocompatible molecules with a hydrophilic outer surface and lipophilic cavity are one of the most popular drug careers that are widely used in the pharmaceutical industry. 11,12 They have the shape of a truncated cone because of the chair conformation of glucopyranose units among natural cyclodextrins (α-, β-, and γ-CDs, consisting of 6, 7, and 8 glucopyranose units, respectively). 13 β-Cyclodextrin (BCD) (Figure 1b) is the best drug career candidate for a variety of drugs due to the perfect cavity size, efficient drug complexation and loading, availability, and relatively low cost preparation. 14 The complexation of LVDP with BCD has been experimentally investigated by several research groups, indicating the formation of a 1:1 complex with favorable free binding energy (ΔG b ) that enhances the solubility, bioavailability, and stability of LVDP. 15,16 These results indicate that BCDs are promising drug delivery systems for this study.
Computational methods such as quantum methods, molecular dynamics (MD), and hybrid quantum mechanics/ molecular mechanics (QM/MM) provide reasonable predictions for the researcher in examining the chemical systems and reveal some of the molecular details that could not be assigned in the experimental methods. 17,18 There are many studies that evaluated the performance of these methods. 15,19−25 Among these methods, QM methods based on solving the Schrodinger equation are the most accurate. However, due to the high computational cost, they can be applied for systems with a small number of atoms. 26 These methods are essential in studying the properties of the electronic structure. MD methods are based on solving Newton's second law equation. They have a relatively high speed and a low computational cost that make them a good fit for evaluating systems with many atoms. This method provides good predictions about the system's mass characteristics but could not be used to study processes involving electronic transfers. However, QM/MM, as a hybrid method, has the most advantages of QM and MD methods but still has a high computational cost. In drug-BCD studies, as many-atom systems, QM methods are rarely used, and most MD methods and sometimes hybrid methods are usually applied. 19,23 Although there is no comprehensive computational study on BCD−LVDP, we have performed such a study in this endeavor. For this purpose, several computational methods such as QM, MD/steered molecular dynamics (SMD), and QM/MM/Poisson−Boltzmann surface area (PBSA) have been used in parallel and comparatively. The complete results of this study reveal the molecular and energetic details of the BCD− LVDP complex and show how BCD could improve the bioavailability and solubility of LVDP. Moreover, according to the experimental results, the accuracy of the applied methods has been assessed and compared confidently.

COMPUTATIONAL METHODS
2.1. Quantum Mechanics. Initial geometries of LVDP (CID: 6047) and BCD (CID: 444041) are obtained from the PubChem site. 27 All structures are fully optimized using the modified density functional theory with the CAM-B3LYP function. 28 CAM-B3LYP is a hybrid exchange−correlation functional that combines B3LYP at short range with an increasing level of exact Hartree−Fock exchange at long range. 29 The double zeta split valence polarization with diffuse basis set functions (aug-ccpVDZ) was applied for all atoms. 30 To perform calculations in the aqueous phase, the conductorlike polarizable continuum model (CPCM) is used as a common water model. This is a kind of PCM that considers the charge-transfer effect as a part of solvation. 31 Some parameters such as structural properties, charge transfer, binding energy, electrostatic potential, normal vibrational frequency, molecular orbitals, and natural bond orbital (NBO) analysis related to LVDP, BCD, and the complex were investigated by the QM method. For quantum calculations, Gaussian09 32 software has been used, and for visualization, Avogadro, 33 Gausview06, 34 and Chimera 35 software have been used. At first, according to the drug and NBO data's optimized structure, the electronic properties of LVDP, such as partial charge distribution, bonding lengths and angles, dipole momentum, and energies, are estimated. Then, by considering the optimized structure of BCD in the vicinity of LVDP, the binding energy is calculated. The scan calculations were performed to confirm the optimization step. The interacted points in LVDP and BCD are explored using the electrostatic potential map. Frequency calculations for LVDP identification in the BCD−LVDP spectrum were applied in both phases in accordance with the optimization method. 18,19 Finally, the nature of the interaction is evaluated through NBO calculations.
2.2. Molecular Dynamics. All MD simulations are carried out with the GROMACS-2016 package 36 using AM-BER99SB 37 as a force field and SPC/E as a flexible water model. 38 The GROMACS package analysis tools are used to interpret the dynamics data where all simulation snapshots are visualized by VMD software. 39 Initial parameters of LVDP and BCD for the AMBER99SB force field were determined using the TLEAP utility implanted in AmberTools18. 40 In general, we used a cubic box and periodic boundary conditions in all simulations. The box lengths are chosen from 5 to 7 nm corresponding to the number of species. All species are randomly placed inside the simulation box. The following steps were considered for all simulations: (1) energy minimization (EM), (2) NVT equilibration, (3) NPT equilibration, and (4) MD production. For EM, we used the steepest descent algorithm with 50,000 steps. A middle-temperature coupling using a velocity-rescale 41 thermostat with a coupling constant of 0.1 ps and a middle-pressure coupling using a Berendsen barostat 42 with a coupling constant of 1.0 ps were applied for NVT, NPT, and MD runs. The reference temperature and pressure are 298 K and 1 atm, respectively. A cutoff radius of 1.4 nm is adopted for nonbonded interactions, short-range electrostatics, and short-range neighbor lists. However, the long-range electrostatics is handled by the particle mesh Ewald (PME) method. 43 Vibrations of bonds containing hydrogen atoms are constrained using a LINCS algorithm. 44 The leapfrog algorithm is employed for integrating Newton's equations of motion with a 2.0 fs time step. The neighbor lists are updated automatically (due to CUDA processing) every 40 times in the equilibration and production runs. The NVT and NPT equilibrium stages consist of 1.0 and 10.0 ns run times, respectively. The MD production stage consists of a 100 ns run time. The trajectories are stored every 5 ps. The quantities such as radius of gyration (R g ), solvent accessible surface area (SASA), radial distribution function (RDF) g(r), hydrogen bond, and principal axes of inertia are calculated based on the output data (trajectory). The binding energy in the previous step is obtained as the energy difference between the product and the raw materials. However, in this step, the bonding energy is derived from the potential of mean force (PMF), which is obtained from a series of sequence simulations, which is known as Umbrella Sampling or SMD. 45,46 In this method, after the system reaches equilibrium, we allow the center of mass (COM) of the molecule (LVDP) to move through the center of BCD to obtain a series of distances related to the potential energy curve. Next, we extract about 400 frames from the obtained trajectory. Considering the calculated BCD− LVDP distance for each frame, only 32 frames that show appropriate variation in the distance are selected. Then, for each of the 32 selected frames, we run a 10 ns simulation. Finally, we use the weighted histogram analysis method (WHAM) 47 to extract the PMF and ΔG b .
2.3. Hybrid Method (QM/MM). After MD simulation for 100 ns, the QM/MM computations were carried out to ensure that the system was in equilibrium. For QM/MM computation, Gromacs2016 and ORCA software were compiled together using a modified script. 48 A 10 ns QM/MM simulation was performed, which obtained 200 snapshots from the last 5 ns to calculate the binding energy and investigate the system dynamics.
In the QM/MM method, LVDP is described quantum mechanically (B3LYP/6-31G*), while its environment, which consists of BCD and the water solvent, is described molecular mechanically (Amber99SB force field). A QM calculation is performed for LVDP using electronic embedding that describes the electrostatic interactions between the electrons of the QM region and the MM atoms and between the QM nuclei and the MM atoms included in the Hamiltonian for the QM subsystem. 49 Nonelectrostatic interactions between the QM and MM subsystems are calculated through a Lennard-Jones potential using parameters from the classical force field. 50,51 Single-point interaction energy calculations based on the time averaged and changes in the internal properties of the LVDP structure were examined for each of the snapshots obtained from the QM/MM MD calculations. The binding energy of BCD−LVDP was calculated using a QM/MM/ PBSA method that was externally implemented in GRO-MACS. 52 As the entropy contribution is not considered in this method, the normal mode frequency, which requires the calculation of eigenvectors in the Hessian matrix, is used to estimate the entropy. As usual, the binding energy was defined as follows: 53 where E complex , E solvation , and TΔS are the interaction energy, solvation energy, and entropy contribution for forming BCD− LVDP from its separated constituents, respectively, and using the QM/MM level for estimation of G binding leads to a more accurate value.

RESULTS AND DISCUSSION
In this study, we used CPCM as an implicit water model to evaluate the behavior of drugs and drug carriers in an aqueous environment. The LVDP molecule is optimized in both gas and aqueous media (pH ≈ 7), and the schematic representation of the optimized structures and the relevant electronic structural properties are shown in Figure 2 and Tables 1, and 2, respectively. It should be mentioned that the predicted structure and ionization state of LVDP in the aqueous phase by ChemAxon software has been used as the initial structure for final optimization by Gaussian. According to Table 1, the most significant change in the estimated partial atomic charge from NBO calculations, as an index of charge transfer from the gas phase to the aqueous medium, refers to the oxygen atom of the carboxyl group (O2) and the nitrogen atom of the amino group (N5). However, the lowest charge transfer belongs to alpha carbon in the amino acid group (C7). C12−O4 shows the most significant change in the carboxylic acid group regarding bond lengths, which represents the increase of bond order in the aqueous media. The lower value of LVDP's chemical potential, as an index of electronic structure stability, in the aqueous phase (−3.03 kcal/mol) in comparison to the gas phase (−2.88 kcal/mol), represents the higher stability of LVDP in the aqueous phase. The higher value of the energy gap in the aqueous phase (4.54 kcal/mol) concerning the gas phase (4.48 kcal/mol) also represents the higher stability in the aqueous phase, in agreement with the chemical potential results.
In the next step, the interaction of LVDP with BCD is investigated in both gas and aqueous phases using the optimized LVDP structure. There is a significant difference between the binding energies in the gas (0.03 kcal/mol) and aqueous (−1.57 kcal/mol) phases, which could be attributed to electrostatic interactions and hydrogen bonding in the aqueous medium, which are discussed below. Usually, for a more detailed study of the optimized structure and energy, the two components of the complex are placed at different distances, and the complex's energy changes are evaluated. As shown in Figure 3, the potential energy curve obtained from the displacement of the LVDP molecule in the Z-axis has the lowest amount of energy at the optimal point ( Figure 3b). As the phenolic groups of the drug approach the narrow head of BCD, the energy level increases sharply due to the increase of the electrostatic repulsion ( Figure 3a). Furthermore, as the phenolic groups move away from the wide head of BCD, the strength of the interactions decreases, resulting in zero energy ( Figure 3c). On the other hand, the drug molecule rotation in BCD does not cause more change in energy.
The molecular electrostatic potential method (MEP) is an efficient computational tool used to estimate the potential energy surface as an essential indicator in molecular interactions and chemical reactions. 54,55 Figure S1 is the molecular electrostatic potential map of BCD and LVDP molecules. As shown in Figure S1a, the inner surface of BCD has hydrophobic properties due to its low electrostatic potential (neutral charge distribution) and its outer surface due to its high electrostatic potential (positive and negative charge distribution) has hydrophilic properties. Also, in the LVDP molecule, the amide group has hydrophilicity due to its high electrostatic potential, and the phenolic group shows hydrophobicity due to benzene ring groups. Therefore, the drug molecule will approach the inner surface of BCD from the phenolic group.
Harmonic vibrational frequencies corresponding to the optimized BCD structure were calculated in the presence and absence of LVDP in the gas phase. Figure 4, as IR spectra of BCD and BCD−LVDP, shows the results. The red, blue, and green spectra are related to BCD, BCD−LVDP, and LVDP in the aqueous phase, respectively. As can be seen in the blue spectrum, the specific peaks in points a, b, c, d, e, and f are   Figure 4]. The other specified peaks are also listed in Table 4. The difference between the observed values and the expected values is higher for the more polar bonds that could be related to not considering intermolecular interactions in the gas phase. 57 The transferred charge quantity from a donor−acceptor orbital and its stability are estimated using the second-order perturbation theory, which is implemented in NBO calculations. The results are listed in Table S1. The low values of intermolecular charge transfer (Nos. 5, 6, and 7) represent the weak hydrogen bonds between BCD and LVDP and the predominant role of electrostatic interactions in the binding energy.
In the next step, the interaction of LVDP with BCD is investigated in MD simulation. Two main parameters are used to evaluate the system equilibrium conditions, which are rootmean-square deviation (RMSD) and radius of gyration (R g ). rmsd is a measure of the amount of displacement between atoms or different groups over time and is usually used as the main parameter in the equilibrium of molecular dynamic systems. R g is also used as a measure of collapse in some  structures such as proteins, polymers, and micelle formation. 58,59 The 100 ns time regime of the calculated rmsd and R g are shown in Figure 5a. The fluctuations have been substantially decreased for both parameters after the 8 ns simulation that guarantees the achievement of equilibrium after 100 ns. The mean values of rmsd and R g in the last 20 ns simulation are 0.52 and 0.56 nm, respectively. The effective radius (R s ) that is related to R g by R s ≈ 1.3R g 60,61 equals 0.73 nm (7.3 Å).
Another parameter attributed to solubility is the SASA; 62,63 as the SASA goes up, this means that a compound solubility increases, or if something is complexed by BCD, we would expect the higher SASA. The time variation of SASA for LVDP, BCD, and BCD−LVDP is given in Figure 5b. The least average value of 1.9 nm 2 belongs to LVDP, indicating the least aqueous solubility (≈3.3 mg/mL), while the highest average value of 11.9 nm 2 is for BCD, representing the highest aqueous solubility (≈18.2 mg/mL). The higher value of SASA for BCD−LVDP (11.1 nm 2 ) with respect to LVDP, representing the solubility of LVDP, goes up due to the complexation with BCD. The lower value of SASA for BCD−LVDP (11.1 nm 2 ) concerning the free BCD (11.9 nm 2 ) could be due to the water replacement by LVDP in the inner cavity of BCD. The water displacement will be further investigated by calculating the normal RDF of water in the subsequent section. However, these results show up the appropriate BCD for the complexation of LVDO and enhancement of its aqueous solubility.
The normal distribution of LVDP around BCD and water molecules around BCD and LVDP have been examined in the next step. For this purpose, the normal RDF module, which considers the COM of the molecules, has been used. The variation of the mean RDF values in the last 10 ns versus distance from COM of BCD is plotted in Figure 6. According to this figure, the maximum distribution of LVDP is observed at 0.57 nm from COM−BCD (blue curve). The distribution of water molecules around BCD in the absence (green graph) and presence (purple graph) of LVDP represents a significant amount of water in the hydrophobic cavity of BCD in the absence of LVDP. However, most of these water molecules have been displaced by LVDP due to the insertion of LVDP into the hydrophobic cavity of BCD and the formation of the complex.
The comparison of green and purple diagrams in Figure 6 represents the increase of water distribution in the presence of LVDP in the cavity of BCD. This observation could be related to the amino group of LVDP that interacts with water molecules via hydrogen bonds. The R s of BCD−LVDP can be estimated from the intersection point of blue and green diagrams. 64 The intersection point appears at 0.73 nm, which agrees with R s ′ obtained value from the previous section. According to the water distribution diagram around BCD (green), the radius of the outer water layer (R o ) is about 1.49 nm. The comparison of R s and R o values of BCD in the presence and absence of LVDP (Table 3)  Another critical analysis between LVDP and its environment is the hydrogen bond. By setting 0.35 nm and 10°for the cutoff radius and cutoff angle, respectively, the values of 3 and   65 7 are obtained for the average number of hydrogen bonds between LVDP−BCD and LVDP−water molecules, respectively. The two oxygens of the carboxylic groups and the oxygen of the phenolic group in LVDP, as hydrogen acceptors, and three OH groups of BCD, as hydrogen donors, formed these hydrogen bonds between LVDP and BCD. The locations of these hydrogen bonds are shown in Figure 7. Figure  As mentioned before, one of the most common analyses, conducted by umbrella sampling simulations, is the PMF extraction, which yields the ΔG for the binding and unbinding processes. The optimized structure of the BCD−LVDP complex, from the previous stage by 100 ns MD production, was considered the starting point. LVDP is moved from its equilibrium position in the Z-axis direction at a speed of 0.01 nm/ps during a 200 ps time evolution, and 400 frames were extracted. As the X-and Y-axis displacements do not significantly change the potential, the movement in these directions has not been considered. The distance between COM−LVDP and COM−BCD was measured for these 400 frames. Subsequently, the 38 frames were selected by a 0.05 nm changing distance for the umbrella sampling stage. Then, each frame was simulated using 1 ns in the NPT ensemble for equilibration and 10 ns in MD production. The calculated PME values using WHAM analysis are shown in Figure 9. The difference between the lowest and highest values of PME is about −3.67 ± 0.05 kcal/mol, which is a reasonable estimation for the BCD−LVDP binding energy obtained. 44 In the last part of this work, we focus on the results of the QM/MM method for the geometrical structure, dynamics properties, IR spectra, and energy parameters of the LVDP− BCD complex. The time evolution of distance between COM of LVDP and BCD, as a structural parameter, was calculated by MD and QM/MM methods, and the results are shown in Figure 10. The average distances between the COM of LVDP and BCD in the last 20 ns of MD and the last 5 ns of QM/MM are about 0.12 and 0.09 nm, respectively. By considering the experimental value, 0.10 nm, obtained from the X-ray structure, 19 QM/MM leads to a more accurate result.
The time evolution of θ and φ dihedral angles, as other structure properties, has been evaluated for 200 snapshots of the last 5 ns of QM/MM calculations, and the results are shown in Figure 11, representing the higher degree of freedom. In order to estimate the most probable range of θ and φ, their normal distribution functions are evaluated and shown in Figure 12. The most probable value for θ is observed at ±24°, while that for φ is observed at 28°. The θ value dispersion is much higher, representing its higher degree of freedom.
As mentioned before, the study and assay of spectroscopic spectra such as IR is one of the most compelling cases in computational studies and covers a wide range of studies. In this study, we intend to investigate the effects of the solvent  and mainly hydrogen bonding on LVDP in the presence of BCD. Figure 4 shows the characteristic peaks of the gas-phase LVDP in BCD, labeled with a to f, and for free LVDP in the aqueous phase, labeled with a' and d'. The IR spectrum of the LVDP−BCD complex has been calculated based on QM/MM, and the integrated nonscaled frequencies are shown in Table 4. The differences between the obtained values from the gasphase calculations and QM/MM could be related to the solvation and hydrogen bonds. For instance, the estimated frequency of CO bound in LVDP from QM/MM is 1670 cm −1 , while the value in the gas phase is 1846 cm −1 , which is due to the formation of a hydrogen bond with the water molecule and interactions with BCD (peak d). The difference between d (1670 cm −1 ) and d′ (1658 cm −1 ) values represents the less hydration of CO in the BCD cavity. We could provide the same interpretation for the phenolic O−H group of LVDP from QM/MM (3552 cm −1 ) and gas-phase calculations (3801 cm −1 ). However, the relatively high difference between a (3552 cm −1 ) and a′ (3795 cm −1 ) represents the weakness of the QM method in accounting for the solvation effect.
The estimation of the binding energy is the last part of this endeavor. There are two approaches for this purpose; in the first, the average energy of free BCD (⟨E QM/   7 kcal/mol, respectively. By considering these quantities and using eq 2, the value of −3.75 ± 0.01 kcal/mol is obtained. As explained in the section on Computational Methods, the main shortcoming of this method is ignoring entropy contribution.
In the second approach, the binding energy has been estimated using the QM/MM/PBSA method, where the entropy contribution is brought up. The estimated binding energy from this approach is −4.14 kcal/mol. According to eq 1, the estimated values of the energy components from the QM/MM/PBSA method are listed in Table 5.
A comparison among the estimated binding energy values obtained from the considered methods and their deviation from the experimental value is provided in Table 6. According to the results, the following order is extracted for the accuracy of the methods: QM Gas phase < QM Implicit water < SMD ≈ ≈QM/MM ≪ QM/ MM/PBSA. While the highest accuracy refers to QM/MM/  PBSA, it has the highest computational cost. The more accuracy of QM/MM/PBSA concerning QM/MM represents that the addition of PBSA analysis leads to more realistic results that could be due to the importance of solvation. Interestingly, the difference between SMD and QM/MM results is not significant, while the computational cost of QM/ MM is much higher.
Different computational methods can complement each other well. In other words, each of these methods can be useful in such a way. For example, in this study, a quantum method was performed to extract drug and carrier parameters (such as length, angle, and partial charge) and the initial guess of the interaction orientation using an MEP. According to Figure S1a, it can be predicted that the inner surface of BCD has hydrophobic properties, and the phenolic group of LVDP shows hydrophobicity. Therefore, it is more probable that the drug molecule will approach the inner surface of BCD from the phenolic group, which is consistent with the results obtained from the MD simulation. The NBO analysis shows three significant charge transfers between BCD and LVDP, representing the minor role of hydrogen bonds and the predominant role of electrostatic interactions in the binding energy. At the same time, the MD calculations lead to the result of more contribution of hydrogen bonds. The obtained structural properties (the distance between COM and dihedral angles) from MD and the hybrid method are very close to the experimental values, representing the close accuracy of these methods for estimating the structural parameters. However, the QM method does not provide a reasonable estimate. While SMD does not provide the MD features, it is considered a complementary method (MD/SMD) in calculating the binding energy.

CONCLUSIONS
This work focuses on the molecular aspects of BCD−LVDP interactions using various computational techniques with different advantages and disadvantages. This strategy provided a comprehensive perspective of this system. The comparison of properties such as chemical potential, energy gap, and binding energy in the QM method for both gas and aqueous phases represents the higher stability of LVDP in the aqueous phase due to electrostatic interactions. The QM method also predicts the presence of only two weak hydrogen bonds between LVDP and BCD in the aqueous phase, which belong to the LVDP phenolic groups, and does not provide any information about the other possible hydrogen bonds that may be formed by amine and carboxylic groups. Moreover, this method does not assign the apparent arrangement of the amino and carboxylic groups of LVDP in the BCD cavity, could not estimate the accurate IR frequency bands, and could not provide a realistic picture of the solvent effects and hydrogen bonds. The MD calculations with the explicit water model ran in the next part of this endeavor to overcome these deficiencies. The results of MD calculations indicate the formation of a 1:1 complex with a suitable internal and external radius that certifies the ability of the complex to cross the blood−brain barrier. The location of LVDP at the center of BCD, the significant role of hydrogen bonds in the complex formation, the highest distribution of water molecules around the amino and phenolic groups of LVDP, and the enhancement of water solubility of LVDP in the presence of the BCD are the other results of the MD calculations.
The RDF analysis represents the water molecule's removal from the cavity and the redistribution around BCD due to the insertion of LVDO. This observation led us to account for the effects of hydrophilicity and hydrophobicity on the binding energy using QM/MM/PBSA as a more precise and accurate method. The value of −4.14 kcal/mol is estimated for the binding energy from QM/MM/PBSA calculations. Moreover, the results represent the higher contribution of the hydrophobic effect concerning the hydrophilic effect in the binding energy.   As the QM/MM method reasonably considers the dipolar effects of water and hydrogen bonding, a more accurate IR spectrum of the complex is estimated from QM/MM calculations. Moreover, by considering the precise arrangement of carboxylic and amine functional groups, QM/MM assigns a higher number of strong hydrogen bonds between LVDP and BCD. However, the MD method's estimated values are not reliable due to the high fluctuations.