# Fileset

[molecules30(2025)2732.pdf](https://mdr.nims.go.jp/filesets/a325e22b-b9b6-4427-950e-a84d854e8550/download)

## Creator

Natarajan Sathiyamoorthy Venkataramanan, Ambigapathy Suvitha, [Ryoji Sahara](https://orcid.org/0000-0003-0788-2985)

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Intramolecular Versus Intermolecular Bonding in Drug Gemcitabine and Nucleobases: A Computational Study](https://mdr.nims.go.jp/datasets/34ac21b5-f0ff-4f34-bca1-54f5febdbc74)

## Fulltext

Intramolecular Versus Intermolecular Bonding in Drug Gemcitabine and Nucleobases: A Computational StudyAcademic Editor: Spyridon MourtasReceived: 30 May 2025Revised: 16 June 2025Accepted: 21 June 2025Published: 25 June 2025Citation: Venkataramanan, N.S.;Suvitha, A.; Sahara, R. IntramolecularVersus Intermolecular Bonding inDrug Gemcitabine and Nucleobases:A Computational Study. Molecules2025, 30, 2732. https://doi.org/10.3390/molecules30132732Copyright: © 2025 by the authors.Licensee MDPI, Basel, Switzerland.This article is an open access articledistributed under the terms andconditions of the Creative CommonsAttribution (CC BY) license(https://creativecommons.org/licenses/by/4.0/).ArticleIntramolecular Versus Intermolecular Bonding in DrugGemcitabine and Nucleobases: A Computational StudyNatarajan Sathiyamoorthy Venkataramanan 1,* , Ambigapathy Suvitha 2 and Ryoji Sahara 3,*1 Department of Chemistry, School of Engineering, Dayanada Sagar University, Bangalore 562 112, India2 Department of Chemistry, SKP Engineering College, Tiruvannamalai 606 601, India; suvithaa@gmail.com3 Computational Structural Materials Group, Research Center for Structural Materials, National Institute forMaterials Science (NIMS), 1-2-1, Sengen, Tsukuba 305-0047, Japan* Correspondence: venkataramanan-chem@dsu.edu.in or nsvenkataramanan@gmail.com (N.S.V.);sahara.ryoji@nims.go.jp (R.S.)AbstractThe adsorption of the drug gemcitabine on nucleobases was investigated using a dispersion-corrected density functional theory (DFT) study. The planar structure of complexes is morestable than those with stacked and buckle-angled configurations. The complexes werefound to possess at least two intermolecular hydrogen bonds. The binding energy andinteraction energy are both negative, with the highest values observed for the gemcitabine–guanine and the lowest in the gemcitabine–thymine complex. The complex formation wasfound to be an enthalpy-driven process. Pyrimidine nucleobases have a lower enthalpyof formation than purine nucleobases. The computed HOMA and NICS values on thegemcitabine–nucleobase complexes show a substantial increase compared to the pristinenucleobases. An MESP analysis of the complexes shows a directional interaction andelectron density shift between the gemcitabine and the nucleobases. A QTAIM analysisindicates that the intermolecular hydrogen bonds have a partial covalent character. Thecomputed bond energy demonstrates that intermolecular NH···N bonds are more potentthan other bonds. An energy decomposition analysis using the DLPNO−CCSD(T) methodindicates that the complexes exhibit a substantial electrostatic attraction, and dispersioncontributes the least towards the system stability. The intermolecular bonds are strongerthan the intramolecular bonds in the drug–nucleobase complexes. The strength of in-tramolecular bonds is determined by the deformation of the gemcitabine ring during thecomplex formation.Keywords: drug; nucleobases; non-covalent interactions; DFT calculations; energydecomposition analysis1. IntroductionCancer is the second leading cause of mortality globally, and one of the primarydiseases affecting millions of people. Unlike normal cells, cancer cells continue to growand divide without regulation. Cancer cells spread from the growth site to other parts ofthe body through the blood or lymphatic system, forming new cancer cells. According tothe World Cancer Research Fund International, 510,992 new pancreatic cancer cases werereported in the year 2022. The most commonly used treatments involve chemotherapy,surgery, and radiotherapy, or their combinations [1–3]. Chemotherapy remains a funda-mental tool because of its simplicity. Among the chemotherapy drugs used, gemcitabinehas been an active agent against colon, pancreatic, ovarian, breast, head, neck, and lungMolecules 2025, 30, 2732 https://doi.org/10.3390/molecules30132732https://doi.org/10.3390/molecules30132732https://doi.org/10.3390/molecules30132732https://creativecommons.org/licenses/by/4.0/https://creativecommons.org/licenses/by/4.0/https://www.mdpi.com/journal/moleculeshttps://www.mdpi.comhttps://orcid.org/0000-0001-6228-0894https://orcid.org/0000-0002-6952-9784https://orcid.org/0000-0003-0788-2985https://doi.org/10.3390/molecules30132732https://www.mdpi.com/article/10.3390/molecules30132732?type=check_update&version=1Molecules 2025, 30, 2732 2 of 21cancers [4,5]. Gemcitabine (2′2′-difluorodeoxycytidine, dFdC) is a novel nucleoside analogof deoxycytidine. Originally, gemcitabine (GEM) was produced as an antiviral agent beforebeing established as an anticancer medication based on the remarkable in vitro and in vivoantitumoral activity [6].Gemcitabine is regarded as a gold standard and the first FDA-approved drug as amonotherapy for pancreatic cancer, which has a long-term survival rate close to zero. Themode of action of gemcitabine arises in the following two pathways: (i) by incorpora-tion into the DNA/RNA inhibiting excision enzyme resulting in replication retardation;(ii) blocks DNA synthesis through inhibition of ribonucleotide reductase enzyme ham-pering the synthesis of new strands [7]. Gemcitabine is identified as a prodrug, withits activation requiring a series of phosphorylation events. The drug is phosphorylatedmainly to gemcitabine monophosphate (dFdCMP) by deoxycytidine kinase, then turningto gemcitabine diphosphate (dFdCDP) and metabolite gemcitabine triphosphate (dFdCTP)by catalysis of other pyrimidine kinases before acting on DNA/RNA [8–10]. The half-lifeof gemcitabine in the human body ranges from 32 to 94 min, during which time it is deami-nated to 2′,2′-difluoro-2′-deoxyuridine (dFdU). Cytidine deaminases potentially inactivategemcitabine by deamination, which can contribute to resistance against the drug in cancertherapy [11,12]. Hence, understanding the mode of interaction of gemcitabine with thenucleobases at the molecular level is essential for improving chemotherapeutic efficacy.Chemotherapy medications are often classified as either groove binders or DNAintercalators, which can alter the duplex structure’s properties and have therapeutic effects.Historically, computational studies have concentrated on predicting the interaction betweenDNA nucleobases and pharmaceuticals [13]. Hence, several theoretical studies have beenconducted to know the mode and nature of the interaction between 5-fluorouracil andnucleobases [14–17]. We recently found that 5-fluorouracil bind with Watson–Crick basepairs more firmly than Hoogsteen base pairs due to higher hydrogen bond stability in theformer complex [18]. Mazumdar et al. examined the base pairing of two adenine analognucleobases utilizing various orientations, considering that numerous antiviral medicinesare recognized nucleobase analogs [19]. Among the orientations studied, the Watson–Crickbase pair possess NH-O/N interaction and contributes more to the interaction energy.Density functional studies on gemcitabine are limited, and most of them focus on theadsorption and delivery of gemcitabine using nanosystems as drug carriers [20–24]. Tounderstand the mechanisms of the formation of C3′ and C2′ sugar radicals, the one-electronoxidation of gemcitabine has been studied by Adhikary et al. [25]. Theoretical studieshave shown that gemcitabine to C3′ conversion needs a substantial barrier, which are inagreement with the experimental studies. Very recently, Malekzadeh et al., using the DFTmethod B3LYP/6-31+G*, studied the electronic structure of DNA nucleosides, 2′-diflouroroanalogues, and gemcitabine [26]. Structural parameters show only a slight variation inthe structure, while sugar puckering was noticed in gemcitabine. Ali et al. reported theexperimental and computational investigation of gemcitabine with bovine serum albumin(BSA) [27]. Gemcitabine was found to form a 1:1 stable complex with BSA by the MDsimulation study. The secondary structure of BSA remained intact while the drug wasfound to bind at the primary binding site 2.In the present work, we have carried out a DFT study to understand the bindingsite of gemcitabine with purine and pyrimidine nucleobases, and the nature of hydrogenbonding that exists between them. To identify the ground state structures, we conducted athorough search by generating initial structures using the artificial bee colony algorithm.The obtained structures were fully optimized using the M06-2X-D3/6-311+G(d,p) level oftheory. Their stability, along with the thermodynamic parameters of their formation, werealso computed. The nature of the interaction between the drug and the nucleobases wasMolecules 2025, 30, 2732 3 of 21studied using quantitative molecular electrostatic potential (MESP), aromaticity descriptors,quantum theory of atoms in a molecule (QTAIM), and non-covalent interaction analysis.This work provides a new insight into the intercalation mode of the drug gemcitabine withthe nucleobases, and it will be helpful to reveal the mechanism of the drug action, and canfurther the design of new anticancer drugs.2. Computational MethodsThe Gaussian G16 Rev. C.01 software package was used for all density functionaltheory computations [28]. To create initial geometries, we used the artificial bee colonytechnique for cluster global optimization (ABCluster) [29,30]. The M06-2X/6-311+G(d,p)approach with Grimme’s D3 dispersion correction was used for the initial optimization [31].For the thermochemistry of main group elements, the Minnesota functional M06-2X hasbeen tailored to consider dispersion interactions, and performs exceedingly well [32].Thanthiriwatte et al. showed that meta-GGA functionals M06-2x and functionals withdispersion corrections provide very good potential energy curves for the prototype bimolec-ular complexes of formic acid, formamide, and foramidine [33]. Furthermore, intermediateand long-range correlation and dispersion must be accounted for when handling hydrogen-bonded systems. To implement implicit solvation effects, the polarizable continuum model(PCM) was used. To ensure that the optimized geometries are at a minimum on the poten-tial surface and to meet the thermodynamic criteria, harmonic vibrational frequencies arecalculated for them. At 298 K and 1 atm of pressure, the thermodynamic parameters forthe most stable model were computed for both the gas and solution phases. The electronicenergy difference between the optimized gemcitabine–nucleobases complex and the sumof the optimized geometries of the gemcitabine molecule and the nucleobases was used tocompute the binding energy (BE) of the gemcitabine molecule with the nucleobases.The wave functions were developed at the M06-2X-D3/6-311+G(d,p) level of theory fortopological analysis using the AIMALL program in order to obtain a fuller understandingof the nature of the bonding [34]. The electron density (ρ) and its Laplacian (▽2ρ), the totalelectron energy density (H), the potential electron energy density (V), and the Lagrangiankinetic energy (G) were used to determine the properties of the bond critical point (BCP).Our earlier articles [35] provide more information on the AIM analysis. The Wave FunctionAnalysis–Surface Analysis Suite (WFA-SAS) was utilized to compute and visualize the three-dimensional surface of the quantitative molecular electrostatic potential (MESP) [36–38].Any region’s sign is dictated by the impact of the local electrons (negative) or nuclei(positive). The local greatest positive and most negative values (of which there may beseveral) of V(r) are denoted as VS,max and VS,min, respectively, when plotted on a molecularsurface. Using the gauge-independent atomic orbital (GIAO) approach, the 1H and 13CNMR chemical shift, as well as the nuclear independent chemical shift (NICS) valuesof stable complexes were calculated in water using the M06-2X/6-311+G(d,p) level oftheory [39]. The Multiwfn program was used to perform non-covalent interaction (NCI)and reduced density gradient (RDG) analyses, and the Chemcraft program was used tovisualize the results [40,41].The type of non-covalent interactions that exist between the nucleobases and thegemcitabine molecule was ascertained by means of local energy decomposition (LED)and dispersion interaction density (DID) studies. The optimized shape of B3LYP-D3/6-31G(d,p) was subjected to LED analysis. With the ‘TightPNO’ setting cc-pvdz basis set, weemployed the ORCA 6.0 program with the single-point domain-based local pair naturalorbital (DLPNO)-based coupled-cluster technique with single, double, and triple excitations(CCSD(T)) computation [42].Molecules 2025, 30, 2732 4 of 213. Results and Discussion3.1. Structures and Geometric Parameters of the Gemcitabine–Nucleobase ComplexesTo know the best mode of interaction of the drug with the nucleobase pairs, the druggemcitabine was allowed to interact with the nucleobases in all possible modes by tak-ing account of the active electrophilic and nucleophilic sites of the drug gemcitabine andthe nucleobases [43,44]. In addition, we employed the artificial bee colony approach forglobal optimization to ascertain the potential global minimum [45]. Initial conformationalsampling revealed ≥16 distinct binding modes per nucleobase. In the Supporting Infor-mation in Figures S1–S5, we have provided all the geometries. The most stable geometryof the drug–nucleobase complexes is shown in Figure 1. The purine nucleobase adenine(A) binds with gemcitabine through its N9-H and N3 bonds. The low-lying isomers areshown in the Supporting Information in Figure S1 in their order with decreasing stability.The binding site of gemcitabine was found to vary along with the bond distance to theadenosine. However, all the structures were found to possess at least one intermolecularhydrogen bond. Furthermore, the planar structures have higher stability than those withstacked configurations. In the most stable geometry shown in Figure 1a, the intermolecularhydrogen bonds (HBs) amongst the A and drug molecule were 1.939 and 1.931 Å distancefor N3···H-N-g and N9-H···N-g bonding, respectively. It is interesting to notice that theintramolecular hydrogen bond in the gemcitabine drug between the carbonyl group andthe hydroxyl group of the sugar moiety remains unchanged during the complex formation,which indicates that neither the carbonyl group nor the sugar moiety is involved in thecomplex formation. Figure 1. The most stable complexes of gemcitabine with (a) adenine, (b) guanine, (c) cytosine,(d) thymine, and (e) uracil. Inter- and intramolecular H-bonds are marked with blue lines along withbond length in angstroms.The drug gemcitabine interacts with the guanine nucleobase in 18 different modes. Themost stable mode of interaction exists in planar geometry, and is shown in Figure 1b. In theabove geometry, we notice three intermolecular HBs with 1.757, 1.939, and 1.929 Å distancesfor C6=O6···H-N, N2-H···N-g, and the N1-H···N-g bonding, respectively. The pristine GCbase pair C6=O6···H-N distance was longer by 0.017 Å, while the N2-H···N and the N1-H···N distances are shorter by 0.011 and 0.023 Å, which implies the strength of N3···H-Nand N9-H···O=C is increased and reduced, respectively [18]. Thus, the introduction ofsugar units increases the intermolecular HBs distance. This was further supported by theMolecules 2025, 30, 2732 5 of 21increased intramolecular HBs distance from 1.962 to 1.991 Å in pristine gemcitabine andthe guanine–gemcitabine complex, respectively. The low-lying isomers for the guanine–gemcitabine complex are shown in Figure S2, in which most of the isomers have two orfewer intermolecular HBs or buckle angles. Previous studies of DNA base pairs haveshown that planar geometries have higher interaction energy than geometries with buckleangles [46].The most stable geometries for the pyrimidine nucleobase–gemcitabine complexes areshown in Figure 1c,d. The most stable geometries of the pyrimidine nucleobase–gemcitabineall have two intermolecular HBs and one intramolecular hydrogen bond in the sugar-likemoiety of gemcitabine. The intermolecular HB distance in these complexes varies bynucleobase. It is worth pointing out that the thymine and uracil nucleobases differ by amethyl group at the C5 carbon. Furthermore, marginal changes in the intramolecular bonddistance are noticed in the gemcitabine drug after complex formation. In the complexes ofthymine and uracil with gemcitabine, an elongation in the intermolecular distance by 0.011and 0.001 Å was observed. By contrast, in the case of the cysteine–gemcitabine complex, thedistance was shortened by 0.020 Å. Furthermore, we also notice a change in glycosidic bondangle in all these complexes. Thus, the intermolecular HBs modulate the intramolecularHBs in these complexes.3.2. EnergeticsThe binding energy and the interaction energy of the gemcitabine drug with thenucleobases are computed at the M06-2X/6-311+G(d,p) level for all studied complexes.The stable geometry, binding energy, and interaction energies are listed in Table 1. Thebinding energy is derived from the total energy difference between the optimized complexand optimized gemcitabine and the nucleobase. While the interaction energy is calculatedas the energy difference between the optimized complex and the single-point energy ofgemcitabine and the nucleobase in their complex geometry. Both the binding energy andthe interaction energy are negative, indicating that their formation is feasible [47,48].Table 1. Computed binding energy, interaction energy, deformation energy, and thermodynamicparameters for the formation of the gemcitabine–nucleobase complexes for stable complexes.NucleobaseBinding Energy *(kcal mol−1)Interaction Energy(kcal mol−1)Deformation Energy(kcal mol−1) ∆H(kcal mol−1)∆G(kcal mol−1)Gas Solution Gas Solution Nucleobase GemcitabineAdenine −18.83(−19.05) −10.79 −21.33 −5.28 0.44 0.91 −10.37 0.33Guanine −26.51(−27.26) −14.84 −31.25 −17.49 1.18 1.54 −14.39 −3.14Cytosine −21.82(−22.25) −11.69 −25.53 −13.44 1.04 1.19 −11.27 −0.61Thymine −19.93(−20.15) −11.56 −23.39 −13.27 0.94 1.14 −11.14 −0.34Uracil −20.68(−20.88) −11.63 −23.66 −13.35 0.84 0.89 −11.21 −0.55* Values in the parentheses are BSSE corrected values.The other complexes’ relative energies with respect to the most stable geometries arelisted in the Supporting Information in Figures S1–S5. Among the complexes, the bindingenergy is highest for the gemcitabine–guanine complex and least for the gemcitabine–thymine complex. The same trend has been noted in the calculated interaction energy forthe 5FU–guanine complexes [49]. The higher binding energy and interaction energy for theMolecules 2025, 30, 2732 6 of 21gemcitabine–guanine complex indicate that gemcitabine binds with the guanine base morefirmly than the other nucleobases. Furthermore, the order of stability was found to matchtheir Gibbs free energy of formation. Incorporation of counterpoise correction tends toincrease the binding energy, while the trend remains the same. The basis set superpositionerror contribution is found to be about 0.009–0.027 kcal mol−1 higher in energy. While thepercentage of change in contribution was found to change by a maximum of 2.74%, whichsuggests that the basis set used in the present study would be good enough to study thebiomolecules with multiple hydrogen bonds.A comparison between the binding and interaction energy shows that the interactionenergy is low for all complexes except for the gemcitabine–guanine complex. Furthermore,in the solution phase, the binding energy and the interaction energy are lower in the gasphase, likely due to the solvation of the intermolecular HBs, which results in the weakeningof the HBs, leading to deformation in the geometry [50]. The computed deformationenergy of the nucleobase and the gemcitabine drug during the complexation is provided inTable 1. It lies in the 0.44–1.54 kcal mol−1 range, which supports the weakening of the HBsduring solvation.The deformation energy is generally higher for the gemcitabine drug than for thenucleobases, probably due to the presence of flexible sugar units. Among the complexes, thedeformation energy is higher for the gemcitabine–guanine complex. The thermodynamicquantities for the interaction among gemcitabine and the nucleobases were calculated, andare shown in Table 1. It is apparent from the table that the thermodynamic parameters,such as enthalpy, were negative for all the studied complexes, while the Gibbs free energyis negative for all the complexes except the gemcitabine–adenine complex. The negativevalues for enthalpy and Gibbs free energy indicate that gemcitabine complexation withthe nucleobase is exothermic and spontaneous. Furthermore, the change in enthalpy andGibbs free energy of formation of the complexes are higher for the gemcitabine–guaninecomplex. The enthalpy of formation for the pyrimidine complexes is lower than thatof the purine complexes. In general, the change in enthalpy is higher for all complexes,indicating that the complex formation is an enthalpy-driven process. Thus, the gemcitabinedrug complexation with the nucleobases is more facile and spontaneous, except for thegemcitabine–adenine complex, which is an enthalpy-driven process. Earlier, Rezaei-Sametiand Borojeni reported that the cytosine and guanine bases pair has the highest enthalpyfor the formation of the 5-Fluorouracil–nucleobase complex [14]. In the present study,gemcitabine shows higher selectivity towards the guanine base, as it is a cytidine analogue.3.3. Frontier Molecular Orbitals and Conceptual DFTThe two most important orbitals that promote a chemical reaction are the lowest unoc-cupied molecular orbital (LUMO) and the highest occupied molecular orbital (HOMO). Theability to contribute an electron is represented by the HOMO orbital, whereas the abilityto receive the electron is determined by the LUMO [51–53]. Consequently, the electricalstability of the system can be predicted using the frontier molecular orbital (FMO) gap,specifically the HOMO–LUMO gap, in a compound. Additionally, Koopmans’ theorem canbe used to estimate the quantum chemical descriptors, such as chemical hardness (η), chemi-cal potential (µ), electrophilicity index (ω), and softness (σ). The following Equations (1)–(5)can be used to calculate the aforementioned quantum chemical descriptors:η = IP − EA ≈ ELUMO − EHOMO (1)µ = −χ = (IP + EA)/2 ≈ (EHOMO + ELUMO)/2 (2)ω = µ2/2η ≈ (EHOMO + ELUMO)2/(4 (ELUMO − EHOMO)) (3)Molecules 2025, 30, 2732 7 of 21σ = 1/(IP − EA) = 1/η ≈ 1/(ELUMO − EHOMO) (4)∆N = −µ/η (5)The computed chemical descriptors for the gemcitabine–nucleobase complexes areprovided in Table 2. The energy gap for the nucleobases was in the range of 6.92–8.12 eV,while gemcitabine has an Egap of 7.82 eV. The drug complex lies in the range of 5.96–7.15 eV,which is less than the pristine nucleobase and the drug gemcitabine, which reveals that drugcomplexes are more chemically reactive towards further reaction with another nucleobaseor the gemcitabine drug, which is less likely. The chemical hardness values are moreprominent for the complexes than the nucleobases, and follow the trend gemcitabine–uracil> gemcitabine–thymine > gemcitabine–cytosine > gemcitabine–adenine > gemcitabine–adenine. Consequently, the gemcitabine–uracil complexes are less likely to undergo electrontransfer, and are less sensitive than the nucleobases. Also, the term softness supports theprior conclusion on reactivity.Table 2. The calculated band gap (Egap (in eV)) and the reactive descriptors (in eV) for gemcitabine,the nucleobase, and the gemcitabine–nucleobase complexes at the M06-2X-D3/6-311+G(d,p) levelof theory.System Egap µ η ω σ ∆NGemcitabine 7.827 4.558 7.827 1.327 0.128 −0.582Adenine 7.670 3.787 7.670 0.935 0.130 −0.494Guanine 6.929 3.791 6.929 1.037 0.144 −0.547Cytosine 7.772 4.150 7.772 1.108 0.129 −0.534Thymine 7.862 4.394 7.862 1.228 0.127 −0.559Uracil 8.129 4.641 8.129 1.325 0.123 −0.571Gemcitabine–Adenine 6.658 4.132 6.658 1.282 0.150 −0.621Gemcitabine–Guanine 5.959 3.888 5.959 1.268 0.168 −0.652Gemcitabine–Cytosine 7.025 4.234 7.025 1.276 0.142 −0.603Gemcitabine–Thymine 6.845 4.390 6.845 1.408 0.146 −0.641Gemcitabine–Uracil 7.150 4.579 7.150 1.466 0.140 −0.640The electrophilicity index (ω) is a positive measure, and high values are features ofmost electrophilic systems [54,55] The electrophilicity value of adenine is the lowest amongthe nucleobases studied. The gemcitabine–thymine and gemcitabine–uracil complexeshave the highest values in the drug–nucleobase complexes. Thus, one can anticipate thatthe above complexes are more electrophilic than the other drug–nucleobase complexes.The chemical potential (µ) in a system designates the escaping tendency of an electroncloud [56]. When complex formation occurs, the electrons flow from the higher chemicalpotential to the lower species until the electronic chemical potential becomes equal. The freedrug gemcitabine has a chemical potential of 4.55 eV, while the nucleobases have chemicalpotential values in the range of 3.78–4.64 eV.Compared to the pristine nucleobase and the gemcitabine drugs, the drug–nucleobasecomplexes have lower chemical potential, making the latter less proficient at attracting elec-trons. We have calculated the charge transfer based on electrophilicity using Equation (6)to ascertain the direction of charge transfer.ECT = (∆Ngemcitabine) − (∆Nnucleobase) (6)In the case of ECT < 0, gemcitabine acts as an electron donor; if ECT > 0, then gemc-itabine acts as an electron acceptor. For the gemcitabine–nucleobase complexes, all valuesare negative, indicating that gemcitabine gains electrons during complex formation. Theabove results are further proved by the existence of the HOMO and LUMO orbitals ob-Molecules 2025, 30, 2732 8 of 21served on the complex. The HOMO and LUMO orbitals of the gemcitabine–nucleobasecomplexes and the pristine nucleobase and gemcitabine are shown in Figure S6. It can beseen that, in all the studied complexes, the HOMO is concentrated on the surface of thenucleobases, and the LUMO is localized on the surface of the gemcitabine drug. Thus, onecan anticipate that the nucleobase of the gemcitabine–nucleobase complex can undergoan electrophilic attack, and the gemcitabine drug can undergo a nucleophilic attack whensubjected to further chemical reactions.3.4. Nature of InteractionsThe interactions between the molecules help one predict and manipulate the behaviorof the substances at the molecular level. To understand the nature of the intermolec-ular interactions between the drug and the nucleobases, we performed a quantitativemolecular electrostatic potential (MESP) analysis [57], an atoms-in-molecules (AIM) analy-sis [58], and a non-covalent interaction–reduced density gradient (NCI–RDG) analysis [59].Hydrogen bonding interactions have been extensively demonstrated to enhance cyclic4n + 2π electron delocalization and augment aromaticity [60]. Therefore, we computed theharmonic oscillator model of aromaticity (HOMA) and the nucleus independent chemicalshifts (NICS) analyses for the gemcitabine–nucleobases and compared them with the drug–nucleobase complexes to comprehend the hydrogen bonding between the drug gemcitabineand the nucleobases.3.4.1. π-Electron Delocalization Indicators and Molecular Electrostatic PotentialIt is well known that hydrogen bonds in DNA base pairs are strongly affected by π-electron delocalization (π-ED), which is not directly measurable but can be evaluated withother molecular properties, including energetic, structural, and magnetic properties [61]. Inthe present study, we have calculated the harmonic oscillator model of aromaticity (HOMA)and the nucleus independent chemical shifts (NICS) analyses for the nucleobases and thedrug gemcitabine and related them with the gemcitabine–nucleobase complexes.The computed HOMA, NICS(0), and NICS(1) values are provided in Table 3. Ingeneral, the HOMA value of gemcitabine was found to increase upon complexation withthe nucleobases. In the gemcitabine–nucleobase complex, the HOMA readings of both 6-membered and 5-membered rings were observed to rise. However, the highest increase wasobserved on the 6-membered ring of the gemcitabine–guanine complex. In the pyrimidinebase complex, the cytosine ring showed a significant increase in the HOMA value aftercomplexation with cytosine, which is determined by the observed binding energy valuesof the system. Previous studies have reported that binding energies of complexes areproposed to increase the HOMA values [62].The NICS is an essential parameter for assessing the π-ED of aromatic and pseudo-aromatic molecules. The magnitude of electron delocalization is directly related to thetotal value of the NICS(1). It is well documented that, when the NICS(1) value ofa molecule increases, the π-ED increases [16]. The computed NICS(1) values on thegemcitabine–nucleobase complexes show a substantial increase compared to the pristine nu-cleobases. The increase in the NICS(1) values was highest for the gemcitabine–guanine andgemcitabine–cytosine complexes, which are by the binding energy and the HOMA values.Molecules 2025, 30, 2732 9 of 21Table 3. Computed HOMA and NICS aromaticity indices for gemcitabine, the nucleobases, and themost stable gemcitabine–nucleobase complexes.ComplexGemcitabine A/G C/T/UHOMA NICS(0) NICS(1) HOMA(6/5)NICS(0)(6/5)NICS(1)(6/5) * HOMA NICS(0) NICS(1)Gemcitabine 0.676 −1.488 −3.668 - - - - - -Adenine - - - 0.981/0.881−6.304/−11.455−8.946/−10.199 - - -Guanine - - - 0.700/0.883−2.628/−11.481−3.719/−9.721 - - -Cytosine - - - - - - 0.672 −0.541 −3.010Thymine - - - - - - 0.466 −1.079 −2.226Uracil - - - - - - 0.499 −0.661 −1.814Gemcitabine–Adenine 0.693 −1.376 −3.254 0.983/0.912−6.346/−11.727−8.703/−10.675 - - -Gemcitabine–Guanine 0.703 −1.568 −3.396 0.824/0.875−3.166/−11.001−4.141/−9.385 - - -Gemcitabine–Cytosine 0.693 −1.549 −3.377 - - - 0.794 −1.441 −3.833Gemcitabine–Thymine 0.701 −1.616 −3.372 - - - 0.540 −1.315 −2.422Gemcitabine–Uracil 0.706 −1.658 −3.480 - - - 0.574 −0.951 −2.187* represents values on the 6- and 5-membered rings.The molecular electrostatic potential (MESP) helps us to understand the attractionbetween individual molecules, as charge transfer/electrostatic interactions are known tobe one of the significant driving forces during complex formation [63,64]. In Figure 2, thequantitative electrostatic potentials for the bases, gemcitabine, and gemcitabine–nucleobasecomplexes are provided and mapped on a 0.001 electrons/Bohr3 isodensity surface. Inthe above figures, the red color represents the electrophilic area, while the blue and greencolor depicts nucleophilic and neutral regions in the molecules. These electrophilic andnucleophilic regions are designated as Vs,max and Vs,min [65]. In the nucleobases and thegemcitabine drug, the Vs,max is observed on the hydrogen atom of the pyrrole-like ringand the hydrogen atoms of the amino group. As gemcitabine resembles that of cytosine,the MESP regions of gemcitabine are similar to those of cytosine; however, the Vs,max, andVs,min values differ a little due to the presence of an extra sugar moiety in gemcitabine.In all the nucleobases and the gemcitabine drugs, the Vs,max are observed on the iminenitrogen and the oxygen atom of the carbonyl group.It is commonly known that the system becomes unstable when sections of the same po-tential are brought together [66]. A directional interaction with opposite potentials leads tothe establishment of the most stable complex. Owing to the directional interaction, the pos-itive potential regions will interact with the negative region, leading to a decrease/increasein positive/negative potential values. The MESP of the gemcitabine–nucleobase complexesare visualized in Figure 2g–k. In the gemcitabine–adenine complex, the most positive valueobserved on the pyrrole-like hydrogen atom of adenine matches with the most negative po-tential value observed on the imide nitrogen atom of the gemcitabine drug. Furthermore, inthe complex, the Vs,min value observed on the carbonyl oxygen was reduced to −26.1 from−45.6 kcal mol−1 in pristine gemcitabine. This indicates the carbonyl group’s participationin complex formation, leading to an electron density shift. A similar supposition can bearrived at while observing the MESP potential values of all other stable complexes. Thus,during complex formation, a directional interaction along with an electron density shiftoccurs to help form a stable complex.Molecules 2025, 30, 2732 10 of 21Figure 2. Molecular electrostatic surface potentials mapped on the corresponding 0.001 a.u. elec-tron density isosurface: (a) A, (b) G, (c) C, (d) T, (e) U, (f) gemcitabine, (g) gemcitabine–A,(h) gemcitabine–G, (i) gemcitabine–C, (j) gemcitabine–T and (k) gemcitabine–U complexes.3.4.2. Atom in Molecule Analysis and Non-Covalent Interaction AnalysisThe quantum theory of atoms in molecules (QTAIM) is used to examine the natureof interactions in many molecular systems and to categorize and understand bondinginteractions [67]. The QTAIM is also used to compute the individual hydrogen bondstrength if the topological parameter potential energy at the bond critical point (BCP) isknown [68]. The various topological parameters at the BCPs, namely electron density(ρ(r)), its Laplacian (∇2ρ(r)), and the total electron density (H(r)), allow us to classify theinteractions. For hydrogen bonds, the Laplacian (∇2ρ(r)) is positive, the potential energy(V(r)) is negative, and the kinetic energy density G(r) and H(r) are positive values. If abond has positive ∇2ρ(r) and H(r), a weak noncovalent hydrogen bond is proposed. Thetomographs obtained for all of the stable gemcitabine–nucleobase complexes are shown inFigure 3.It can be observed that the total number of intramolecular BCPs on gemcitabineincreases from two to three in both pristine gemcitabine and the gemcitabine–nucleobasecomplexes. During the complexation, an additional intramolecular BCP between thefluorine and oxygen atoms arises due to the structural deformation of gemcitabine duringMolecules 2025, 30, 2732 11 of 21complexation. It is worth pointing out that the computed deformation energies are higheron gemcitabine than on the nucleobase.Figure 3. Molecular topography analysis for (a) gemcitabine–adenine, (b) gemcitabine–guanine,(c) gemcitabine–cytosine, (d) gemcitabine–thymine, (e) gemcitabine–uracil, and (f) gemcitabine.The topological values for the intermolecular interactions between gemcitabine andthe nucleobases are provided in Table 4. The ρ(r) values for all the intermolecular BCPswere in the range of 0.004–0.037 a.u., which is close to the value proposed for hydrogenbonding. The ρ(r) value is highest for the N-H···O=C type of interaction and lowest forthe C=O···H-C interaction. The ∇2ρ(r) values are positive for all the intermolecular BCPs,while the H(r) values are negative for most bonds. The negative H(r) values proposethe supremacy of potential energy over kinetic energy. H(r) with negative values forintermolecular BCPs shows a partial covalent character.The gemcitabine molecule shows intramolecular BCPs, with the oxygen atom ofthe carbonyl in pyrimidine, bonding with the hydroxyl, alkyl, and fluoro groups. Theintramolecular BCPs have ρ(r) values of 0.009–0.023 a.u. Among the intramolecular bondsin gemcitabine, the hydroxyl–carbonyl oxygen interaction has the highest ρ(r) values andremains significantly unchanged after complexation with the nucleobases. The ρ(r) valueson the alkyl hydrogen–carbonyl oxygen in complexes are lower compared to the freegemcitabine molecule. The intramolecular F···O bond has ρ(r) values in the range of0.009–0.014 a.u., which is closer to the value reported for the Cl···O intramolecular bond. Itis interesting to see that the ρ(r) values of F···O are higher for the complexes with largerbinding energy, the gemcitabine–guanine complex, while the other intramolecular bondsremain nearly unchanged. The existence of weak interactions between electronegativeatoms is not uncommon. A comparison between the inter- and intramolecular ρ(r) valuesshows that the latter are weak. The higher deformation energy of gemcitabine and higherρ(r) values of F···O in the gemcitabine–guanine complex suggest that the alkyl chain playsa significant role in the complex formation process.Molecules 2025, 30, 2732 12 of 21Table 4. Topological parameters of intermolecular bonds obtained from the AIM analysis for thestable gemcitabine–nucleobase complexes. The values of topological parameters are provided in a.u.and bond length in Å; The value of EHB is provided kcal/mol.BCP BondLength ρ ∇2ρ λ1 λ2 λ3 V G H EHB EllipticityGemcitabine–Adenine complexHB1 1.939 0.02981 0.08956 −0.04274 −0.04004 0.17233 −0.02203 0.02221 −0.00018 6.91 0.07HB2 1.931 0.02988 0.09520 −0.04237 −0.04013 0.17770 −0.02290 0.02335 −0.00045 7.18 0.06Gemcitabine–Guanine complexHB1 1.757 0.03704 0.13148 −0.05969 −0.05781 0.24897 −0.03392 0.03340 0.00053 10.64 0.03HB2 1.939 0.02992 0.09284 −0.04263 −0.03997 0.17543 −0.02244 0.02282 −0.00039 7.04 0.07HB3 1.929 0.02475 0.09970 −0.03412 −0.03222 0.16603 −0.01977 0.02235 −0.00258 6.20 0.06Gemcitabine–Cytosine complexHB1 1.778 0.03512 0.12700 −0.05529 −0.05339 0.23568 −0.03151 0.03163 −0.00012 9.88 0.04HB2 1.834 0.03815 0.10525 −0.06008 −0.05632 0.22164 −0.03214 0.02922 0.00291 10.08 0.07Gemcitabine–Thymine complexHB1 1.886 0.02698 0.10566 −0.03861 −0.03682 0.18109 −0.02193 0.02417 −0.00224 6.87 0.05HB2 1.849 0.03667 0.10455 −0.05674 −0.05278 0.21407 −0.03045 0.02829 0.00216 9.55 0.08Gemcitabine–Uracil complexHB1 1.879 0.02719 0.10691 −0.03916 −0.03753 0.18360 −0.02224 0.02448 −0.00225 6.98 0.04HB2 1.828 0.03849 0.10526 −0.06085 −0.05704 0.22315 −0.03259 0.02945 0.00314 10.23 0.07The ratio of |G/V| value helps one to decide the electrostatic or covalent nature of theHB. The computed values for the intermolecular HBs lie in the range of 0.90–1.26. Thesevalues confirm that the bonds have a partially covalent character. The ellipticity at the BCPis used as a measure of bond stability as it provides the extent of charge accumulation at theBCP. A BCP with a higher ellipticity suggests that the bond can undergo rupture. Hence,electrostatic interactions with the van der Waals interactions have higher ellipticity thanhydrogen bonds [69]. From Table 4, it is evident that the C=O···H-C6 bonds have higherellipticity values and hence are weaker in strength than the N3···H-N and the NH···N3bonds. Computing the individual hydrogen bond strength is intimidating; however, usingthe AIM analysis and the Espinosa–Molins–Lecomte (EML) semi-empirical equation, wecomputed the hydrogen bond energy (EHB), which is approximately half the value ofthe potential energy density at the BCPs. It should be pointed out that potential energydensity describes the attractive component of the total energy density at the BCP. Thecomputed intermolecular EHB values in kcal mol−1 are provided in Table 4. The NH···Nbond strength ranges from 6.87 to 10.64 kcal mol−1. An NCI analysis is a significant methodfor identifying non-covalent interactions within a molecule [70]. It involves hydrogenbonding, steric effects, and the presence of van der Waals forces. The attractive, repulsive,and van der Waals inter- and intramolecular interactions are deciphered from the sign ofλ2. In Figure 4a–e, the RDG isosurface for the stable drug nucleobases is provided. Thered regions depict the repulsive interactions. The green and brown regions portray theattractive van der Waals interaction, with the green region being more substantial than thebrown one. The blue color patches confirm the strong hydrogen bond interactions. Thepresence of brown regions is due to the rings possessing π-electrons. In the pyrimidinebases of the drug, the presence of additional H-bonding is confirmed by the green and blueisosurfaces between the H and carbonyl O atoms of the gemcitabine drug. The occurrenceMolecules 2025, 30, 2732 13 of 21of several regions with red, green, and brown colors implies the coexistence of repulsiveand attractive interactions, which can lead to anti-cooperative effects [71].Figure 4. NCI plots of (a) gemcitabine-adenine (b) gemcitabine-guanine (c) gemcitabine-cytosine(d) gemcitabine-thymine and (e) gemcitabine-uracil complex and RDG of (f) gemcitabine-adenine(g) gemcitabine-guanine (h) gemcitabine-cytosine (i) gemcitabine-thymine and (j) gemcitabine-uracilcomplex depicting non-covalent interaction for the most stable complexes.Figure 4f–j, depicts the plots of the RDG vs. the electron density multiplied by thesign of λ2 for the stable drug–nucleobase complexes. For the pristine nucleobases andthe gemcitabine drug, the same is shown in the Supporting Information Figure S7. TheMolecules 2025, 30, 2732 14 of 21spike regions corresponding to the steric effect exhibit no noticeable change in the pattern,confirming that there is no apparent change in the steric effect after the adsorption ofgemcitabine to the nucleobase. However, in the spike regions involving the H-bondinginteractions (−0.010–−0.025 a.u.), many new spikes were produced along with an increasein their strength. This supports that the H-bonding interactions are stronger and play adominant role in stabilizing the complexes [72].3.4.3. Dispersion Interaction Density (DID) and Local Energy Decomposition (LED) AnalysesA dispersion interaction density (DID) analysis helps for visualizing weak van derWaals (vdW) interactions in complexes [73]. The DID values are projected onto an isodensitysurface, obtained from the total molecular density, and are used as a simple technique forvisualizing the most critical contacts that contribute to the van der Waals (vdW) interactionof any complex. The DID profile obtained for the complexes is shown in Figure 5, wherethe red zones indicate strong dispersion interactions and the blue zones indicate weakerdispersion interactions.Figure 5. Dispersion interaction density profile for the complexes (a) gemcitabine–adenine,(b) gemcitabine–guanine, (c) gemcitabine–cytosine, (d) gemcitabine–thymine, and (e) gemcitabine–uracil.The DID profile revealed that attractive dispersion interactions occur along the inter-face between gemcitabine and the nucleobases, with the strongest interactions observedbetween the hydrogen atoms of the amine group of gemcitabine and the nitrogen and car-bonyl oxygen atoms of the nucleobases. In addition, the imide group of the nucleobases hasthe strongest interaction with the nitrogen atom of gemcitabine. Furthermore, the presenceof green zones indicates weak attractive dispersion interactions between the alcohol groupof gemcitabine and the C-H/N-H group of the nucleobases. A comparison between thepurine and pyrimidine nucleobases reveals that the red zones are more concentrated in thepurine bases; whereas, in the pyrimidine bases, they are more dispersed.In order to gain a quantitative measurement on the nature of interactions in thesecomplexes, we carried out local energy decomposition (LED) calculations. The LED analysisprovides a breakdown of the domain-based local pair natural orbital CCSD(T) [DLPNO–CCSD(T)] energy into additive contributions representing the interaction between pairs ofMolecules 2025, 30, 2732 15 of 21user-defined fragments [74,75]. In LED, the total binding energy between the fragmentscan be written as in Equation (7):∆E = ∆Eint + ∆Egeo-prep (7)where ∆Eint is the interaction energy between the fragments considered, and ∆Egeo-prep isthe deformation energy.The ∆Eint can be decomposed into Hartree–Fock (∆EHFint ) and correlation contributions(∆ECint) as follows:∆Eint = ∆EHFint +∆ECint (8)The Hartree–Fock and correlation contributions can be further disintegrated as follows:∆EHFint = ∆EHFel−prep+∆E f ragselstat+∆E f ragsexch (9)∆ECint= ∆EC−CCSDnon−dispersion+∆EC−CCSDdis +∆EC−(T)int (10)In Equation (9), ∆EHFel−prep is the Hartree–Fock contribution to the electronic preparationof each fragment in the complex, and is always a positive term. The terms ∆E f ragselstat + ∆E f ragsexchdefine the electrostatic and exchange interactions, while ∆EC−CCSDnon−dispersion is a summationterm of weak pair and strong pair dispersion interactions between fragments, and ∆EC−(T)intis the triples correction to the interaction energy, which includes the electronic preparationand interfragment interaction energies.The DLPNO–CCSD(T)–LED calculations for the gemcitabine–nucleobase complexeswere performed using the cc-pVDZ basis set on the B3LYP–D3/6–31G(d,p) optimizedgeometry. The results obtained in the LED calculations are presented in Table 5. The∆Egeo-prep values for the complexes follow the same trend as the computed summationdeformation energy of the nucleobase and gemcitabine [76]. The ∆Egeo-prep values aremore significant for the gemcitabine–guanine complex, and least for the gemcitabine–adenine complex. Expectedly, the ∆Egeo-prep values are far smaller than the ∆Eint values.Furthermore, the ∆Eint, which is interaction energies computed at the HF level, is far moresignificant than the correlation interaction contribution obtained for the complexes. It isevident from Table 5 that the sum of the attractive electrostatic preparation energy ∆Ere felstatand the exchange energy ∆Ere fexch is almost fully compensated by its repulsive electronicpreparation term ∆Ere fel−prep. Dispersion contributes the least to interaction energy, andelectrostatic contributes the most. Furthermore, the larger electrostatic term ∆Ere felstat of thegemcitabine–guanine complex, accounts for its higher stability among the complexes.Table 5. DLPNO–CCSD(T) based local energy decomposition of the binding energy of gemcitabineon the nucleobases into various energy components (kcal mol−1).System ∆Eint ∆Egeo-prep ∆Erefel−prep ∆Erefelstat ∆Erefexch ∆EC−CCSDnon−disp ∆EC−CCSDdisp ∆EC−(T)intGemcitabine–Adenine −17.68 0.92 76.44 −81.51 −12.62 5.45 −5.33 −0.93Gemcitabine–Guanine −27.59 2.00 109.34 −119.48 −17.44 7.53 −6.50 −1.03Gemcitabine–Cytosine −22.59 1.57 93.33 −101.03 −14.89 6.75 −5.38 −0.82Gemcitabine–Thymine −19.96 1.44 79.74 −86.60 −13.11 6.60 −5.40 −0.81Gemcitabine–Uracil −20.65 1.27 82.13 −89.41 −13.36 6.63 −5.29 −0.81Molecules 2025, 30, 2732 16 of 213.5. Role of Intramolecular Hydrogen Bond of Gemcitabine in Stabilizing the ComplexTo understand the role of the flexible sugar moiety of gemcitabine in the complexformation, the intramolecular interactions that exist between the hydroxyl group/alkylhydrogen/fluorine and the carbonyl oxygen were studied using NBO analysis and QTAIMtheory. The selected interactions are listed in Table 6. The gemcitabine molecule is stabilizedby the transfer of a lone pair of electrons located on the oxygen atom of the carbonyl groupto the σ orbital of the hydroxyl group with a stabilization energy of 3.03 kcal mol−1. Thecomputed EML-based hydrogen bond energy for the interaction was 5.46 kcal mol−1. Inaddition, the carbonyl oxygen atom donates electrons to the σ orbital of the alkyl hydrogenatom, stabilizing by 0.53 kcal mol−1, while its EHB was 3.08 kcal mol−1. Furthermore, thecomputed EHB values are nearly equal in strength to that of a CH···O bond, and the com-bined values are found to be nearly equivalent to the conventional intramolecular hydrogenbond observed between the –OH and the carbonyl oxygen atom. Thus, charge transferalone does not act as a predominant factor in deciding the strength of the intramolecularbonds. The fixable sugar moiety undergoes a more significant structural deformation, andthis cooperative effect improves the intramolecular bond strength.Table 6. NBO analysis of the drug gemcitabine–nucleobase complexes and the hydrogen bond energy(EHB) computed using QTAIM theory.Donor NBO (i) Acceptor NBO (j) E(2) kcal mol−1 F(I,j) a.u. EHBGEMLP(1) O29 BD*(1) O22-H23 3.03 0.06 5.46LP(1) O29 BD*(1) C15-H16 0.53 0.02 3.10Gemcitabine–Adenine complexLP(1) O29 BD*(1) O22-H23 2.92 0.06 5.54LP(1) O29 BD*(1) C15-H16 0.50 0.22 2.97Gemcitabine–Guanine complexLP(1) O29 RY*(1) H23 0.84 0.04 5.39LP(1) O29 RY*(1) H23 0.60 0.04 -LP(1) O29 BD*(1) C15-H16 - - 2.98Gemcitabine–Cytosine complexLP(1) O29 BD*(1) O22-H23 3.21 0.06 5.92LP(1) O29 BD*(1) C15-H16 - - 3.01Gemcitabine–Thymine complexLP(1) O29 BD*(1) O22-H23 2.69 0.05 5.39LP(1) O29 BD*(1) C15-H16 - - 2.74Gemcitabine–Uracil complexBD*(2) C28-29 BD*(1) O22-H23 0.69 0.04 5.62LP(1) O29 BD*(1) C15-H16 - - 2.904. ConclusionsIn summary, the interaction between the drug gemcitabine and the nucleobases hasbeen explored by means of the dispersion-corrected density functional theory study. Thebinding site of gemcitabine was found to vary along with the nucleobase, and all the moststable structures were found to possess at least two intermolecular hydrogen bonds. Inthe pyrimidine nucleobase–gemcitabine complexes, a change in glycosidic bond angle wasobserved. The binding energy and the interaction energy are all negative and are maximumfor the gemcitabine–guanine complex, and minimum for the gemcitabine–thymine complex.Molecules 2025, 30, 2732 17 of 21Due to the presence of flexible sugar units, the deformation energy is higher for thegemcitabine drug than for the nucleobases. The enthalpy for the formation of all thecomplexes was negative, while the Gibbs free energy was positive for the gemcitabine–adenine complex. The complex formation was found to be an enthalpy-driven process.Pyrimidine nucleobases have a lower enthalpy of formation than purine nucleobases. Theorder of stability was found to match their Gibbs free energy of formation.The nature of the interaction between the drug gemcitabine and the nucleobases wasidentified using the MESP, AIM, HOMA, NICS, and NCI–RDG analyses. In the gemcitabine–nucleobase complex, the HOMA values of 6- and 5-membered rings increased comparedto the pristine nucleobase. The highest increase was observed on the 6-membered ringof the gemcitabine–guanine complex. The computed NICS(1) values on the gemcitabine–nucleobase complexes showed a substantial increase compared to the pristine nucleobases.In the gemcitabine–nucleobase complex, the most positive value was observed on thepyrrole-like hydrogen atom of the nucleobase matches, with the most negative potentialvalue observed on the imide nitrogen atom of the gemcitabine drug. Furthermore, in thecomplex, the Vs,min value observed on the carbonyl oxygen was reduced substantiallycompared to pristine gemcitabine. Thus, a directional interaction and electron density shiftoccurred during the complex formation to help form a stable complex. The QTAIM analysisshowed the presence of two/three intramolecular BCPs in the gemcitabine–nucleobasecomplexes, and an additional intramolecular BCP between the fluorine and the oxygenatoms arose due to the structural deformation of gemcitabine during complexation. Theintermolecular hydrogen bonds showed a partial covalent character. The computed bondenergy showed that the intermolecular NH···N bond is stronger than the other bonds.Among the intramolecular bonds, the electrostatic interaction between F and carbonyloxygen was the weakest bond observed. The above results are supported by the NCI–RDG analysis. The intermolecular bond strength between the gemcitabine–nucleobasecomplexes is about two-times higher in energy than the intramolecular bond strengthobserved in the gemcitabine drug. Though the present work cannot be used to understandthe drug interactions with DNA and enzymes, this work provides a new insight into theintercalation mode of the drug gemcitabine with the nucleobases, and it will be helpful toreveal the mechanism of drug action and the further design of new anticancer drugs.Supplementary Materials: The following supporting information can be downloaded at https://www.mdpi.com/article/10.3390/molecules30132732/s1: Figures S1–S5: Optimized geometries of thelow-lying isomers of gemcitabine–adenine, gemcitabine–guanine, gemcitabine–cytosine, gemcitabine–thymine. and gemcitabine–uracil; Figure S6: HOMO and LUMO orbitals of the nucleobase, thegemcitabine drug, and the most stable gemcitabine–nucleobase complexes; Figure S7: NCI and RDGplots of the nucleobases and gemcitabine.Author Contributions: N.S.V.: conceptualization, writing—original draft, methodology, formalanalysis, investigation. A.S.: methodology, formal analysis. R.S.: resources, methodology, formalanalysis. All authors have read and agreed to the published version of the manuscript.Funding: This research received no external funding.Institutional Review Board Statement: Not applicable.Informed Consent Statement: Not applicable.Data Availability Statement: Data are contained within the article and Supplementary Materials.Acknowledgments: The author thank the staff of the Center for Computational Materials Science,Institute for Materials Research, Tohoku University, and the supercomputer resources through theHPCI System Research Project (Project ID: hp200040).https://www.mdpi.com/article/10.3390/molecules30132732/s1https://www.mdpi.com/article/10.3390/molecules30132732/s1Molecules 2025, 30, 2732 18 of 21Conflicts of Interest: The authors declare no conflicts of interest.AbbreviationsThe following abbreviations are used in this manuscript:CCSD(T) Coupled Cluster Singles Doubles with perturbative TriplesDFT Density Functional TheoryDLPNO Domain-Based Local Pair Natural OrbitalHOMA Harmonic Oscillator Measure of AromaticityNICS Nucleus-Independent Chemical ShiftNCI–RDG Non-Covalent Interaction–Reduced Density GradientQTAIM Quantum Theory of Atoms in MoleculesReferences1. Bray, F.; Laversanne, M.; Sung, H.; Ferlay, J.; Siegel, R.L.; Soerjomataram, I.; Jemal, A. Global cancer statistics 2022: GLOBOCANestimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA. Cancer J. Clin. 2024, 74, 229–263. [CrossRef][PubMed]2. Arnold, M.; Morgan, E.; Rumgay, H.; Mafra, A.; Singh, D.; Laversanne, M.; Vignat, J.; Gralow, J.R.; Cardoso, F.; Siesling, S.; et al.Current and future burden of breast cancer: Global statistics for 2020 and 2040. Breast 2022, 66, 15–23. [CrossRef] [PubMed]3. Gottesman, M.M. Mechanisms of cancer drug resistance. Annu. Rev. Med. 2022, 53, 615–627. [CrossRef] [PubMed]4. Kroemer, G.; Galluzzi, L.; Kepp, O.; Zitvogel, L. Immunogenic cell death in cancer therapy. Annu. Rev. Immunol. 2013, 31, 51–72.[CrossRef]5. Toschi, L.; Finocchiaro, G.; Bartolini, S.; Cappuzzo, F. Role of gemcitabine in cancer therapy. Future Oncol. 2005, 1, 7–17. [CrossRef]6. Cha, H.-M.; Kim, U.I.; Ahn, S.B.; Lee, M.K.; Lee, H.; Bang, H.; Jang, Y.; Kim, S.S.; Bae, M.A.; Kim, K.; et al. Evaluation of AntiviralActivity of Gemcitabine Derivatives against Influenza Virus and Severe Acute Respiratory Syndrome Coronavirus 2. ACS. Infect.Dis. 2023, 9, 1033–1045. [CrossRef]7. Larson, A.C.; Doty, K.R.; Solheim, J.C. The double life of a chemotherapy drug: Immunomodulatory functions of gemcitabine incancer. Cancer Med. 2024, 13, e7287. [CrossRef]8. Alsarayreh, N.; Abdelghany, S.; Alqudah, D.; Dbuarqoub, D.; Alshaer, W. Gemcitabine-loaded niosomes: Optimization, character-ization, and in vitro efficacy against invasive malignancies. J. Drug. Deliv. Sci. Technol. 2024, 95, 105617. [CrossRef]9. Hruba, L.; Das, V.; Hajduch, M.; Dzubak, P. Nucleoside-based anticancer drugs: Mechanism of action and drug resistance.Biochem. Pharmacol. 2023, 215, 115741. [CrossRef]10. Chai, X.; Meng, Y.; Ge, W.; Wang, J.; Li, F.; Wang, X.J.; Wang, X. A novel synthesized prodrug of gemcitabine based on oxygen-freeradical sensitivity inhibited the growth of lung cancer cells. J. Biomed. Res. 2023, 37, 355–366. [CrossRef]11. Natu, J.; Nagaraju, G.P. Gemcitabine effects on tumor microenvironment of pancreatic ductal adenocarcinoma: Special focus onresistance mechanisms and metronomic therapies. Cancer Lett. 2023, 573, 216382. [CrossRef] [PubMed]12. Abbaspour, A.; Dehghani, M.; Setayesh, M.; Tavakkoli, M.; Rostamipour, H.A.; Ghorbani, M.; Ramzi, M.; Omidvari, S.; Moosavi,F.; Firuzi, O. Cytidine deaminase enzyme activity is a predictive biomarker in gemcitabine-treated cancer patients. CancerChemother. Pharmacol. 2023, 92, 475–483. [CrossRef] [PubMed]13. Alniss, H.Y.; Al-Jubeh, H.M.; Msallam, Y.A.; Siddiqui, R.; Makhlouf, Z.; Ravi, A.; Hamdy, R.; Soliman, S.S.M.; Khan, N.A.Structure-based drug design of DNA minor groove binders and evaluation of their antibacterial and anticancer properties. Eur. J.Med. Chem. 2024, 271, 116440. [CrossRef] [PubMed]14. Rezaei-Sameti, M.; Borojeni, Z.I. Interaction of 5-fluorouracil anticancer drug with nucleobases: Insight from DFT, TD-DFT, andAIM calculations. J. Biomol. Struct. Dyn. 2023, 41, 5882–5893. [CrossRef]15. Hokmabady, L.; Raissi, H.; Khanmohammadi, A. Interactions of the 5-fluorouracil anticancer drug with DNA pyrimidine bases:A detailed computational approach. Struct. Chem. 2016, 27, 487–504. [CrossRef]16. Nakhaei, E.; Nowroozi, A.; Ravari, F. The influence of 5-fluorouracil anticancer drug on the DNA base pairs; a quantum chemicalstudy. J. Biomol. Struct. Dyn. 2018, 37, 1–19. [CrossRef]17. Qiu, Z.M.; Wang, G.L.; Wang, H.L.; Xi, H.P.; Hou, D.N. MP2 study on the hydrogen-bonding interaction between 5-fluorouraciland DNA bases: A, C, G, T. Struct. Chem. 2014, 25, 1465–1474. [CrossRef]18. Venkataramanan, N.S.; Suvitha, A.; Sahara, R. Unveiling the Intermolecular Interactions between Drug 5-Fluorouracil andWatson–Crick/Hoogsteen Base Pairs: A Computational Analysis. ACS Omega 2024, 9, 24831–24844. [CrossRef]19. Mazumdar, P.; Kashyap, A.; Choudhury, D. Investigation of hydrogen bonding in small nucleobases using DFT, AIM, NCI andNBO technique. Comput. Theor. Chem. 2023, 1226, 114188. [CrossRef]https://doi.org/10.3322/caac.21834https://www.ncbi.nlm.nih.gov/pubmed/38572751https://doi.org/10.1016/j.breast.2022.08.010https://www.ncbi.nlm.nih.gov/pubmed/36084384https://doi.org/10.1146/annurev.med.53.082901.103929https://www.ncbi.nlm.nih.gov/pubmed/11818492https://doi.org/10.1146/annurev-immunol-032712-100008https://doi.org/10.1517/14796694.1.1.7https://doi.org/10.1021/acsinfecdis.3c00034https://doi.org/10.1002/cam4.7287https://doi.org/10.1016/j.jddst.2024.105617https://doi.org/10.1016/j.bcp.2023.115741https://doi.org/10.7555/JBR.37.20230022https://doi.org/10.1016/j.canlet.2023.216382https://www.ncbi.nlm.nih.gov/pubmed/37666293https://doi.org/10.1007/s00280-023-04579-8https://www.ncbi.nlm.nih.gov/pubmed/37668680https://doi.org/10.1016/j.ejmech.2024.116440https://www.ncbi.nlm.nih.gov/pubmed/38678825https://doi.org/10.1080/07391102.2022.2099976https://doi.org/10.1007/s11224-015-0578-8https://doi.org/10.1080/07391102.2017.1417912https://doi.org/10.1007/s11224-014-0427-1https://doi.org/10.1021/acsomega.4c01545https://doi.org/10.1016/j.comptc.2023.114188Molecules 2025, 30, 2732 19 of 2120. Sukumaran, S.; Zochedh, A.; Chandran, K.; Sultan, A.B.; Karthiresan, T. Exploring the co-activity of FDA approved drug gemc-itabine and docetaxel for enhanced anti-breast cancer activity: DFT, docking, molecular dynamics simulation and pharmacophorestudies. Int. J. Quantum Chem. 2024, 124, e27359. [CrossRef]21. Ajeel, F.N.; Bardan, K.H.; Kareem, S.H.; Khudhair, A.M. Pd doped carbon nanotubes as a drug carrier for Gemcitabine anticancerdrug: DFT studies. Chem. Phys. Impact 2023, 7, 100298. [CrossRef]22. Venkataramanan, N.S.; Suvitha, A.; Sahara, R.; Kawazoe, Y. Unveiling the gemcitabine drug complexation with cucurbit[n]urils(n = 6–8): A computational analysis. Struct. Chem. 2023, 34, 1869–1882. [CrossRef]23. Roohi, H.; Rouh, M.; Facehi, A. Assessing the performance of Al- and Ga-doped BNNTs for sensing and delivering Cytarabineand Gemcitabine anticancer drugs: A M06-2X study. Mol. Phys. 2023, 121, e2233639. [CrossRef]24. Bibi, S.; Ur-rehman, S.; Khalid, L.; Bhatti, I.A.; Bhatti, H.N.; Iqbal, J.; Bai, F.Q.; Zhang, H.-X. Investigation of the adsorptionproperties of gemcitabine anticancer drug with metal-doped boron nitride fullerenes as a drug-delivery carrier: A DFT study.RSC Adv. 2022, 12, 2873–2887. [CrossRef] [PubMed]25. Adhikary, A.; Kumar, A.; Rayala, R.; Hindi, R.M.; Adhikary, A.; Wnuk, S.F.; Sevilla, M.D. One-Electron Oxidation of Gemcitabineand Analogs: Mechanism of Formation of C3′ and C2′ Sugar Radicals. J. Am. Chem. Soc. 2014, 136, 15646–15653. [CrossRef]26. Malekzadeh, M.; Heshmati, E.; Badalkhani-Khamseh, F.; Nojoumi, S.A. Investigation of conformational structures of gemcitabineand its 2′,2′-difluoro 2′-deoxy derivatives: A computational study. Comput. Theor. Chem. 2020, 1177, 112727. [CrossRef]27. Ali, M.S.; Muthukumaran, J.; Jain, M.; Al-Lohedan, H.A.; Farah, M.A.; Alsowilem, O.I. Experimental and computationalinvestigation on the binding of anticancer drug gemcitabine with bovine serum albumin. J. Biomol. Struct. Dyn. 2022, 40,9144–9157. [CrossRef]28. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Petersson, G.A.;Nakatsuji, H.; et al. Gaussian 16, Rev. C.01; Gaussian, Inc.: Wallingford, CT, USA, 2016.29. Zhang, J.; Dolg, M. Global optimization of clusters of rigid molecules using the artificial bee colony algorithm. Phys. Chem. Chem.Phys. 2016, 18, 3003–3010. [CrossRef]30. Zhang, J.; Dolg, M. ABCluster: The artificial bee colony algorithm for cluster global optimization. Phys. Chem. Chem. Phys. 2015,17, 24173–24181. [CrossRef]31. Walker, M.; Harvey, A.J.A.; Sen, A.; Dessent, C.E.H. Performance of M06, M06-2X, and M06-HF Density Functionals forConformationally Flexible Anionic Clusters: M06 Functionals Perform Better than B3LYP for a Model System with Dispersionand Ionic Hydrogen-Bonding Interactions. J. Phys. Chem. A. 2013, 117, 12590–12600. [CrossRef]32. Hohenstein, E.G.; Chill, S.T.; Sherrill, C.D. Assessment of the Performance of the M05−2X and M06−2X Exchange-CorrelationFunctionals for Noncovalent Interactions in Biomolecules. J. Chem. Theory Comput. 2008, 4, 1996–2000. [CrossRef] [PubMed]33. Thanthiriwatte, K.S.; Hohenstein, E.G.; Burns, L.A.; Sherrill, C.D. Assessment of the performance of DFT and DFT-D methods fordescribing distance dependence of hydrogen-bonded interactions. J. Chem. Theory Comput. 2011, 7, 88–96. [CrossRef] [PubMed]34. Keith, T.A. AIMAll (Version 19.10.12), TK Gristmill Software, Overland Park KS. 2016. Available online: http://aim.tkgristmill.com (accessed on 25 October 2024).35. Venkataramanan, N.S.; Suvitha, A.; Vijayaraghavan, A.; Thamotharan, S. Investigation of inclusion complexation of ac-etaminophen with pillar[5]arene: UV–Vis, NMR and quantum chemical study. J. Mol. Liq. 2017, 241, 782–791. [CrossRef]36. Bulat, F.A.; Toro-Labbé, A.; Brinck, T.; Murray, J.S.; Politzer, P. Quantitative analysis of molecular surfaces: Areas, volumes,electrostatic potentials and average local ionization energies. J. Mol. Model. 2016, 16, 1679–1691. [CrossRef]37. Politzer, P.; Murray, J. The fundamental nature and role of the electrostatic potential in atoms and molecules. Theor. Chem. Acc.2002, 108, 134–142. [CrossRef]38. Venkataramanan, N.S.; Suvitha, A.; Sirajunnisa, A.R.; Sahara, R. Intermolecular hydrogen bond interactions in water clusters ofzwitterionic glycine: DFT, MESP, AIM, RDG, and molecular dynamics analysis. J. Mol. Liq. 2024, 396, 123932.39. Varadwaj, A.; Varadwaj, P.R.; Yamashita, K. Do surfaces of positive electrostatic potential on different halogen derivatives inmolecules attract? like attracting like! J. Comput. Chem. 2018, 39, 343–350. [CrossRef]40. Schrenckenbach, G.; Ziegler, T. Calculation of NMR Shielding Tensors Using Gauge-Including Atomic Orbitals and ModernDensity Functional Theory. J. Phys. Chem. 1995, 99, 606–611. [CrossRef]41. Lu, T.; Chen, F. Multiwfn: A multifunctional wavefunction analyzer. J. Comput. Chem. 2012, 33, 580–592. [CrossRef]42. Chemcraft—Graphical Software for Visualization of Quantum Chemistry Computations. Version 1.8, Build 682. Available online:https://www.chemcraftprog.com (accessed on 10 February 2025).43. Nesse, F. A perspective on the future of quantum chemical software: The example of the ORCA program package. FaradayDiscuss. 2024, 254, 295–314. [CrossRef]44. Sathiyamoorthy, V.N.; Suvitha, A.; Sahara, R.; Kawazoe, Y. Intermolecular interactions in microhydrated ribonucleoside anddeoxyribonucleoside: A computational study. Comput. Theor. Chem. 2021, 1204, 113422. [CrossRef]45. Suvitha, A.; Venkataramanan, N.S.; Sahara, R. The structure, stability, thermochemistry, and bonding in SO3-(H2O)n (n = 1–7)clusters: A computational analysis. Struct. Chem. 2023, 34, 225–237. [CrossRef]https://doi.org/10.1002/qua.27359https://doi.org/10.1016/j.chphi.2023.100298https://doi.org/10.1007/s11224-023-02133-zhttps://doi.org/10.1080/00268976.2023.2233639https://doi.org/10.1039/D1RA09319Chttps://www.ncbi.nlm.nih.gov/pubmed/35425316https://doi.org/10.1021/ja5083156https://doi.org/10.1016/j.comptc.2020.112727https://doi.org/10.1080/07391102.2021.1924270https://doi.org/10.1039/C5CP06313Bhttps://doi.org/10.1039/C5CP04060Dhttps://doi.org/10.1021/jp408166mhttps://doi.org/10.1021/ct800308khttps://www.ncbi.nlm.nih.gov/pubmed/26620472https://doi.org/10.1021/ct100469bhttps://www.ncbi.nlm.nih.gov/pubmed/26606221http://aim.tkgristmill.comhttp://aim.tkgristmill.comhttps://doi.org/10.1016/j.molliq.2017.06.095https://doi.org/10.1007/s00894-010-0692-xhttps://doi.org/10.1007/s00214-002-0363-9https://doi.org/10.1002/jcc.25125https://doi.org/10.1021/j100002a024https://doi.org/10.1002/jcc.22885https://www.chemcraftprog.comhttps://doi.org/10.1039/D4FD00056Khttps://doi.org/10.1016/j.comptc.2021.113422https://doi.org/10.1007/s11224-022-02085-wMolecules 2025, 30, 2732 20 of 2146. Balci, F.M.; Uras-Aytemiz, N. Isomeric forms of 1,4-dioxane in a microsolvation environment. Comput. Theor. Chem. 2024, 1236, 114587.[CrossRef]47. Jurečka, P.; Šponer, J.; Černy, J.; Hobza, P. Benchmark database of accurate (MP2 and CCSD(T) complete basis set limit) interactionenergies of small model complexes, DNA base pairs, and amino acid pairs. Phys. Chem. Chem. Phys. 2006, 8, 1985–1993. [CrossRef]48. Chen, Z.; Shao, H. Exploring the in-situ conversion of intermolecular hydrogen bonds into static electricity through experimentaland theoretical methods. J. Mol. Liq. 2024, 408, 125374. [CrossRef]49. Venkataramanan, N.S.; Suvitha, A.; Sahara, R. Structure, stability, and nature of bonding between high energy water clustersconfined inside cucurbituril: A computational study. Comput. Theor. Chem. 2019, 1148, 44–54. [CrossRef]50. Nakhaei, E.; Nowroozi, A.; Ravari, F. The hydrogen-bonded complexes of the 5-fluorouracil with the DNA purine bases: Acomprehensive quantum chemical study. Struct. Chem. 2018, 29, 69–80. [CrossRef]51. Meredith, N.Y.; Borsley, S.; Smolyar, I.V.; Nichol, G.S.; Baker, C.M.; Ling, K.B.; Cockroft, S.L. Dissecting Solvent Effects onHydrogen Bonding. Angew. Chem. 2022, 61, e202206604. [CrossRef]52. Liu, C.; Gong, L.; Gong, B. Theoretical investigation of the electronic and second-order non-linear optical properties of [n] helicenederivatives. Comput. Theor. Chem. 2024, 1238, 114699. [CrossRef]53. Anan, J.; Adu Fosu, E.; Obuah, C.; Kojo Ainooson, M.; Aniagyei, A.; Hamenu, L.; Oppong, A.; Muller, A. A DFT and TD-DFTstudies of the photosensitizing capabilities of thiophene-based dyes. Comput. Theor. Chem. 2024, 1237, 114633. [CrossRef]54. Pal, R.; Chattaraj, P.K. Electrophilicity index revisited. J. Comput. Chem. 2023, 44, 278–297. [CrossRef] [PubMed]55. Edet, H.O.; Louis, H.; Benjamin, I.; Gideon, M.; Unimuke, T.O.; Adalikwu, S.A.; Nwagu, A.D.; Adeyinka, A.S. Hydrogen storagecapacity of C12X12 (X = N, P, and Si). Chem. Phys. Impact. 2022, 5, 100107. [CrossRef]56. Tegegn, D.F.; Zewude Belachew, H.; Salau, A.O. Investigation of substituent effects on the electronic structure and antiviralactivity of favipiravir derivatives for COVID-19 treatment using DFT and molecular docking. Sci. Rep. 2024, 14, 8146. [CrossRef]57. Rajimon, K.J.; Sreelakshmi, N.; Rajendran Nair, D.S.; Fazulunnisa Begum, N.; Thomas, R. An in-depth study of the synthesis,electronic framework, and pharmacological profiling of 1-(anthracen-9-yl)-N-(4-nitrophenyl) methanimine: In vitro and in silicoinvestigations on molecular docking, dynamics simulation, BSA/DNA binding and toxicity studies. Chem. Phys. Impact 2024, 8, 100462.58. Suvitha, A.; Venkataramanan, N.S.; Sahara, R.; Kawazoe, Y. A theoretical exploration of the intermolecular interactions betweenresveratrol and water: A DFT and AIM analysis. J. Mol. Model. 2019, 25, 56. [CrossRef]59. Saleh, G.; Gatti, C.; Lo Presti, L. Non-covalent interaction via the reduced density gradient: Independent atom model vsexperimental multipolar electron densities. Comput. Theor. Chem. 2012, 998, 148–163. [CrossRef]60. Baranac-Stojanović, M.; Aleksić, J.; Stojanović, M. Theoretical investigation of tautomerism of 2-and 4-pyridones: Origin,substituent and solvent effects. Org. Biomol. Chem. 2024, 22, 144–158. [CrossRef]61. Jiang, X.; Zhang, H.; Wu, W.; Mo, Y. Critical Check for the Role of Resonance in Intramolecular Hydrogen Bonding. Chem.-Eur.J.2017, 23, 16885–16891. [CrossRef]62. Karuppasamy, J.; Zochedh, A.; Al-Asbahi, B.A.; Kumar, Y.A.; Shanmuganarayanan, A.; Sultan, A.B. Combined experimental andquantum computational studies on 5-Fluorouracil salicylic acid cocrystal: Assessment on anti-breast cancer potential throughdocking simulation. J. Mol. Struct. 2024, 1311, 138406. [CrossRef]63. Al-Wahaibi, L.H.; Kumar, N.S.; El-Emam, A.A.; Venkataramanan, N.S.; Ghabbour, H.A.; Tamimi, A.-M.S.; Percino, J.; Thamotharan,S. Invesigation of potential anti-malarial lead candidate 2-(4-fluorobenzylthio)-5-(5-bormothiophene-2-yl)-1,3,4-oxadiazole:Insights from crystal structure, DFT, QTAIM and hybrid QM/MM binding energy analysis. J. Mol. Struct. 2019, 1175, 230–240.[CrossRef]64. Varadwaj, P.R. Halogen Bond via an Electrophilic π-Hole on Halogen in Molecules: Does It Exist? Int. J. Mol. Sci. 2024, 25, 4587.[CrossRef] [PubMed]65. Murray, J.S.; Politzer, P. The electrostatic potential: An overview. WIREs Comput. Mol. Sci. 2011, 1, 153–163. [CrossRef]66. Roos, G.; Murray, J.S. Probing intramolecular interactions using molecular electrostatic potentials: Changing electron densitycontours to unveil both attractive and repulsive interactions. Phys. Chem. Chem. Phys. 2024, 26, 7592–7601. [CrossRef]67. Thaynara Guimarães, M.; Nicolas Nascimento, C.; Murielly Fernanda Ribeiro, B.; Anna Karla dos Santos, P.; Grasiele Soares, C.;Douglas Henrique, P. Interactions between DNA and the acridine intercalator: A computational study. Comput. Biol. Chem. 2024,109, 108029.68. Wolstenholme, D.J.; Stanley Cameron, T. Comparative study of weak interactions in molecular crystals: H− H bonds vs hydrogenbonds. J. Phys. Chem. A. 2006, 110, 8970–8978. [CrossRef]69. Shishkin, O.V.; Palamarchuk, G.V.; Gorb, L.; Leszczynski, J. Intramolecular Hydrogen Bonds in Canonical 2 ‘-Deoxyribonucleotides: An Atoms in Molecules Study. J. Phys. Chem. B. 2006, 110, 4413–4422. [CrossRef]70. Gassoumi, B.; Mehri, A.; Hammami, H.; Castro, M.E.; Karayel, A.; Özkιnal, S.; Melendez, F.J.; Nouar, L.; Madi, F.; Ghalla, H.; et al.Spectroscopic characterization, host-guest charge transfer, Hirshfeld surfaces, AIM-RDG and ELF study of adsorption and chemicalsensing of heavy metals with new derivative of Calix[4]quinone: A DFT-D3 computation. Mater. Chem. Phys. 2022, 278, 125555.[CrossRef]https://doi.org/10.1016/j.comptc.2024.114587https://doi.org/10.1039/B600027Dhttps://doi.org/10.1016/j.molliq.2024.125374https://doi.org/10.1016/j.comptc.2018.12.015https://doi.org/10.1007/s11224-017-1001-4https://doi.org/10.1002/anie.202206604https://doi.org/10.1016/j.comptc.2024.114699https://doi.org/10.1016/j.comptc.2024.114633https://doi.org/10.1002/jcc.26886https://www.ncbi.nlm.nih.gov/pubmed/35546516https://doi.org/10.1016/j.chphi.2022.100107https://doi.org/10.1038/s41598-024-68712-0https://doi.org/10.1007/s00894-019-3941-7https://doi.org/10.1016/j.comptc.2012.07.014https://doi.org/10.1039/D3OB01588Bhttps://doi.org/10.1002/chem.201703952https://doi.org/10.1016/j.molstruc.2024.138406https://doi.org/10.1016/j.molstruc.2018.07.102https://doi.org/10.3390/ijms25094587https://www.ncbi.nlm.nih.gov/pubmed/38731806https://doi.org/10.1002/wcms.19https://doi.org/10.1039/D3CP06005Ehttps://doi.org/10.1021/jp061205ihttps://doi.org/10.1021/jp056902+https://doi.org/10.1016/j.matchemphys.2021.125555Molecules 2025, 30, 2732 21 of 2171. Bouchemela, H.; Madi, F.; Nouar, L. DFT investigation of host–guest interactions between α-Terpineol and β-cyclodextrin. J. Incl.Phenom. Macrocyclic Chem. 2022, 95, 247–258. [CrossRef]72. Venkataramanan, N.S.; Suvitha, A. Encapsulation of sulfur, oxygen, and nitrogen mustards by cucurbiturils: A DFT study. J. Incl.Phenom. Macrocyclic Chem. 2015, 83, 387–400. [CrossRef]73. Wuttke, A.; Matta, R.A. Visualizing dispersion interactions through the use of local orbital spaces. J. Comput. Chem. 2017, 38,15–23. [CrossRef]74. Bistoni, G. Finding chemical concepts in the Hilbert space: Coupled cluster analyses of noncovalent interactions. Wiley Interdiscip.Rev. Comput. Mol. Sci. 2020, 10, e1442. [CrossRef]75. Altun, A.; Izsák, R.; Bistoni, G. Local energy decomposition of coupled-cluster interaction energies: Interpretation, benchmarks,and comparison with symmetry-adapted perturbation theory. Int. J. Quantum Chem. 2021, 121, e26339. [CrossRef]76. Venkataramanan, N.S. Evaluating the drug delivery and sensing performance of XB23N24,(X = B, Al, Ga) nanocages forgemcitabine anticancer drug. Comput. Theor. Chem. 2025, 1248, 115222. [CrossRef]Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individualauthor(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury topeople or property resulting from any ideas, methods, instructions or products referred to in the content.https://doi.org/10.1007/s10847-019-00940-8https://doi.org/10.1007/s10847-015-0575-yhttps://doi.org/10.1002/jcc.24508https://doi.org/10.1002/wcms.1442https://doi.org/10.1002/qua.26339https://doi.org/10.1016/j.comptc.2025.115222 Introduction  Computational Methods  Results and Discussion  Structures and Geometric Parameters of the Gemcitabine–Nucleobase Complexes  Energetics  Frontier Molecular Orbitals and Conceptual DFT  Nature of Interactions  -Electron Delocalization Indicators and Molecular Electrostatic Potential  Atom in Molecule Analysis and Non-Covalent Interaction Analysis  Dispersion Interaction Density (DID) and Local Energy Decomposition (LED) Analyses  Role of Intramolecular Hydrogen Bond of Gemcitabine in Stabilizing the Complex  Conclusions  References