Introduction

Classical molecular mechanics (MM) and molecular dynamics (MD) simulations have become powerful tools for studying large molecular systems. They can deal with systems consisting of hundreds of thousands of atoms, and can simulate their dynamic behavior over a time scale of nanoseconds. They are especially important when studying biomolecules and condensed phases, as the results of such simulations complement experimental results and provide insight into microscopic properties that cannot be probed experimentally [1, 2]. They are capable of effectively treating these large-scale systems through the utilization of simplified functions that express the system’s potential energy U(R) in terms of the coordinates of the atoms, \( {\bf{R}} = \left( {{{{\bf{R}}}_{1}},\:{{{\bf{R}}}_{2}}, \ldots ,\:{{{\bf{R}}}_{N}}} \right) \), where N is the total number of atoms. In the most popular additive force fields [3–9], electrostatic contributions to U(R) are described using the Coulomb potential between static atomic partial charges:

$$ {{U}_{{{\rm{elst}}}}}({\bf{R}}) = \sum\limits_{i} {\sum\limits_{{j \ne i}} {\frac{{{{q}_{i}}{{q}_{j}}}}{{\varepsilon \cdot {{R}_{{ij}}}}}} } , $$
(1)

where q i and q j are the charges on atoms i and j, respectively, ε is the dielectric constant, and \( {{{\bf{R}}}_{{ij}}} = \left| {{{{\bf{R}}}_{i}} - {{{\bf{R}}}_{j}}} \right| \) is the distance between atoms i and j.

In additive force fields, many-body effects—particularly polarization—are neglected (i.e., electronic degrees of freedom are ignored, and the charge distribution is assumed to remain constant during the simulation). To implicitly account for the polarization that occurs in condensed phases, partial atomic charges determined for the gas phase are systematically overestimated by 10–20 % [10, 11]. This approach takes into account the potential energy contributions neglected in U(R) using an average. There are, however, a variety of systems for which such an approximation is insufficient. In general, these systems contain inhomogeneities where polarizable molecules experience changes in the dielectric environment, for example in interfacial regions: liquid–vapor, liquid–liquid, or protein–solvent surfaces, or solutions of ions. Polarization also plays an important role when various thermodynamic processes have to be considered.

Currently models that explicitly include polarization in MD simulations include polarizable dipoles [12–14], the classical Drude oscillator [15, 16], and the fluctuating charge (FQ) model [17–23]. The first two models allow for the occurrence of induced dipole moments at various sites in the system (atoms, bonds, lone pairs). In the FQ model, the values of atomic partial charges are allowed to vary. This is based on the electronegativity equalization (EE) principle [24–26]. The Drude and FQ models have the advantage of preserving the Coulombic potential for electrostatic interactions without introducing an additional dipole–dipole term. In MD simulations, all three models can be solved iteratively or with extended Lagrangian mechanics [27]. The FQ model can also be solved using matrix inversion.

Charge sensitivity analysis (CSA) was first formulated in the 1990s in the framework of density functional theory (DFT) [28, 29]. Although it was originally designed as a tool to complement ab initio calculations, it can—after suitable parameterization—be used as a semi-empirical method. It provides information about the charge distribution in a molecular system and about the system’s reactivity (Fukui function and softness/hardness data). It is based on the second-order Taylor expansion of the system’s energy with respect to population variables and on the EE principle. From this point of view, the FQ model is similar to CSA. However, this is only one aspect of CSA. CSA undoubtedly provides a wider perspective. The reactivity information built into the CSA formalism can be used in qualitative structure–activity relationship (QSAR) and qualitative structure–property relationship (QSPR) models. The huge potential of conceptual DFT in such models has been demonstrated in the literature [30–36]. In our previous paper [37], we parameterized CSA for the AMBER ff99 [38, 39] force field on the basis of static calculations performed in the global equilibrium state. In this paper, we present the utilization of CSA for other types of equilibria. In MD simulations, it is typical to employ a reactant perspective—each reactant preserves its own identity, manifested as a unique set of bonding parameters. In this paper, in addition to global equilibrium, we will use the approach of imposing a partial equilibrium on the system (i.e., the reactants are mutually polarized but there is no charge transfer between them). Moreover, we will explore the dynamic behavior of partial charges and other CSA parameters along a simulated trajectory of a molecular system.

The paper is organized as follows. We first present a short survey of CSA, before introducing a model system and providing computational details. We then describe the obtained results, and finally present our conclusions and future prospects in this field.

Charge sensitivity analysis

The basis of CSA is the Taylor expansion of the energy of a molecular system M with respect to atomic charges, truncated at the quadratic term:

$$ \begin{array}{*{20}c} \begin{array}{*{20}c} {d{E_{\text{M}}}={d^1}{E_{\text{M}}}+{d^2}{E_{\text{M}}}=\sum\limits_i {\left( {\frac{{\partial {E_{\text{M}}}}}{{\partial {q_i}}}} \right)d{q_i}} +\tfrac{1}{2}\sum\limits_i {\sum\limits_j {\left( {\frac{{{\partial^2}{E_{\text{M}}}}}{{\partial {q_i}\partial {q_j}}}} \right)d{q_i}d{q_j}} } } \\ {=\boldsymbol{\upchi}_0^{\dag}d\mathbf{q}+\tfrac{1}{2}d{{\mathbf{q}}^{\dag}}{{\boldsymbol{\upeta}}_0}d\mathbf{q}} \\ \end{array}, \hfill \\ \quad {\,}{\,}\quad \hfill \\\end{array}. $$
(2)

Here, q is the vector grouping atomic partial charges, \( {{\mathbf{q}}^{\dag}}=\left( {{q_1},\;{q_2},\ldots,\;{q_N}} \right) \), \( \boldsymbol{\upchi}_0^{\dag}=\left( {\chi_1^0,\;\chi_2^0,\ldots,\;\chi_N^0} \right)=\left\{ {\chi_i^0={{{\partial {E_{\text{M}}}}} \left/ {{\partial {q_i}}} \right.}} \right\} \) is the vector of atomic electronegativities, \( {{\boldsymbol{\upeta}}_0}=\left\{ {\eta_{ij}^0={{{{\partial^2}{E_{\text{M}}}}} \left/ {{\partial {q_i}\partial {q_j}}} \right.}} \right\} \) is the hardness matrix, and N is the total number of atoms in the system. An index of “0” indicates that derivatives were evaluated for q = 0. Electronegativity represents the ability of an atom to attract electrons and, in the framework of DFT, is identified as the negative chemical potential [40]. It fulfills Sanderson’s electronegativity equalization principle, which states that for the equilibrium charge distribution, electronegativity is equalized throughout space:

$$ {\chi_1}={\chi_2}=\ldots ={\chi_N}={\chi^{\mathrm{eq}}}\equiv \chi . $$
(3)

In other words, electrons flow between atoms until the atomic electronegativities are equalized. Note that the electronegativity in an open system (global equilibrium) is equal to \( {\chi_i}=\chi_i^0+\sum\nolimits_j {\eta_{ij}^0{q_j}} \). In MD simulations, one often deals with molecular systems constructed from smaller subsystems (e.g., individual molecules). It is convenient to impose a charge conservation constraint on each of these subsystems. Electronegativity is then equalized in each of the fragments but not in the system as a whole.

By taking the derivative of Eq. 2 with respect to the i-th atomic charge and applying the EE principle (Eq. 3), one arrives at a set of linear equations for the equilibrium charge distribution in M. If M is constructed from two subsystems A and B, M = (A|B) (the solid line indicates no charge transfer between A and B), the CSA equations can be summarized as the following single matrix equation:

$$ \left( {\begin{array}{*{20}c} 0 & 0 & {{{\mathbf{1}}_{\text{A}}}} & {{{\mathbf{0}}_{\text{B}}}} \\ 0 & 0 & {{{\mathbf{0}}_{\text{A}}}} & {{{\mathbf{1}}_{\text{B}}}} \\ {\mathbf{1}_{\text{A}}^{\dag}} & {\mathbf{0}_{\text{A}}^{\dag}} & {\boldsymbol{\upeta}_{{\text{AA}}}^0} & {\boldsymbol{\upeta}_{{\text{AB}}}^0} \\ {\mathbf{0}_{\text{B}}^{\dag}} & {\mathbf{1}_{\text{B}}^{\dag}} & {\boldsymbol{\upeta}_{{\text{BA}}}^0} & {\boldsymbol{\upeta}_{{\text{BB}}}^0} \\ \end{array}} \right)\left( {\begin{array}{*{20}c} {-\chi_{\text{A}}^{{\text{eq}} }} \\ {-\chi_{\text{B}}^{{\text{eq}} }} \\ {{{\mathbf{q}}_{\text{A}}}} \\ {{{\mathbf{q}}_{\text{B}}}} \\ \end{array}} \right)=\left( {\begin{array}{*{20}c} {{q_{\text{A}}}} \\ {{q_{\text{B}}}} \\ {-\boldsymbol{\upchi}_{\text{A}}^0} \\ {-\boldsymbol{\upchi}_{\text{B}}^0} \\ \end{array}} \right), $$
(4)

where 1 A (1 B) and 0 A (0 B) are unity and zero vectors of dimension N A (N B), and \( \chi_{\text{A}}^{\text{eq} } \) and \( \chi_{\text{B}}^{\text{eq} } \) are the fragments’ equalized electronegativities. The hardness matrix η 0 and the vector that groups the atomic electronegativities χ 0 are divided into sub-blocks. The first two rows in Eq. 4 represent the charge conservation constraints for A and B on q A and q B, respectively. The remaining equations describe electronegativity equalization.

The equilibrium charge distribution in M is obtained by inverting the matrix in Eq. 4:

$$ \left( {\begin{array}{*{20}c} {-\chi_{\text{A}}^{{\text{eq}} }} \\ {-\chi_{\text{B}}^{{\text{eq}} }} \\ {{{\mathbf{q}}_{\text{A}}}} \\ {{{\mathbf{q}}_{\text{B}}}} \\ \end{array}} \right)=\left( {\begin{array}{*{20}c} {-{\eta_{{\text{AA}} }}} & {-{\eta_{{\text{AB}} }}} & {{{\mathbf{f}}_{{\text{AA}} }}} & {{{\mathbf{f}}_{{\text{AB}} }}} \\ {-{\eta_{{\text{BA}} }}} & {-{\eta_{{\text{BB}} }}} & {{{\mathbf{f}}_{{\text{BA}} }}} & {{{\mathbf{f}}_{{\text{BB}} }}} \\ {\mathbf{f}_{{\text{AA}}}^{\dag}} & {\mathbf{f}_{{\text{AB}}}^{\dag}} & {-{{\boldsymbol{\upbeta}}_{{\text{AA}} }}} & {-{{\boldsymbol{\upbeta}}_{{\text{AB}} }}} \\ {\mathbf{f}_{{\text{BA}}}^{\dag}} & {\mathbf{f}_{{\text{BB}}}^{\dag}} & {-{{\boldsymbol{\upbeta}}_{{\text{BA}} }}} & {-\boldsymbol{\upbeta} {}_{{\text{BB}} }} \\ \end{array}} \right)\left( {\begin{array}{*{20}c} {{q_{\text{A}}}} \\ {{q_{\text{B}}}} \\ {-\boldsymbol{\upchi}_{\text{A}}^0} \\ {-\boldsymbol{\upchi}_{\text{B}}^0} \\ \end{array}} \right) $$
(5)

The elements of the inverse matrix are: diagonal Fukui function (FF) indices, \( {{\mathbf{f}}_{{\text{AA}} }}=\left\{ {{{{\partial {q_{{i\in {\text{A}}}}}}} \left/ {{\partial {q_{\text{A}}}}} \right.}} \right\} \) and \( {{\mathbf{f}}_{{\text{BB}} }}=\left\{ {{{{\partial {q_{{i\in {\text{B}}}}}}} \left/ {{\partial {q_{\text{B}}}}} \right.}} \right\} \); off-diagonal Fukui function indices, \( {{\mathbf{f}}_{{\text{AB}} }}=\left\{ {{{{\partial {q_{{i\in {\text{A}}}}}}} \left/ {{\partial {q_{\text{B}}}}} \right.}} \right\} \) and \( {{\mathbf{f}}_{{\text{BA}} }}=\left\{ {{{{\partial {q_{{i\in {\text{B}}}}}}} \left/ {{\partial {q_{\text{A}}}}} \right.}} \right\} \); the condensed hardness matrix \( {{\boldsymbol{\upeta}}_{{\text{frg}} }}=\left\{ {{{{{\eta_{XY }}=\partial {\chi_X}}} \left/ {{\partial {q_Y},}} \right.}\;X,\,\,Y=\left( {{\text{A}},\;{\text{B}}} \right)} \right\} \), and the polarization matrix \( {{\boldsymbol{\upbeta}}_{XY }}=\left\{ {{{{\partial {q_{{i\in X}}}}} \left/ {{\partial {v_{{j\in Y}}}}} \right.}} \right\}, \) \( X,\,Y=\left. {\left( {{\text{A}},\;{\text{B}}} \right)} \right\} \). The Fukui function describes the response of a given atom to a perturbation in the total number of electrons in its fragment (diagonal FF) or another fragment (off-diagonal FF). Diagonal FF vectors are normalized to unity, f AA 1 A † = f BB 1 B † = 1, while off-diagonal FF vectors are normalized to zero, f AB 1 B † = f BA 1 A † = 0. The block (β XY ) of the polarization matrix represents the charge reorganization within X induced by perturbation in the external potential on the j-th atom of Y (v j∈Y ). Both the diagonal and off-diagonal blocks of the polarization matrix are normalized to zero:

$$ \sum\nolimits_{{i\in X}}^{{{N_X}}} {{{{\left( {{\beta_{XY }}} \right)}}_{ij }}} =0,\;\left( {X,Y} \right)\in \left\{ {{\text{A}},{\text{B}}} \right\}. $$
(6)

Atomic (diagonal) hardnesses describe an atom’s strength to block electron flow. They can be related to the hard and soft acids and bases (HSAB) principle [41]. In the frozen core approximation, the off-diagonal elements of the hardness matrix can be related to the Coulomb exchange integrals between valence shell s electrons:

$$ \eta_{\mathrm{ij}}^0\approx \left( {{s_i}{s_i}\left\| {{s_j}{s_j}} \right.} \right)=\left( {{s_i}{s_i}\left| {{s_j}{s_j}} \right.} \right)-\frac{1}{2}\left( {{s_i}{s_j}\left| {{s_j}{s_i}} \right.} \right)\equiv {J_{ij }}-\frac{1}{2}{K_{ij }} $$
(7)

They can be further approximated by neglecting the exchange integrals, K ij , and reducing η ij to the Coulomb integrals J ij . Therefore, one can determine the crudest approximation of \( \eta_{{ij}}^0 \) as 1/R ij . However, for short separations, the utilization of 1/R ij results in incorrect behavior, and semi-empirical combination formulae should be used instead. In accordance with our previous studies [42], we adopted Ohno’s formula [43]:

$$ \eta_{ij}^0=\frac{1}{{\sqrt{{a_{ij}^2+R_{ij}^2}}}}, $$
(8)

where \( {a_{ij }}={2 \left/ {{\left( {\eta_{ii}^0+\eta_{jj}^0} \right)}} \right.}\equiv {2 \left/ {{\left( {\eta_i^0+\eta_j^0} \right)}} \right.} \). Therefore, only diagonal hardnesses are needed to construct the hardness matrix.

The empirical parameters of CSA are the effective electronegativities and diagonal hardnesses of atoms in molecules. These depend on the atomic number, hybridization, and chemical environment of a given atom, and are determined by reproducing the gas-phase charge distribution calculated with an ab initio method for a set of training molecules. In our previous work [37, 42], we determined CSA parameters based on Mulliken population analysis (MPA) [44] and a hybridization/chemical environment classification inherited from the AMBER ff99 force field [38].

Model systems and computational details

An inclusion complex of β-cyclodextrin (BCD) with salicylic acid (SAL) was used to illustrate the performance of CSA. The BCD molecule contains seven glucopyranose residues which form a shallow truncated cone with primary and secondary hydroxyl groups exposed to the solvent. The primary hydroxyl groups form the smaller rim of the cone. The secondary hydroxyl groups form the bigger rim. The molecule contains a hydrophobic cavity and is able to incorporate many organic molecules. Inclusion complexes of BCD with different “guest” molecules are stabilized by weak interactions (van der Waals and hydrogen bonds). At the classical MD level of theory, each reactant preserves its identity, so complexes of BCD are good candidates to use when describing mutual polarization. In addition, the guest molecule breaks the high symmetry of the host (BCD) molecule. This leads to huge charge distribution inhomogeneities in glucopyranose residues and allows us to check the parameterization of CSA in reactant resolution.

BCD has several conformers; more information about their structures can be found in [45–49]. The conformer with a counterclockwise orientation of hydrogen bonds in both rims was chosen. The guest molecule was assumed to be in a planar conformation. The most stable conformer of SAL contains an intramolecular hydrogen bond between the phenyl hydroxyl hydrogen atom and the carbonyl oxygen atom from the carboxylic group. A less stable conformer was chosen for the MD calculation since, at first glance, it appears better suited to inter-reactant hydrogen-bond formation. The structures of the host and guest molecules are shown in Scheme 1.

Scheme 1
scheme 1

The host BCD molecule (a) together with the glucopyranose residue (b) and the guest SAL molecule (c)

Full size image

The initial structures of both molecules were obtained at the B3LYP/6-31G(d) level of theory using the Gaussian 09 [50] suite of programs. The starting geometry of the complex was obtained by placing the SAL molecule on the C7 symmetry axis of BCD, with the line joining the C1 and C4 aromatic carbons coinciding with the symmetry axis. This initial structure was the starting point in Born–Oppenheimer MD (BOMD) calculations. The TeraChem [51] general purpose quantum chemistry software package was used. A 10 ps MD run with 1 fs time steps at 300 K was performed at the B3LYP/6-31G level of theory. An empirical dispersion correction was utilized. The trajectory obtained was recalculated at the molecular mechanics (MM) level of theory. The CHARMM force field for hydrocarbons was employed. Data concerning pyranose monosaccharides and glycoside linkages were taken from [52] and [53]. In our previous work [37], we parameterized CSA using a diverse training set consisting of small organic molecules and AMBER ff99 all-atom force-field resolution. In the CHARMM force field, however, the sp 3 carbon atoms are additionally differentiated into secondary, primary, tertiary, and quaternary, as well as according to their chemical environment. None of these features were observed when force-field resolution was achieved with AMBER ff99. The classification of aliphatic hydrogen atoms also differed between the two force fields. In order to preserve the consistency of the CHARMM carbohydrate force field, we corrected the CSA parameters for BCD using the same protocol as before and a training set containing mono- and polyhydroxylated alcohols, simple linear and cyclic ethers, carbohydrates (aldoses, ketoses, furanoses, pyranoses, and disaccharides), and other model compounds. For the SAL molecule, defined using the CHARMM general force field for small molecules, the atom types were consistent with the previous CSA resolution, so the original parameter set was retained. All parameters are collected in Table 1. Three types of MM runs were performed using the NAMD package [54]. In the first MM run, the charge distribution was frozen. The Mulliken charges obtained at the B3LYP/6-31G(d) level were used. In the second and third MM runs, charges were derived using Eq. 5 for each point on the trajectory. In the second MM run, the charges in each subsystem were constrained (see Eq. 4). In the third run, CT between subsystems was allowed, and only the global charge of the BCD/SAL complex was constrained. In other words, two first rows/columns of Eq. 4 were replaced with one row/column: (0, 1 N ). These three types of calculations are termed the electrostatic (ES), polarization (P), and charge-transfer (CT) runs throughout the paper; these names reflect the physical effects that are additionally included in the MM calculations. The ES, P, and CT runs are characterized by the following equilibria: \( {{\text{M}}^{{\text{ES}} }}=\left( {\,{1_{\text{A}}}\left| {\,{2_{\text{A}}}} \right.\left| {\,\ldots } \right.\left| {\,{N_{\text{A}}}} \right.\left| {\,{1_{\text{B}}}\left| {\,{2_{\text{B}}}} \right.\left| {\,\ldots } \right.\left| {\,{N_{\text{B}}}} \right.} \right.} \right) \), \( {{\text{M}}^{\text{P}}}=\left( {\,{1_{\text{A}}}\vdots\,{2_{\text{A}}}\vdots\,\ldots \vdots\,{N_{\text{A}}}\left| {\,{1_{\text{B}}}\vdots\,{2_{\text{B}}}\vdots\,\ldots \vdots\,{N_{\text{B}}}} \right.} \right) \), and \( {{\text{M}}^{{\text{CT}} }}=\left( {{1_{\text{A}}}\vdots\,{2_{\text{A}}}\vdots\,\ldots \vdots\,{N_{\text{A}}}\vdots {1_{\text{B}}}\vdots\,{2_{\text{B}}}\vdots\,\ldots \vdots\,{N_{\text{B}}}} \right) \). Solid/dotted lines indicate that charge flow between atoms is forbidden/allowed.

Table 1 Effective hardness (in eV/e2; second column) and electronegativity (in eV/e; third column) data for CHARMM force-field atoms (first column)
Full size table

Results and discussion

A detailed analysis of the trajectory obtained allowed us to distinguish the following steps in the inclusion process: (i) the O7′ and H7′ atoms of SAL form a hydrogen bond with the O6 atom of BCD (inter-reactant hydrogen bond); (ii) the O2 and H2 atoms of SAL form a hydrogen bond with O2 of BCD (inter-reactant hydrogen bond); (iii) the O2 and H2 atoms of SAL form a hydrogen bond with O3 of BCD (inter-reactant hydrogen bond); (iv) O2, H2, and O7 of SAL form a hydrogen bond (intra-reactant hydrogen bond). Steps (ii) and (iii) have intermediate character. These inter-reactant hydrogen bonds appear to consecutively disappear after a conformational change in the SAL molecule and the formation of the intra-reactant hydrogen bond (step (iv)). The hydrogen bonds formed in steps (i) and (iv) are preserved for the rest of the simulation. The first is formed for t ∈ (700,900), while the second is formed for t ∈ (6400, 6600). Therefore, two main domains can be distinguished: one for t ∈ (900, 6400) and the other for t ∈ (6600, 10000).

Energy and charge profiles

Figure 1 shows the time evolution of the system’s energy. Note that the B3LYP/6-31G total energy (red curve) was shifted by a constant value in order to allow all of the curves to be plotted together. In general, the remaining curves should be shifted, since none of the force-field terms resulting from Taylor expansions include the zero-order contribution (e.g., bond stretching term). Two main domains can be distinguished along the trajectory. The average energy in the first domain (after step (i)) is higher than that in the second domain (after step (iv)). The relative positions of the curves on the energy scale agree with expectation: mutual polarization introduces additional stabilization into the system. It can be seen that the P (purple) and CT (blue, dotted line) curves are shifted down by about 100 kcal/mol with respect to the ES (green) curve. By removing the barrier to charge flow between BCD and the guest molecule, further stabilization is introduced. The CT stabilization doesn’t exceed 0.5 kcal/mol, so the P and CT curves are practically indistinguishable. Such energy lowering is typical in the hydrogen-bonded systems [55–57].

Fig. 1
figure 1

Variations in total energy versus time in the ab initio BOMD run (red curve), ES run (green curve), CT run (blue curve), and P run (purple curve) for the BCD/SAL complex

Full size image

The fluctuations in BOMD energy are smaller than those seen in the MM runs (ES, P, and CT). The MM curves differ only in the electrostatic contribution. The geometries at a given time step are exactly the same for all runs. Therefore, the bonding and van der Waals energies are also the same. This indicates that the bigger amplitudes seen for the P, CT, and ES runs are due to the bonding and van der Waals terms, which can be corrected for by adjusting some of the force-field parameters.

The shapes of ES, P, and CT curves are almost identical, and reproduce the main topological features of the BOMD curve. Figures S1 and S2 of the “Electronic supplementary material” (ESM) clearly demonstrate this, as they present the first and second derivatives of the system’s energy with respect to time for a narrow region of the trajectory, 6400 ≤ t ≤ 6600. This part of the trajectory was chosen since it corresponds to the hydrogen-bond formation that restores the most stable conformer of the SAL molecule (step (iv)). It is evident that the BOMD energy profile is very well reproduced. The positions of the extrema and the curvatures of all of the curves are almost identical.

A diagram showing the correlation between the B3LYP/6-31 and CSA charges is shown in Fig. 2. This figure includes all atoms for t ∈ (6400, 6600). The correlation coefficient is higher than 0.99. Our charges are slightly underestimated with respect to the B3LYP charges. Some of these discrepancies are due to basis set inconsistency: namely, the CSA model was parameterized for 6-31G(d) charges while, due to a lack of d-functions in the TeraChem package, they were compared with the 6-31G charges. The correlation diagram shows the average characteristic. It is important to determine whether the trends in charge reorganization accompanying hydrogen-bond formation are the same in both sets. Figures S3a and S3b of the ESM relate to the MPA and CSA charges, respectively. Charges from the CT run and the P run are indistinguishable, so only one of these sets (that from the CT run) was plotted. The main features of MPA charge evolution are preserved by CSA. Namely, O7 and H2 (the hydroxyl hydrogen atom) exhibit “step-like” population changes.

Fig. 2
figure 2

Correlation between CSA and MPA charges for the BCD/SAL complex, obtained for simulation frames corresponding to 6400 < t < 6600 fs

Full size image

CSA parameter profiles

Inclusion in the BCD cavity can modify the physicochemical properties of the guest molecules (solubility, diffusion, sublimation, volatility, chromatographic mobility, hydrophobicity, reactivity, etc.) [58]. The influence of BCD on the physicochemical properties of the guest molecule can be seen in Figs. 3 and 4, which show electronegativity and hardness data. Figure 3 presents the reactant electronegativities χ SAL (green curve) and χ BCD (red curve), together with the global electronegativity χ (blue curve). The most important observation that can be made from this figure is the fact that the electronegativity of the whole inclusion complex is almost the same as the electronegativity of the BCD molecule. A similar observation follows from Fig. 4. The hardness of the inclusion complex η (global hardness) is the same as the diagonal hardness η BCD,BCD. Such behavior demonstrates that the electronic properties of the inclusion complex are closer to those of the host BCD molecule. In other words, the presence of the guest molecule is masked and the SAL/BCD complex resembles the host molecule. Taking into account the bioadaptability of BCD [59], it is not surprising that BCD has found widespread application as an encapsulating material for medicaments and food ingredients, since the weak interactions stabilizing the complex allows the guest molecule to be released at the target location.

Fig. 3
figure 3

Time evolutions of the electronegativities χ SAL (green), χ BCD (red), and χ

Full size image
Fig. 4
figure 4

Time evolutions of η SAL,SAL (green), η BCD,BCD (red), η SAL,BCD (purple), and η (blue)

Full size image

The history of the inclusion process is reflected in the CSA data. The parameters most sensitive to bond formation are the polarization matrix elements and the diagonal FF indices. The diagonal FF indices and polarization matrix elements correspond to the P run [intra-reactant equilibrium: \( {{{\rm{M}}}^{{\rm{P}}}} = ({{1}_{{\rm{A}}}} \vdots \,{{2}_{{\rm{A}}}} \vdots \, \ldots \vdots \,{{N}_{{\rm{A}}}}\left| {{{1}_{{\rm{B}}}} \vdots \,{{2}_{{\rm{B}}}} \vdots \, \ldots \vdots \,{{N}_{{\rm{B}}}}} \right.) \)]. One can expect that the atoms involved in hydrogen-bond formation are the most sensitive to a perturbation in the external potential or in the number of electrons.

In Fig. 5, the population responses of the hydrogen-accepting oxygen to an external potential perturbation on the hydrogen atom (black curve) and the hydrogen-donating oxygen atom (red curve) are plotted. The polarization matrix is symmetric, so these elements also describe the population responses of the hydrogen and hydrogen-donating oxygen atoms to an external potential perturbation on the hydrogen-accepting oxygen atom (Maxwell relations). Parts A and B correspond to inter-reactant (step (i)) and intra-reactant (step (iv)) hydrogen bonds, respectively. One can observe step-like behavior by \( \beta_{\mathrm{H}7\text{',}\mathrm{O}6}^{\mathrm{SAL}\text{,}\mathrm{BCD}} \), \( \beta_{\mathrm{O}7\text{',}\mathrm{O}6}^{\mathrm{SAL}\text{,}\mathrm{BCD}} \) (part A) and \( \beta_{\mathrm{H}\text{2,}\mathrm{O}7}^{\mathrm{SAL}\text{,}\mathrm{SAL}} \), \( \beta_{\mathrm{O}\text{2,}\mathrm{O}7}^{\mathrm{SAL}\text{,}\mathrm{SAL}} \) (part B). The jump in the average value of \( \beta_{\mathrm{H}7\text{',}\mathrm{O}6}^{\mathrm{SAL}\text{,}\mathrm{BCD}} \) or \( \beta_{\mathrm{H}\text{2,}\mathrm{O}7}^{\mathrm{SAL}\text{,}\mathrm{SAL}} \) is greater than that of \( \beta_{\mathrm{O}7\text{',}\mathrm{O}6}^{\mathrm{SAL}\text{,}\mathrm{BCD}} \) or \( \beta_{\mathrm{O}\text{2,}\mathrm{O}7}^{\mathrm{SAL}\text{,}\mathrm{SAL}} \).

Fig. 5
figure 5

Time evolutions of selected elements of the polarization matrix. Part A presents \( \beta_{\mathrm{H}7\text{',}\mathrm{O}6}^{\mathrm{SAL}\text{,}\mathrm{BCD}} \) and \( \beta_{\mathrm{O}7\text{',}\mathrm{O}6}^{\mathrm{SAL}\text{,}\mathrm{BCD}}, \) while part B presents \( \beta_{\mathrm{H}\text{2,}\mathrm{O}7}^{\mathrm{SAL}\text{,}\mathrm{SAL}} \) and \( \beta_{\mathrm{O}\text{2,}\mathrm{O}7}^{\mathrm{SAL}\text{,}\mathrm{SAL}} \) (see Scheme 1 for atom labels)

Full size image

Figure 6 presents the evolutions of diagonal FF indices over time. Atoms involved in the formation of the intra-fragment hydrogen bond (namely O2, H2, and O7 of SAL) were taken into account. The diagonal FF index f i AA= ∂q i A/∂q A, measures the response of the i-th atom population of A (SAL) to the oxidation/reduction of A as a whole. One can observe step-like behavior in all of the curves (in an averaged sense). This behavior is especially apparent for H2. After hydrogen-bond formation, the atom becomes harder and the initial average FF value of 0.10 is halved to 0.05 (FF can be considered as the normalized softness and hardness is the inverse of the softness). For the remaining atoms, the jump in FF index is less pronounced. The hydrogen-donating atom (O2) becomes harder after bond formation, while the hydrogen-accepting oxygen (O7) becomes softer. The variations in the off-diagonal FF indices are shown in Fig. 7. The same atoms involved in the intra-fragment hydrogen bond of the SAL molecule are considered. Off-diagonal effects are rather subtle in comparison with those for the diagonal FF indices. This is not surprising, since they measure the response of the i-th atom of SAL to the oxidation/reduction of BCD. None of the curves show step-like behavior. The values of \( f_{\mathrm{O}2}^{\mathrm{SAL}\text{,}\mathrm{BCD}} \) and \( f_{\mathrm{O}7}^{\mathrm{SAL}\text{,}\mathrm{BCD}} \) gradually increase from one domain (its average value) to the other domain. \( f_{\mathrm{H}2}^{\mathrm{SAL}\text{,}\mathrm{BCD}} \) exhibits huge fluctuations that are damped after intra-fragment hydrogen-bond formation. These fluctuations were connected with inter-reactant hydrogen-bond formation (steps (ii) and (iii)). The off-diagonal FF indices carry the “chemical” information (hybridization, rehybridization, polarization, etc.), especially in local or molecular orbital resolution [60, 61].

Fig. 6
figure 6

Time evolutions of the diagonal FF indices of the O2 (cyan), H2 (black), and O7 (green) atoms for partitioning of the system, where A = SAL and B = BCD (see Scheme 1 for atom labels)

Full size image
Fig. 7
figure 7

Time evolutions of the off-diagonal FF indices of the O2 (cyan), H2 (black), and O7 (green) atoms for partitioning of the system, where A = SAL and B = BCD (see Scheme 1 for atom labels)

Full size image

The clear and sharp character of the changes in FF indices and polarization matrix elements indicates that they should be a very good descriptors for detecting bond formation during MD simulations. In particular, they are better indicators than atomic charges, which have a smoother transition between the two domains. Of course, atomic charges cannot be used as indicators in static force fields.

Conclusions and future prospects

In this paper, polarization and charge-transfer contributions were incorporated into MM calculations using a CSA approach. CSA is based on electronegativity equalization equations supplemented by charge conservation equations. The results obtained for the BCD/SAL inclusion complex indicate that our model properly describes the energetics and charge reorganization that accompany the insertion of a guest molecule into the BCD cavity. In addition, the reactivity descriptors built into the CSA formalism can be used as indicators of specific chemical processes that occur in the system. This method can be easily extended to systems comprising a huge number of molecules. Work aiming at incorporating CSA in MD packages is currently being carried out. It appears that CSA parameters can be applied to developed reactive force fields which allow to investigate chemical reactions.