The chemical reactivity of a compound can be interpreted as the resistance or ease with which it attracts or gives away electrons under the action of an external potential v(r). In this sense, there are global parameters that allow us to characterize this electron transfer from a theoretical point of view, such as the electronic chemical potential (μ)  , the molecular hardness (η)  and the electrophilicity index (ω)  . Since local descriptors of reactivity, such as the Fukui function (FF)    and local softness (s), are properties that depend on the position (r), they may explain selectivity in certain parts of a molecule.
Felodipine (FD), (methyl-4-(2,3-dichlorophenyl)-1,4-dihydro-2,6-dimethyl- 3,5-pyridinedicarboxylate), shown in Figure 1, is a calcium channel blocker and an effective and widely prescribed medication in the treatment of hypertension and angina  .
Four crystalline forms of FD have been reported  , all shaping an orthorhombic system. Form I, used in some commercialized products, is the most stable, where its crystalline structure, determined in 1986, has a space group P2(1)/c, and a unit cell of a = 12.086 Å, b = 12.077 Å, c = 13.425 Å and β = 116.13˚ . Many examples have been described in which the compound polymorphs show variation in reactivity    . Furthermore, in the same polymorph, the reactivity can be anisotropic with respect to the crystalline faces where remaining chemicals are different, whereby the reaction rates along particular crystallographic directions can vary significantly. Paul and Curtin were pioneers in conducting solid state reaction studies, specifically the acid-base character of gas-solid reactions, such as with ammonia gas and p-chlorobenzoic anhydride crystals  . It was observed that, often, certain crystal faces are attacked by gases preferentially and the reactions propagate along specific directions. However, there are still few reported studies on the solid state reactivity of organic crystals at the electronic level. Currently, the use of computational methods as an indispensable tool for the understanding of physical and chemical properties has led to important advances in the molecular sciences, and in the last decades, in particular, there has been an increased interest in the concepts of Density Functional Theory (DFT)   - . Important examples are the recent study by Luty et al. based on Fukui’s nuclear functions using DFT to investigate the explosive mechanism of RDX (hexahydro-1,3,5-trinitro-1,3,5-triazine)  , and the study conducted by Shaoxin Feng and Tonglei Li, determining that the Fukui nuclear function can be used to characterize the difference in the chemical reactivity of two polymorphs of flufenamic acid . Since a chemical reaction is driven by the change in energy of the system and it is accompanied by electron transfer and atomic displacement, the Fukui nuclear function, a local function for describing the sensitivity
Figure 1. Molecular structure of felodipine (FD), showing the atom-labelling scheme.
of the system to a simultaneous disturbance in the number of electrons N, may be useful to characterize the reactivity of crystals with respect to crystal packing .
The compounds of the group of 1,4-DHP drugs undergo hepatic metabolization (by which the oxidation of 1,4-dihydropyridines in pyridines occurs) catalyzed by the CYP3A enzyme of cytochrome P450  . The metabolization reaction takes place in a two-step mechanism: in the first step the hydrogen is abstracted from the dihydropyridine, where as in the second a hydroxyl group is added to the substrate.
There are many computational approaches that attempt to explain the process of metabolism. Some of them use statistical models in molecular descriptors   , other works use ligands to predict the preferred sites of hydrogen abstraction by quantum chemistry  . Another, simpler, yet computationally expensive, way is to chemically model the complete process of hydroxylation catalyzed by cytochrome   . The hydroxylation reaction has been modeled by coupled quantum/molecular mechanics (QM/MM) . These calculations strongly support the “two-state reactivity” (TSR) model, in which, the active site of the cytochrome―where the oxidation takes place―is a ferro-oxyl species called “Compound I” (CpdI) from the current knowledge . Cpd I can be seen as an “electrophilic oxidant” . Therefore, the Fukui function ( ) should help to pinpoint sites in a molecule that prefer to be attacked by Cpd I. On the contrary, the Fukui function , evaluated for Cpd I prefer to be attacked by dihydropyridine .
Felodipine shows a particularly rich variety of metabolites  . Most involve the oxidation of 1,4-dihydropyridine to pyridine, sometimes together with the loss of one of the esters, which is well reflected in . For FD and other compounds, Singh et al.  estimated the abstraction energies of the H in a semi-empirical way and its correlation with quite successful hydroxylation sites. Several experiments were carried out to study the structure and properties of felodipine. However, few investigations on the crystal structure faces are available to date.
In this report we investigate how crystal morphology can relate to properties of FD. The novelty of this work lies in studying the surfaces (100), (010) and (001), especially the first one, with supercells to model the periodic systems, based on DFT with the Generalized Gradient Approximation (GGA) and a dispersion correction as formulated by Tkatchenko and Scheffler (TS). By using the FF the most reactive sites of the surfaces to electrophilic species are described.
We will lay the basic definitions of the global descriptors according to DFT. The energy E is expressed in terms of the number of electrons N and an external potential v(r), so that , where is the electronic density    . The derivatives of with respect to N and v(r) lead to a set of global and local properties that allows us to quantify the concept of reactivity and site selectivity. As is well known in the literature, global descriptors are properties that explain the stability of a molecule, such as, the electronic chemical potential (µ), molecular hardness (η) and electrophilicity index (ω)  - . The variation of energy with respect to the external potential gives local information, i.e., it depends on the position (r), and therefore it is defined as an index of selectivity. Local descriptors such as the Fukui function (FF)    and the local softness (s)  are properties that explain the selectivity of a region in a molecule. The FF is defined as:
Equation (1) presents a discontinuity problem in atoms and molecules when combined with the finite difference approximation, resulting in expressions of the FF   :
when a molecule accepts electrons, they tend to move around places where is large because at these locations the molecule is most likely to stabilize additional electrons, and therefore, it is susceptible to nucleophilic attack at such sites. Likewise, a molecule is susceptible to electrophilic attack at sites where is large, since these are the regions where electron removal least destabilizes the molecule. In chemical density functional theory, the FF are the key of region selectivity indicators for electron-transfer controlled reactions. The function quantification is possible through a system of condensation on an atomic region, where FF can be written by using population analysis techniques . This leads to the following equations:
where denotes the electronic population of atom k of the reference system, more correctly called .
The Fukui functions favourably determinate the reactive sites for most chemical systems. However, the values of the FF rely upon the scheme chosen to calculate the charges and the accuracy of the population analysis used. In this study the Hirshfeld method was used  , which has been proved to be accurate for organic systems .
The evaluation of the surface energies of the crystalline faces can be useful to compare how the resistance of different surfaces affects the kinetics of the reaction. As often observed, a solid state reaction can proceed in a specific crystallographic direction. Therefore, the study of mechanical properties on different faces can give some insight into the solid state reaction. The surface energy can be calculated as  :
where Eslab and Ebulk are total energies of the surface and crystal in bulk, respectively, and n is the thickness (or layers of unit cells) of the surface, and S is the surface area. In the present study, the surface energies, the electronic structures, and Fukui nuclear functions of three surfaces of the crystal form I of FD were calculated.
The crystal structure of FD, form I, was obtained from the Cambridge Structural Database (CSD) (ref code: DONTIJ). FD crystallizes in an orthorhombic lattice with space group P2(1)/c, and cell parameters a = 12.086 Å, b = 12.077 Å, c = 13.425 Å and β = 116.13˚, packed with four molecules .
A periodic solid state program was used, with DFT-D (dispersion correction)      at the GGA level developed by Perdew-Burke-Ernzerhof (PBE)  , with plane waves as a basis set with a cutting energy of 380 eV and a tolerance for the SCF cycle of 10−6 eV/atom. Vanderbilt pseudopotentials  were used to model ion-electron interactions with: H: 1s1, C: 2s22p2, N: 2s22p3, O: 2s22p4, Cl: 3s23p5. Surface (100) (a = 12.077 Å, b = 13.425 Å, c = 32.736 Å and α = β = γ = 90˚), Surface (010) (a = 13.425 Å, b = 12.086 Å, c = 43.082 Å and α = β = 90˚, γ = 116.13˚) and Surface (001) (a = 12.086 Å, b = 12.077 Å, c = 42.005 Å and α = β = γ = 90˚) of form I were modeled (Figure 2) and their Fukui nuclear functions FF were analyzed as they may be linked to their chemical reactivity.
3. Results and Discussion
The electronic structure of the three surfaces of Form I of FD was studied. Also, the FF and the surface energy for each were calculated. The bulk crystal structure was optimized with the same methods that were used to calculate the FF. The network parameters were set during the optimization. Surface models of (100), (010) and (001) faces, of two unitary cells of thickness were thus constructed. On the (100) and (001) faces the rings of the pyridine are exposed on the surface unlike the face (010) where the benzyl ring together with the methyl ether group are more exposed.
We focus our attention on the (100) surface model which turned out to be the most stable, in which there are intermolecular hydrogen bonds formed by N-H-O located on the normal plane of the surface (see Figure 3).
Once the optimization of the geometry of the surface model (100) was carried out, it was observed that the hydrogen bonds lying on the surface tend to move downwards, along the z-axis. The above can be seen in Table 1, where the values corresponding to the midpoint of the hydrogen bonds are given as well as the displacement after optimizing such surface of FD form I. The hydrogen bonds corresponding to those in the layer in between have a displacement along the
Figure 2. Surfaces (100), (010) and (001) of form I of the felodipine. The light blue dotted line shows the hydrogen bonds of the surface models.
Figure 3. Intermolecular hydrogen bonds of the surface (100) of FD.
Table 1. Midpoint coordinates of the hydrogen bonds and their displacement after optimizing surface (100) of felodipine form I.
z-axis lower than those on the surface. From the above it can be observed that surface molecules, having no interaction with other molecules in the upper part, tend to flatten the surface. This displacement of the hydrogen bonds causes an activation of the carbon atoms C2 and C4, which favors the reactivity of these atoms as shown by Fukui functions (Table 2), making them more susceptible to nucleophilic substitutions.
Figure 4 shows the boundary orbitals HOMO (Highest Occupied Molecular Orbital) and LUMO (Lowest Unoccupied Molecular Orbital), from which a high density of unoccupied states in the three surfaces can be seen, located mainly in the carbon atoms 1 and 5 of the 1,4-pyridine ring, which gives us an idea of the places most likely to interact with electron donor atoms.
The FF of each atom was obtained from the calculations of the neutral and anionic forms of the surface models. The molecular structure of the anionic form remained the same as its neutral counterpart. The quantitative results of the FF for faces (100), (010) and (001) are listed in Table 2. The electrophilic attack of the crystalline faces is illustrated in Figure 5. Fukui nuclear functions of
Table 2. Calculated Nuclear Fukui Functions of the (100), (010) and (001) slab models of Form I of felodipine with DFT-D.
Figure 4. HOMO and LUMO orbitals of surfaces (100), (010) and (001) of FD form I.
Figure 5. Illustration of Fukui nuclear function ( blue and red) for (001), (010) and (100) surfaces of felodipine Form I.
the atoms N1, C1 and C4 of the pyridinic ring and the atoms O1 and O3 of the carbonyls groups in the three surfaces are larger than those of other atoms. These results agree with those reported by Michael E. Beck .
In Figure 5 it can be seen that the (100) surface has a distribution of molecules favoring electrophilically reactive sites, thus making this surface susceptible to greater reactivity to agents such as Cpd I, which, as already mentioned, is the active site of the cytochrome, where the oxidation of dihydropyridine to pyridine is carried out.
To elucidate the influence of mechanical resistance on chemical reactivity, the surface energies of the surface models of Form I were calculated by the DFT-D with the functional GGA, obtaining the following Esurf: −0.2303 eV/Å2 for (100), 0.0222 eV/Å2 for (010) and 0.0302 eV/Å2 for (001) surfaces. The latter energies remain significantly above the one for (100) surface, thus indicating a closer bond between the molecules on the faces (010) and (001). As the reaction begins on the surfaces, it is expected that their propagation and penetration in the bulk are limited by forces of intermolecular nature. The surface energy characterizes the intermolecular interactions within a crystallographic plane, specifying the mechanical resistance. Therefore, the lower surface energy of (100) may facilitate a faster reaction rate than the other two surfaces
The reaction capacity of the three analyzed surfaces of form I of felodipine was investigated by studying their electronic structures, nuclear Fukui functions and their surface energies. The present findings show that Fukui nuclear functions constitute a useful tool in the analysis of surface reactivity for a crystal such as FD. In addition, due to the highly heterogeneous nature of a molecular crystal reaction, the surface intermolecular forces ought to be taken into account to elucidate chemical reactions occurring on this type of crystals.
These results can provide information on experimental work in surface catalysis as based on theoretical knowledge of local reactivity of the compound here analyzed, thereby saving efforts when selecting the best sites a priori. Our findings may also be useful in some pharmaceutical applications of felodipine.
One of us (CTC) gratefully acknowledges a PhD fellowship granted by CONACyT (Mexico). JFRS and AFR are grateful to Sistema Nacional de Investigadores (SNI, Mexico) for provision of funds.
 Parr, R.G. and Pearson, R.G. (1983) Absolute Hardness: Companion Parameter to Absolute Electronegativity. Journal of the American Chemical Society, 105, 7512-7516. https://doi.org/10.1021/ja00364a005
 Liu, S.B. and Chattaraj, P.K. (2009) Electrophilicity. In: Chattaraj, P.K., Ed., Chemical Reactivity Theory: A Density Functional View, Taylor and Francis, Boca Raton, 179. https://doi.org/10.1201/9781420065442.ch13
 Ayers, P.W. and Levy, M. (2000) Perspective on “Density Functional Approach to the Frontier-Electron Theory of Chemical Reactivity”. Theoretical Chemistry Accounts, 103, 353-360. https://doi.org/10.1007/s002149900093
 Perdew, J.P., Parr, R.G., Levy, M. and Balduz, J.L. (1982) Density-Functional Theory for Fractional Particle Number: Derivative Discontinuities of the Energy. Physical Review Letters, 49, 1691-1694. https://doi.org/10.1103/PhysRevLett.49.1691
 Ayers, P.W. (2008) The Dependence on and Continuity of the Energy and Other Molecular Properties with Respect to the Number of Electrons. Journal of Mathematical Chemistry, 43, 285-303. https://doi.org/10.1007/s10910-006-9195-5
 Surov, A.O., Solanko, K.A., Bond, A.D., Perlovich, G.L. and Bauer-Brand, A. (2012) Crystallization and Polymorphism of Felodipine. Crystal Growth & Design, 12, 4022-4030. https://doi.org/10.1021/cg300501u
 Fossheim, R. (1986) Crystal Structure of the Dihydropyridine Calcium Antagonist Felodipine. Dihydropyridine Binding Prerequisites Assessed from Crystallographic Data. Medicinal Chemistry, 29, 305-307. https://doi.org/10.1021/jm00152a023
 Byrn, S.R., Sutton, P.A., Tobias, B., Frye, J. and Main, P. (1988) Crystal Structure, Solid-State NMR Spectra, and Oxygen Reactivity of Five Crystal Forms of Prednisolone Tert-Butylacetate. Journal of the American Chemical Society, 101, 1609-1614. https://doi.org/10.1021/ja00213a039
 Chen, X.M., Morris, K.R., Griesser, U.J., Byrn, S.R. and Stowell, J.G. (2002) Reactivity Differences of Indomethacin Solid Forms with Ammonia Gas. Journal of the American Chemical Society, 124, 15012-15019. https://doi.org/10.1021/ja017662o
 Chen, X.M., Li, T.L., Morris, K.R. and Byrn, S.R. (2002) Crystal Packing and Chemical Reactivity of Two Polymorphs of Flufenamic Acid with Ammonia. Molecular Crystals and Liquid Crystals, 381, 121-131. https://doi.org/10.1080/713738743
 Shaoxin, F. and Li, T.L. (2005) Understanding Solid-State Reactions of Organic Crystals with Density Functional Theory-Based Concepts. Journal of Physical Chemistry A, 109, 7258-7263. https://doi.org/10.1021/jp0519666
 Luty, T., Ordon, P. and Eckhardt, C.J.J. (2002) A Model for Mechanochemical Transformations: Applications to Molecular Hardness, Instabilities, and Shock Initiation of Reaction. Chemical Physics, 117, 1775. https://doi.org/10.1063/1.1485968
 Ruiz, V.G., Liu, W., Zojer, E., Scheffler, M. and Tkatchenko, A. (2012) Density-Functional Theory with Screened van der Waals Interactions for the Modeling of Hybrid Inorganic-Organic Systems. Physical Review Letters, 108, Article ID: 146103. https://doi.org/10.1103/PhysRevLett.108.146103
 Klimeš, J. and Michaelides, A. (2012) Perspective: Advances and Challenges in Treating van der Waals Dispersion Forces in Density Functional Theory. The Journal of Chemical Physics, 137, Article ID: 120901. https://doi.org/10.1063/1.4754130
 Katoh, M., Nakajima, M., Shimada, N., Yamazaki, H. and Yokoi, T. (2000) Inhibition of Human Cytochrome P450 Enzymes by 1,4-Dihydropyridine Calcium Antagonists: Prediction of in Vivo Drug-Drug Interactions. European Journal of Clinical Pharmacology, 55, 843-852. https://doi.org/10.1007/s002280050706
 Rodríguez Arcas, M.J., García-Jiménez, E., Martínez-Martínez, F. and Conesa-Zamora, P. (2011) Role of CYP450 in Pharmacokinetics and Pharmacogenetics of Antihypertensive Drugs. Farmacia Hospitalaria, 35, 84-92. https://doi.org/10.1016/j.farma.2010.05.006
 Langowski, J. and Long, A. (2002) Anthony Long, Computer Systems for the Prediction of Xenobiotic Metabolism. Advanced Drug Delivery Reviews, 54, 407-415. https://doi.org/10.1016/S0169-409X(02)00011-X
 Higgins, L., Korzekwa, K.R., Rao, S., Shou, M. and Jones, J.P. (2001) An Assessment of the Reaction Energetics for Cytochrome P450-Mediated Reactions. Archives of Biochemistry and Biophysics, 385, 220-230. https://doi.org/10.1006/abbi.2000.2147
 Singh, S.B., Shen, L.Q., Walker, M.J. and Sheridan, R.P.J. (2003) A Model for Predicting Likely Sites of CYP3A4-Mediated Metabolism on Drug-Like Molecules. Journal of Medicinal Chemistry, 46, 1330-1336. https://doi.org/10.1021/jm020400s
 Shaik, S., Cohen, S., de Visser, S.P., Sharma, P.K., Kumar, D., Kozuch, S., Ogliaro, F. and Danovich, D. (2004) The “Rebound Controversy”: An Overview and Theoretical Modeling of the Rebound Step in C-H Hydroxylation by Cytochrome P450. European Journal of Inorganic Chemistry, 2004, 207-226. https://doi.org/10.1002/ejic.200300448
 Schöneboom, J.C., Lin, H., Reuter, N., Thiel, W., Cohen, S., Ogliaro, F. and Shaik, S. (2002) The Elusive Oxidant Species of Cytochrome P450 Enzymes: Characterization by Combined Quantum Mechanical/Molecular Mechanical (QM/MM) Calculations. Journal of the American Chemical Society, 124, 8142-8151. https://doi.org/10.1021/ja026279w
 Schöneboom, J.C., Cohen, S., Lin, H., Shaik, S. and Thiel, W. (2004) Quantum Mechanical/Molecular Mechanical Investigation of the Mechanism of C-H Hydroxylation of Camphor by Cytochrome P450cam: Theory Supports a Two-State Rebound Mechanism. Journal of the American Chemical Society, 126, 4017-4034. https://doi.org/10.1021/ja039847w
 Groves, J.T. and Watanabe, Y. (1988) Reactive Iron Porphyrin Derivatives Related to the Catalytic Cycles of Cytochrome P-450 and Peroxidase. Studies of the Mechanism of Oxygen Activation. Journal of the American Chemical Society, 110, 8443-8452. https://doi.org/10.1021/ja00233a021
 Wang, S.X., Sutfin, T.A., Bäärnhielm, C. and Regårdh, C.G. (1989) Contribution of the Intestine to the First-Pass Metabolism of Felodipine in the Rat. Journal of Pharmacology and Experimental Therapeutics, 250, 632-636.
 De Proft, F. and Geerlings, P. (1997) Calculation of Ionization Energies, Electron Affinities, Electronegativities, and Hardnesses Using Density Functional Methods. The Journal of Chemical Physics, 106, 3270-3279. https://doi.org/10.1063/1.473796
 Ayers, R.G. and Parr, P. (2000) Variational Principles for Describing Chemical Reactions: The Fukui Function and Chemical Hardness Revisited. Journal of the American Chemical Society, 122, 2010-2018. https://doi.org/10.1021/ja9924039
 Yang, R.G. and Parr, W. (1985) Hardness, Softness, and the Fukui Function in the Electronic Theory of Metals and Catalysis. Proceedings of the National Academy of Sciences, 82, 6723-6726. https://doi.org/10.1073/pnas.82.20.6723
 Contreras, R.R., Fuentealba, P., Galván, M. and Pérez, P. (1999) A Direct Evaluation of Regional Fukui Functions in Molecules. Chemical Physics Letters, 304, 405-413. https://doi.org/10.1016/S0009-2614(99)00325-5
 Padmanabhan, J., Parthasarathi, R., Sarkar, U., Subramanian, V. and Chattaraj, P.K. (2004) Effect of Solvation on the Condensed Fukui Function and the Generalized Philicity Index. Chemical Physics Letters, 383, 122-128. https://doi.org/10.1016/j.cplett.2003.11.013
 Perdew, J.P. and Zunger, A. (1981) Self-Interaction Correction to Density-Functional Approximations for Many-Electron Systems. Physical Review B, 23, 5048-5079. https://doi.org/10.1103/PhysRevB.23.5048
 Grimme, S. (2006) Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. Journal of Computational Chemistry, 27, 1787-1799. https://doi.org/10.1002/jcc.20495
 Jurecka, P., Cerny, J., Hobza, P. and Salahub, D.R. (2007) Density Functional Theory Augmented with an Empirical Dispersion Term. Interaction Energies and Geometries of 80 Noncovalent Complexes Compared with ab Initio Quantum Mechanics Calculations. Journal of Computational Chemistry, 28, 555-569. https://doi.org/10.1002/jcc.20570
 Ortmann, F., Bechstedt, F. and Schmidt, W.G. (2006) Semiempirical van der Waals Correction to the Density Functional Description of Solids and Molecular Structures. Physical Review B, 73, Article ID: 205101. https://doi.org/10.1103/PhysRevB.73.205101