Can we treat ab initio atomic charges and bond orders as conformation-independent electronic structure descriptors?

Abstract

Chemically meaningful descriptors of molecular electronic structure which can be calculated from the results of ab initio calculations in a well-defined manner play an important role in establishing a fundamental link between chemistry and physics. Bond orders (BOs) and effective atomic charges, both expressible in terms of the electron density matrix and the overlap matrix, comprise examples of such descriptors. In the present paper the dependence of BOs (calculated by NPA-Wiberg, Natiello–Medrano and Mayer methods) and atomic charges (calculated by Natural Population Analysis (NPA) and Löwdin schemes) on the conformation of DNA deoxyribonucleoside 2′-deoxycytidine is investigated. Strong correlations between the studied descriptors of molecular structure calculated by different methods are revealed and discussed. The atoms and bonds of 2′-deoxycytidine exhibiting the most pronounced dependence of their effective charges and bond orders on the conformation are identified.

---

Introduction

An exceptionally long-standing challenge in chemical physics is the search for the methods able to link the results of quantum-chemical calculations to ‘traditional’ chemical concepts; it is fair to say that a comment ‘With increasing availability of good…wave functions for molecules, a systematic procedure for obtaining maximum insight from such data has become desirable’ made by R. S. Mulliken more than 60 years ago¹ has not become irrelevant over the years. Some of the most fundamental examples of these ‘traditional’ chemical concepts include the bond order and atomic valence. These two, along with another widely-used intuitive notion of ‘atomic charge’, belong to a more general class of quantities, so-called ‘electronic population descriptors’ (a term first used in ref. 2, although with a slightly different focus), i.e., the quantities which can be derived from the first-order density matrix (1-RDM)(where Ψ is all-electron wavefunction of molecule and σ_i are the spin indices of the electrons), or its finite-basis representation.
where the summation indices μ and ν go over all atomic basis functions χ_μ present in the system (we will omit the spin indices σ, σ′ hereinafter for simplicity; this can be readily done in case of a closed-shell systems – the one considered throughout this paper, whereas in a more general treatment the ‘spin-up’ and ‘spin-down’ parts of 1-RDM are analyzed separately^3–5).

In case when all of the atomic basis functions (atomic orbitals χ_μ(r⃑)) are mutually orthogonal,⁶ diagonal elements of the density matrix sum up to the total number of electrons and one can introduce ‘the (average) number of electrons’ at atom A also known as its electronic population

and thus introduce an effective atomic charge q_A of the atomwhere Z_A is the charge of the nucleus of atom A expressed in atomic units.

Additionally, in the case when atomic basis functions form an orthogonal basis set, it is possible to define a measure for a covalent bond order between any pair of atoms, X and Y, known as the Wiberg's bond index (defined in the footnote on p. 1086 in ref. 7; see also discussion in ref. 6), as

By calculating the sum of the orders of the bonds formed by a given atom A towards all other atoms of the system one can define³ the ab initio valence of the atom in the molecule

(strictly speaking, the latter equality holds only in case of a closed-shell system with a single-determinant wave function – see discussion of eqn (5) in ref. 8 for details).

Although definitions (2)–(4) provide, at first sight, an unambiguous first-principles ground for the concepts of atomic charge, bond order and atomic valence, their practical application requires the use of orthogonal set of atomic basis functions which is rarely the case in practical calculations. Instead, a typical basis set used in ab initio calculations is usually built from the basis functions which are not orthogonal between different atoms so that these functions first need to be orthogonalized before applying any of eqn (2)–(4). Depending on the orthogonalization method used, the resulting set of orthogonal basis functions (which we will refer to as ‘atomic orbitals’ (AOs) for convenience) can be quite different thus leading to different numerical values of the electronic population descriptors (EPDs) derived from the density matrix D. Moreover, since any orthogonalization procedure depends on the elements of the basis set overlap matrix

which, in turn, depend on mutual distance between the atoms at which the μ-th and ν-th basis functions are located, the EPDs are inherently dependent on the conformation of the molecule (and in some cases, even on the orientation of the molecule as a whole⁹).

In the present paper we investigate the dependence of atomic charges, bond orders and atomic valences derived from the first-order density matrix on the conformation of the molecule and compare this dependence for two orbital orthogonalization schemes: Löwdin-orthogonalization^10,11 and the NPA-procedure used in the Natural Bond Orbital (NBO) analysis family of methods^12–14 during the stage known as the natural population analysis¹⁵ (see also ref. 16 for a recent review). When applied to calculation of bond orders, the former of these two strategies of orbital orthogonalization leads to the Natiello–Medrano bond orders¹⁷

(where S^1/2DS^1/2 = D^L is the Löwdin density matrix¹⁸) while the latter leads to NPA-Wiberg bond orders given by eqn (3) with the density matrix D evaluated in the basis of Natural Atomic Orbitals (NAOs).¹⁵ As an alternative approach, we also consider the bond ordersintroduced by I. Mayer^3,8,19 with no explicit relation to any particular orbital orthogonalization schemes; nonetheless, the Mayer bond orders are not independent of the conformation of the molecule since they still depend on the overlap matrix S, right as the other bond orders do. Note also that the Mayer's approach does not suggest any method for atomic charge evaluation.

Although in case when the basis functions are mutually orthogonal (the overlap matrix S is identity matrix) the Natiello–Medrano, NPA-Wiberg and the Mayer's bond orders reduce to the Wiberg's definition given by eqn (3), but as soon as the basis functions have non-zero overlap (off-diagonal elements of S are non-zero) all the three methodologies are expected to give different values.

The investigation of conformational variability of EPDs takes its importance from the fact that atomic charges are essential parameters in nearly all force fields used in molecular dynamics simulations and their proper assignment is often a non-trivial task.^20–22 The situation becomes even more challenging due to the fact that for proper representation of ion–nucleic acid base interactions,²³ experimental temperature dependence of peptide helix formation,²⁴ G-DNA stem–loop junction ions stabilization²⁵ and many other biologically relevant phenomena^26,27 it is necessary to carefully account for polarization effects, i.e., to model the redistribution of the electron density of atomic fragments in response to an external influence or, more generally, in response to the changes in the environment (which vary from almost nonconductive inside protein cavities to very conductive on protein–water interfaces²⁸). The most common, whereas not the only possible, source of such an influence is electric field created by another molecule or atomic group in proximity. When the simulation involves a molecule which can adopt multiple conformations, just as already the simplest biomolecules do,^29–37 both atomic groups, one of which is subjected to an ‘external’ perturbation and the other one which creates the perturbation, can actually belong to the same molecule. In this case the information on the dependence of EPDs on the conformation can help in discriminating between the atoms and/or their groups which are more and less susceptible to the changes in their environment and thus provide some rationale for optimization of computational overhead²⁸ caused by inclusion of polarization effects into the force field.

Even if polarization effects are neglected, special care is still required in the selection of appropriate conformations for calculation of atomic charges during the force field parameterization.^38,39

Thus, the information on conformational dependence of EPDs can be considered as additional guidance in the process of force fields improvement or (re)parameterization.

Quantification of conformational variability of EPDs is also important for achieving better transferability of localized orbital which is necessary for developing new and assessing existing fragment-based electronic structure methods.^40,41

Computation details

As a model compound for the study we consider 2′-deoxycytidine (dC, Fig. 1), a biologically important molecule^29–34 which can be considered as the simplest model^35–37 of the DNA monomeric unit.

Fig. 1 The molecules of 2′-deoxycytidine (dC) and 2′-deoxycytidine-5′-monophosphate (dCMP) in their DNA-like conformations (A- and B-forms of the DNA are denoted by the letter in parentheses). Also shown are the atom numbering (in the upper row for backbone atoms and in the lower row for the nucleotide base) and conformational parameters conventions used through this study.

In addition to its biological importance, this molecule has a rich family of different conformations, as shown by Zhurakivsky and Hovorun by an exhaustive conformational analysis performed at DFT B3LYP/6-31G(d,p) level of theory.³⁴ They have found as many as 89 conformers of the molecule and presented the values of characteristic torsion angles for all of them.

Basing on these parameters we reconstructed the structures of all conformations and verified them by geometry optimization at B3LYP/6-31G(d,p) level of theory using Psi4 (beta 5) software package⁴² with default settings. The final Cartesian coordinates of all 89 structures are available in the ESI.†

Molden-like files obtained from Psi4 were processed by molden2molden program from JANPA package^16,43 to convert them from Cartesian into pure spherical harmonics of basis functions and achieve a proper normalization of basis functions. The final MOLDEN files were analyzed by janpa program^16,43 (version 1.04_r2) and Löwdin and NPA charges for all atoms have been obtained independently for each conformer of the molecule.

Wiberg bond orders and the corresponding atomic valences were calculated according to eqn (3) and (4) using built-in capabilities of the janpa program operating on the density matrix converted to the basis of the natural atomic orbitals by the NPA procedure. Natiello–Medrano and Mayer bond orders were calculated according to eqn (5) and (6) respectively with a home-made program written in Python (its source code is given in ESI†). The developed program was also used for statistical and graphical analyses of the obtained results.

Since deoxyribonucleosides are biologically significant compounds which, being linked by a phosphate group, form the strands of the DNA, it is of certain interest to specifically consider those ‘DNA-like’ geometries of 2′-deoxycytidine molecule which have the conformation similar to the conformation of the nucleosides within the DNA in its A- or B-forms (see Table 1). Interestingly, such DNA-like conformers can indeed be identified among the whole family of 89 conformers of an isolated 2′-deoxycytidine molecule.³⁴ The values of conformational parameters of these conformers of dC are given in Table 1 and Cartesian coordinates are available in ESI† (A-DNA-like conformation – structure dC_007; B-DNA-like conformation – structure dC_053).

Table 1 Conformational parameters of A- and B-DNA-like structures of the model DNA monomers^a

Conformational parameter^b, deg.	dC (A)	dC (B)	dCMP (A)	dCMP (B)	Experimental^c (A)	Experimental^c (B)
a dC(A/B) – 2′-deoxycytidine, dCMP(A/B) – 2′-deoxycytidine-5′-monophosphate; the letter in parentheses distinguishes between A- and B-forms of the DNA; the values of conformational parameters for dC and dCMP were obtained from the structures optimized at DFT B3LYP/6-31G(d,p) level of theory. Cartesian coordinates of the corresponding structures are given in ESI.b The torsion angles α, β, γ, ε, ζ, χ (possible values range from −180° to 180°) are defined in Fig. 1, while the sugar pseudorotation angle^44 P (0° ≤ P ≤ 360°) was defined, along with a ‘puckering amplitude’ νmax, by equations νmax cos P = ν2, νmax sin P = (ν1 − ν0 + ν4 − ν3)/(2(sin(36°) + sin(72°))), where νi are the sugar ring endocyclic torsion angles (defined in Fig. 1).c Typical values of the conformational parameters observed in crystalline DNA in it's A- and B-forms – see ref. 45–48 and references therein.
α	—	—	−44.4	−45.8	−69 ± 14	−62 ± 19
β	178.2	178.8	177.2	−172.4	175 ± 15	168 ± 25
γ	54.7	53.9	58.6	56.2	57 ± 10	51 ± 16
ε	−172.4	170.0	−171.1	171.9	−155 ± 14	−173 ± 20
ζ	—	—	−104.2	−106.1	−73 ± 9	−98 ± 23
χ	−163.5	−145.6	−163.8	−144.9	−161 ± 11	−108 ± 23
P	11.9	163.5	14.4	153.6	18 ± 18 (C3′ endo)	162 ± 18 (C2′ endo)

---

Moreover, not only 2′-deoxycytidine, but an isolated 2′-deoxycytidine-5′-monophosphate (dCMP, Fig. 1), a corresponding 2′-deoxyribonucleotide, was also shown⁴⁹ to have DNA-like conformations within its full family of 613 conformations. Since there exist DNA-like structures of dCMP with the conformational parameters really close to those of dC (see Table 1), it is possible to investigate the influence of the phosphate group, which links the neighbouring DNA monomers, on the charge (re)distribution in DNA-like conformers of 2′-deoxycytidine. To this end, quantum-chemical analysis similar to those outlined above for dC has been carried out for A- and B-DNA-like conformers of dCMP and the obtained results were compared.

Finally, since biological molecules do not function as closed systems, but continuously interact with other molecules, it is reasonable to investigate the possible influence of the environment (which can be considered as another kind of a ‘perturbing factor’) on EPDs.

To this end, additional set of single-point DFT calculations using integral equation formalism polarizable continuum (IEFPCM) model (see ref. 50 for review) were carried out at B3LYP/6-31G(d,p) level of theory in order to investigate the influence of dielectric continuum environment on atomic charges of DNA-like conformations of dC and dCMP molecules. Continuum environment parameters were set to those of water since a dielectric constant of this solvent is high enough (static permittivity ε_s = 78.39) to produce a notable mutual solute–solvent polarization and thus cause a non-negligible redistribution of electron density of the solvated molecule. These calculations were carried out using a PCMSolver library⁵¹ (version 1.1.3) integrated into Psi4 software⁴² (version 1.0.0). Most of PCMlibrary input parameters had their default values except for those defining the solute cavity construction process. The construction was done by GePol method (see ref. 50 and 52 and references therein) with the average tesserae area of the surface partition set to 0.2 a.u. and using UFF radii⁵³ with no additional scaling factor.

Results and discussion

Effective atomic charges

Statistical characteristics of the obtained Löwdin and NPA charges such as the mean values and mean absolute deviation (MAD) from these values have been calculated for each atom and are presented in Table 2.

Table 2 Mean NPA and Löwdin charges of 2′-deoxycytidine molecule atoms averaged over all conformations (q_av), and their mean absolute deviations (MAD)^a

Atom^b	q^NPAav	q^Lowdinav	MAD^NPA	MAD^Lowdin
a 89 conformations^34 optimized at B3LYP/6-31G(d,p) level of theory.b See Fig. 1 for atom numbering scheme.
N1	−0.467	−0.047	0.004	0.004
C2	0.805	0.152	0.007	0.006
O2	−0.644	−0.305	0.007	0.005
N3	−0.590	−0.226	0.002	0.002
C4	0.449	0.124	0.002	0.001
N4	−0.823	−0.232	0.002	0.003
H41	0.430	0.192	0.001	0.001
H42	0.415	0.183	0.001	0.001
C5	−0.403	−0.209	0.004	0.003
H5	0.249	0.108	0.001	0.001
C6	0.070	0.006	0.004	0.003
H6	0.251	0.117	0.009	0.004
C1′	0.280	0.112	0.004	0.004
H1′	0.242	0.098	0.019	0.011
C2′	−0.527	−0.193	0.005	0.006
H2′(1)	0.261	0.108	0.013	0.009
H2′(2)	0.268	0.113	0.011	0.007
C3′	0.070	0.028	0.002	0.003
H3′	0.236	0.093	0.015	0.011
C4′	0.037	0.022	0.006	0.008
H4′	0.232	0.089	0.012	0.008
O4′	−0.595	−0.249	0.008	0.006
O3′	−0.758	−0.308	0.006	0.006
HO3′	0.481	0.205	0.005	0.003
C5′	−0.122	−0.045	0.004	0.004
H5′(1)	0.216	0.089	0.012	0.009
H5′(2)	0.215	0.088	0.013	0.009
O5′	−0.762	−0.315	0.008	0.005
HO5′	0.485	0.205	0.007	0.004

---

To analyze the variance in the atomic charges of particular atoms depending on the conformation of 2′-deoxycytidine, we depicted these charge values on a Löwdin charge–NPA charge plane (Fig. 2).

Fig. 2 Correlation between Löwdin and NPA charges: blue markers correspond to the values averaged over all conformations while the orange markers represent the charge values in all of 89 investigated conformers of 2′-deoxycytidine.

It should be noted that apart from N1, N4 and C2 atoms (perhaps, affected by electron conjugation involving the C2=O2 bond), there is a strong correlation between Löwdin and NPA charges, and NPA charge is typically 2.6 times larger than the Löwdin charge. This can be attributed to the fact that NAO transformation leads to a stronger orbital localization as compared to the basis function transformation used in Löwdin population analysis. At the same time, NPA charges are a little more conformationally sensitive since their MADs are typically about 1.7 times higher than MADs of Löwdin charges. However, only H1′, H2′(1), H2′(2), H3′, H4′, H5′(1) and H5′(2) atoms have NPA MADs above 0.01e in the molecule under study. All of these atoms are hydrogens attached to sp³-hybridized carbons. On the other hand, H41, H42, the hydrogens of cytosine base amino group, have much smaller NPA MADs (typically ∼0.001e), and H_O3′, H_O5′, the hydrogens attached to oxygens, show intermediate NPA MADs (∼0.005e). It is worth noting that there is a distinction between H5 and H6 NPA MADs. While both of these hydrogens are attached to a sp²-hybridized carbon, the latter participates in hydrogen bonding (CH⋯O, where O4′ or O5′ atoms act as a proton acceptor) in 63 of 89 conformers,⁵⁴ while the former does not participate in hydrogen bonding in any of the conformers. Therefore, it can be a charge transfer phenomenon caused by the hydrogen bonding that is responsible for several times higher NPA MAD of H6 as compared to H5.

In addition to the above mentioned general trends in conformational dependence of effective atomic charges, it is of certain interest to focus on biologically important DNA-like conformers of dC and dCMP and consider how do the change in environment (from vacuum to continuum dielectric with permittivity ε_s = 78.39) and the addition of the phosphate group influence the atomic charges of these conformers. To this end, we present in Table 3 (Löwdin charges) and 4 (NPA charges) the values of effective atomic charges of A- and B-DNA-like conformations (shown in Fig. 1; see Table 1 for their conformational parameters and ESI† for Cartesian coordinates) of dC and dCMP calculated in vacuum and with the influence of continuum dielectric environment. For dCMP these tables contain also Σ_Ph – the sum of atomic charges of the phosphate group which substitutes the H_O5′ atom of dC (Tables 3 and 4).

Table 3 Löwdin charges of atoms of A- and B-DNA-like conformers of 2′-deoxycytidine (dC) and 2′-deoxycytidine-5′-monophosphate (dCMP) obtained in vacuum (vac) and under the influence of continuum dielectric environment (sol)

Atom^a	dC-A (vac)	dC-A (sol)	dC-B (vac)	dC-B (sol)	dCMP-A (vac)	dCMP-A (sol)	dCMP-B (vac)	dCMP-B (sol)
a See Fig. 1 for atom numbering scheme.b ΣPh – sum of charges of the dCMP phosphate group atoms P, OP, OP1, HP1, OP2, HP2 (applicable to dCMP only; the dC columns contain the charge of the HO5′ atom).
N1	−0.052	−0.047	−0.044	−0.041	−0.053	−0.047	−0.044	−0.041
C2	0.144	0.142	0.148	0.145	0.144	0.143	0.149	0.146
O2	−0.314	−0.371	−0.307	−0.366	−0.315	−0.370	−0.306	−0.364
N3	−0.230	−0.273	−0.229	−0.273	−0.230	−0.273	−0.229	−0.272
C4	0.122	0.125	0.122	0.127	0.123	0.127	0.124	0.128
N4	−0.239	−0.227	−0.237	−0.224	−0.235	−0.222	−0.234	−0.220
H41	0.189	0.197	0.190	0.197	0.190	0.198	0.191	0.198
H42	0.181	0.202	0.181	0.202	0.182	0.202	0.182	0.203
C5	−0.216	−0.203	−0.212	−0.197	−0.209	−0.200	−0.207	−0.194
H5	0.105	0.123	0.105	0.125	0.111	0.124	0.107	0.126
C6	0.007	0.012	0.012	0.016	0.005	0.012	0.013	0.017
H6	0.124	0.127	0.121	0.127	0.124	0.129	0.122	0.130
C1′	0.116	0.115	0.112	0.108	0.117	0.116	0.111	0.108
H1′	0.108	0.106	0.110	0.102	0.110	0.108	0.112	0.104
C2′	−0.195	−0.194	−0.197	−0.194	−0.194	−0.193	−0.197	−0.193
H2′(1)	0.107	0.109	0.107	0.112	0.108	0.111	0.109	0.114
H2′(2)	0.113	0.112	0.109	0.111	0.114	0.114	0.110	0.112
C3′	0.024	0.031	0.032	0.036	0.024	0.031	0.031	0.037
H3′	0.088	0.092	0.084	0.094	0.089	0.094	0.087	0.096
C4′	0.034	0.032	0.037	0.034	0.033	0.033	0.037	0.037
H4′	0.094	0.097	0.104	0.103	0.097	0.103	0.107	0.110
O4′	−0.243	−0.253	−0.251	−0.265	−0.244	−0.252	−0.252	−0.263
O3′	−0.309	−0.325	−0.314	−0.334	−0.307	−0.322	−0.312	−0.331
HO3′	0.211	0.226	0.206	0.220	0.213	0.227	0.208	0.222
C5′	−0.040	−0.035	−0.045	−0.038	−0.019	−0.015	−0.024	−0.019
H5′(1)	0.090	0.094	0.082	0.092	0.113	0.115	0.106	0.114
H5′(2)	0.088	0.094	0.087	0.091	0.108	0.114	0.107	0.112
O5′	−0.318	−0.330	−0.323	−0.335	−0.364	−0.369	−0.369	−0.374
P	—	—	—	—	1.064	1.080	1.062	1.079
OP1	—	—	—	—	−0.420	−0.421	−0.421	−0.422
OP	—	—	—	—	−0.563	−0.607	−0.560	−0.606
OP2	—	—	—	—	−0.403	−0.409	−0.406	−0.410
HP2	—	—	—	—	0.244	0.261	0.244	0.261
HP1	—	—	—	—	0.243	0.258	0.243	0.258
HO5′/ΣPh^b	0.212	0.227	0.211	0.226	0.165	0.163	0.162	0.159

---

Table 4 NPA charges of atoms of A- and B-DNA-like conformers of 2′-deoxycytidine (dC) and 2′-deoxycytidine-5′-monophosphate (dCMP) obtained in vacuum (vac) and under the influence of continuum dielectric environment (sol)

Atom^a	dC-A (vac)	dC-A (sol)	dC-B (vac)	dC-B (sol)	dCMP-A (vac)	dCMP-A (sol)	dCMP-B (vac)	dCMP-B (sol)
a See Fig. 1 for atom numbering scheme.b ΣPh – sum of charges of the dCMP phosphate group atoms P, OP, OP1, HP1, OP2, HP2 (applicable to dCMP only; the dC columns contain the charge of the HO5′ atom).
N1	−0.473	−0.468	−0.466	−0.464	−0.475	−0.468	−0.468	−0.465
C2	0.797	0.807	0.803	0.812	0.797	0.808	0.804	0.813
O2	−0.645	−0.706	−0.637	−0.701	−0.646	−0.706	−0.637	−0.699
N3	−0.594	−0.642	−0.594	−0.642	−0.594	−0.642	−0.594	−0.642
C4	0.446	0.452	0.447	0.455	0.446	0.453	0.449	0.456
N4	−0.828	−0.825	−0.827	−0.822	−0.825	−0.821	−0.824	−0.819
H41	0.426	0.434	0.427	0.435	0.427	0.436	0.428	0.436
H42	0.412	0.438	0.412	0.439	0.415	0.440	0.414	0.440
C5	−0.410	−0.401	−0.406	−0.395	−0.405	−0.398	−0.401	−0.392
H5	0.245	0.268	0.245	0.269	0.253	0.269	0.248	0.271
C6	0.070	0.072	0.072	0.072	0.068	0.073	0.071	0.072
H6	0.271	0.275	0.262	0.271	0.264	0.270	0.259	0.269
C1′	0.282	0.282	0.275	0.275	0.281	0.281	0.274	0.274
H1′	0.261	0.257	0.264	0.253	0.263	0.260	0.266	0.256
C2′	−0.530	−0.529	−0.530	−0.529	−0.530	−0.529	−0.530	−0.529
H2′(1)	0.257	0.259	0.258	0.265	0.258	0.261	0.261	0.268
H2′(2)	0.269	0.267	0.261	0.264	0.270	0.269	0.260	0.263
C3′	0.067	0.071	0.071	0.072	0.068	0.073	0.070	0.073
H3′	0.229	0.234	0.226	0.238	0.229	0.235	0.228	0.240
C4′	0.047	0.044	0.046	0.045	0.047	0.045	0.046	0.045
H4′	0.237	0.241	0.254	0.253	0.241	0.250	0.258	0.261
O4′	−0.592	−0.602	−0.598	−0.611	−0.592	−0.600	−0.598	−0.609
O3′	−0.761	−0.784	−0.764	−0.792	−0.760	−0.783	−0.764	−0.790
HO3′	0.488	0.509	0.484	0.504	0.490	0.511	0.486	0.506
C5′	−0.121	−0.121	−0.123	−0.122	−0.129	−0.129	−0.131	−0.131
H5′(1)	0.219	0.224	0.208	0.221	0.247	0.250	0.237	0.247
H5′(2)	0.217	0.224	0.217	0.221	0.242	0.249	0.242	0.248
O5′	−0.775	−0.795	−0.777	−0.798	−0.867	−0.873	−0.869	−0.875
P	—	—	—	—	2.568	2.593	2.566	2.591
OP1	—	—	—	—	−1.031	−1.039	−1.031	−1.039
OP	—	—	—	—	−1.063	−1.108	−1.061	−1.107
OP2	—	—	—	—	−1.021	−1.036	−1.024	−1.037
HP2	—	—	—	—	0.532	0.554	0.532	0.554
HP1	—	—	—	—	0.532	0.551	0.532	0.551
HO5′/ΣPh^b	0.491	0.513	0.491	0.512	0.517	0.516	0.514	0.513

---

The presented data indicate that the charge of the O2 atom is the most sensitive to the changes in the environment: it changes by −0.06 (independently of the charge calculation method, Löwdin or NPA) when dC or dCMP in its DNA-like conformation is put from vacuum into continuum dielectric environment. The changes in the charge of N3 and O_P (in dCMP) atoms closely follow the same trend, but are slightly smaller in absolute value. The most notable increases in atomic charges are observed for H42 and H_P2 atoms, but in general it can be said that the total 0.1e (0.15e in dCMP) transferred to O2, N3 and O_P (in dCMP) atoms upon the polarization of the molecule in dielectric environment originate roughly uniformly from many other atoms, i.e., it is hardly possible to distinguish an atom which increases its charge more significantly than the others.

The change of conformation of dC and dCMP from A- to B-DNA-like, which is virtually the flipping of the sugar ring from C3′-endo into C2′-endo form, leaves most of the atomic charges almost unchanged with the exception of H4′ atom. Its NPA charge increases by 0.016e when the molecule is switched from A- to B-DNA-like conformation in vacuum.

Finally, attachment of the phosphate group to dC noticeably influences only the charge of O5′ atom, the one to which it is actually being attached. Its Löwdin charge changes by about −0.04 upon attachment, while NPA charge changes twice as much. Interestingly, the sum of charges of the phosphate group atoms P, O_P, O_P1, H_P1, O_P2, H_P2 in dCMP is almost the same as the NPA charge of H_O5′ atom of dC, so it can be speculated that except for the absence of phosphate group nucleosides and nucleotides are quite similar in terms of their electron density distribution.

Bond orders and atomic valences

A general idea about the distribution of the bond order values calculated by the three different methods can be deduced from the data in Fig. 3. As it could be anticipated in advance, an arbitrary selected pair of atoms is more likely to be not bonded with a covalent bond rather than bonded. This simple idea explains the presence of a strong peak near zero value of the bond order in Fig. 3a (notice a semi-logarithmic scale). However, if the bond orders lower than 0.5 are excluded from the statistics, it becomes clear that the remaining bond orders are mainly concentrated near the value of 1.0 corresponding to a single bond (Fig. 3b). Two features of the latter statistics are worth mentioning. First, independently of the method, the value 0.5 can be taken as a threshold for distinguishing between covalently bound and non-bound pairs of atoms (the averaged bond orders exceeding the mentioned threshold are given in Table 5 along with the mean absolute deviations (MADs) from the averaged values). Second, all the three methods, Mayer, NPA-Wiberg and Natiello–Medrano, give comparable variance of bond order values in the region near 1.0.

Fig. 3 Statistical distribution of Mayer, Natiello–Medrano and NPA-Wiberg bond orders ((a), semi-log scale) and a magnified fragment (in range from 0.5 to 2.5) of the distributions (b).

Table 5 Statistical characteristics of bond orders between all pairs of covalently bound atoms of 2′-deoxycytidine molecule^a calculated by NPA-Wiberg, Natiello–Medrano and Mayer methods: average values and mean absolute deviations (MAD) from the average

Bond^b	Averages			MAD × 100
	NPA-Wiberg	Natiello–Medrano	Mayer	NPA-Wiberg	Natiello–Medrano	Mayer
a 89 conformations^34 optimized at B3LYP/6-31G(d,p) level of theory.b See Fig. 1 for atom numbering scheme.
N1–C2	0.958	1.092	0.904	0.50	0.53	0.72
C2–O2	1.615	2.024	1.708	0.99	1.06	1.50
C2–N3	1.209	1.389	1.191	0.43	0.47	0.53
N3–C4	1.446	1.620	1.436	0.39	0.45	0.45
C4–N4	1.208	1.358	1.097	0.56	0.57	0.44
N4–H41	0.794	0.946	0.882	0.12	0.07	0.07
N4–H42	0.812	0.958	0.899	0.10	0.06	0.05
C4–C5	1.176	1.210	1.176	0.61	0.60	0.59
C5–H5	0.907	0.921	0.937	0.08	0.07	0.13
C5–C6	1.601	1.655	1.603	0.98	0.97	1.46
C6–H6	0.898	0.903	0.920	0.85	0.40	0.95
N1–C6	1.208	1.359	1.160	0.55	0.78	0.63
C1′–N1	0.911	1.014	0.886	0.91	1.18	1.31
C1′–H1′	0.885	0.899	0.929	1.02	0.26	0.54
C1′–C2′	0.991	1.015	0.960	0.18	0.22	0.34
C2′–H2′(1)	0.901	0.927	0.940	0.63	0.26	0.32
C2′–H2′(2)	0.898	0.927	0.938	0.73	0.50	0.88
C2′–C3′	0.988	1.023	0.971	0.26	0.43	0.92
C3′–H3′	0.894	0.908	0.937	0.59	0.26	0.62
C4′–C3′	0.975	0.993	0.942	0.38	0.48	1.02
C4′–H4′	0.892	0.905	0.933	0.41	0.28	0.42
C4′–O4′	0.883	1.089	0.865	0.59	0.92	1.33
O4′–C1′	0.922	1.141	0.922	0.78	1.21	2.49
C3′–O3′	0.932	1.168	0.911	0.69	1.03	1.31
O3′–HO3′	0.753	1.008	0.877	0.41	0.28	0.32
C4′–C5′	0.993	1.039	0.968	0.58	0.80	1.60
C5′–H5′(1)	0.912	0.928	0.937	0.47	0.25	0.64
C5′–H5′(2)	0.912	0.928	0.936	0.43	0.26	0.54
C5′–O5′	0.939	1.178	0.946	0.92	1.21	1.52
O5′–HO5′	0.744	1.005	0.867	1.45	0.98	1.59

---

To get a deeper insight into numerical interconnections between the bond orders (only the ones above 0.5 were considered) obtained by the three methods, the values of bond orders obtained by all possible combinations of calculation methods were plotted in Fig. 4 and approximated with linear functions using the least-squares method.

Fig. 4 Correlations between Mayer, Natiello–Medrano and NPA-Wiberg bond orders.

Above all, these plots reveal that different methods lead to bond orders which are strongly correlated with one another. At the same time, the numerical values of the different bond orders are far from being equal. Indeed, as it has been determined by a direct comparison, the NPA-Wiberg bond orders (B_NPA-Wiberg) are almost always smaller than the Natiello–Medrano ones (B_{Natiello–Medrano}) and inequality

B_NPA-Wiberg < B_{Natiello–Medrano}

holds for more than 98,7% of all pairs of covalently bound atoms in all investigated conformations of 2′-deoxycytidine. A similar inequality holds for Mayer (B_Mayer) and Natiello–Medrano bond orders

B_Mayer < B_{Natiello–Medrano}

in 70% of all of the covalently bound pairs of atoms. Interestingly, no similar inequality between B_Mayer and B_NPA-Wiberg can be established since the values of these bond orders are always quite close and the difference between their values never exceeded 0.132. It is also worth noting that Natiello–Medrano bond orders are typically slightly higher than those expected from chemically intuition.

Although a visual inspection of the plots shown in Fig. 4 can hardly reveal any differences in the degree to which a particular bond order depends on the conformation of the molecule, the data presented in Table 5 reveals that the mean absolute deviation of the bond orders of individual covalent bonds are quite close for Natiello–Medrano and NPA-Wiberg methods, whereas the method by Mayer leads to about 1.5 times higher values of MADs (and thus to more conformationally sensitive bond orders) and there are only four cases (N1–C6, C4–C5, C4–N4 and N4–H42 bonds) when the Mayer bond orders exhibits the smallest conformational dependence among the studied methods. The highest value of MAD both Natiello–Medrano (0.012) and Mayer (0.025) bond orders achieve for O4′–C1′ bond, while the NPA-Wiberg scheme distinguishes the bond O5′–H_O5′ as the one having the most conformationally-sensitive bond order (its MAD equals 0.014).

To get a general idea of the maximum possible extent to which a bond order can change depending on a conformation of the molecule, we also calculated the differences between the largest and the smallest values of each bond order between a covalently bound pair of atoms observed in all investigated conformations of 2′-deoxycytidine. The results indicate that according to the NPA-Wiberg method, the bond order of O5′–H_O5′ bond changes, depending on conformation, from 0.667 (conformation dC_047, see ESI†) to 0.766 (dC_074), and no other bond has higher difference between its maximum and minimum bond orders higher. In contrast to this, two other methods, Natiello–Medrano and Mayer, showed that it is C2=O2 bond whose bond order changes more than anyone else: for this bond Natiello–Medrano method gives the maximum bond order of 2.056 (dC_006) vs. minimum of 1.972 (dC_039), while Mayer's method leads to the maximum value of 1.752 (dC_051) vs. minimum of 1.644 (dC_039) for the same bond. It can be easily noted that these individual, extreme deviations are considerably higher than any of MADs.

We also analyzed some ‘integral characteristics’ deduced from bond orders and obtained a quantitative measure of their dependence on the conformation of the molecule. A natural choice for such characteristics is the atomic valence defined by eqn (4).

We evaluated the atomic valences for each atom of 2′-deoxycytidine molecule in all of its conformations using each of three methods of bond orders calculations. Table 6 contains the results of statistical analysis of these valences and in this respect is conceptually similar to Table 5 mentioned previously. In the presented data the dependence of atomic valences on conformation is characterized by mean absolute deviation of atomic valences of each atom from the corresponding values averaged over all conformations.

Table 6 Statistical characteristics of atomic valences of 2′-deoxycytidine molecule^a atoms derived from bond orders calculated by NPA-Wiberg, Natiello–Medrano and Mayer methods: average values and mean absolute deviations (MAD) from the average

Atom^b	Averages			MAD
	NPA-Wiberg	Natiello–Medrano	Mayer	NPA-Wiberg	Natiello–Medrano	Mayer
a 89 conformations^34 optimized at B3LYP/6-31G(d,p) level of theory.b See Fig. 1 for atom numbering scheme.
N1	3.417	4.072	3.102	0.003	0.003	0.006
C2	3.883	4.810	3.907	0.002	0.006	0.006
O2	1.978	2.515	1.974	0.006	0.014	0.008
N3	3.047	3.568	2.824	0.002	0.002	0.003
C4	3.953	4.579	3.835	0.000	0.001	0.002
N4	3.027	3.571	2.983	0.006	0.007	0.005
H41	0.818	1.047	0.918	0.001	0.001	0.001
H42	0.829	1.046	0.914	0.001	0.001	0.001
C5	3.915	4.209	3.867	0.001	0.001	0.007
H5	0.940	1.036	0.964	0.001	0.000	0.000
C6	3.882	4.318	3.836	0.008	0.004	0.019
H6	0.940	1.043	0.962	0.004	0.008	0.009
C1′	3.810	4.430	3.787	0.016	0.015	0.038
H1′	0.946	1.046	0.956	0.009	0.003	0.003
C2′	3.880	4.205	3.811	0.004	0.003	0.011
H2′(1)	0.935	1.039	0.958	0.007	0.002	0.002
H2′(2)	0.931	1.042	0.962	0.006	0.005	0.004
C3′	3.861	4.396	3.789	0.007	0.006	0.012
H3′	0.948	1.044	0.956	0.007	0.004	0.004
C4′	3.849	4.384	3.739	0.006	0.009	0.017
H4′	0.951	1.043	0.956	0.006	0.002	0.002
O4′	2.000	2.627	1.884	0.009	0.009	0.010
O3′	1.790	2.390	1.853	0.008	0.010	0.014
HO3′	0.772	1.076	0.904	0.005	0.004	0.004
C5′	3.812	4.272	3.811	0.009	0.010	0.019
H5′(1)	0.956	1.041	0.951	0.006	0.003	0.004
H5′(2)	0.957	1.041	0.950	0.006	0.004	0.004
O5′	1.782	2.374	1.880	0.008	0.012	0.013
HO5′	0.767	1.082	0.901	0.007	0.013	0.006

---

It is worth noting that the Natiello–Medrano method gives the values of atomic valence which always exceed the values anticipated based on chemical intuition. At the same time, it can be seen that Natiello–Medrano valences typically demonstrate weaker sensitivity to the molecule conformation than the ones obtained by the Mayer's method and only four atoms, H2′(2), H_O5′, N4 and O2, present an exception from this rule. Apart from that, all three methods perform similarly.

Plots shown in Fig. 5 indicate that the values of atomic valences are interrelated in the same way as their ‘parent’ bond orders. Moreover, it has been found that inequalities V_NPA-Wiberg < V_{Natiello–Medrano} and V_Mayer < V_{Natiello–Medrano} hold with no exception, and that the difference between V_Mayer and V_NPA-Wiberg never exceeds 0.33. Similarly to bond orders, atomic valences obtained by the three methods exhibit strong correlations with one another (see Fig. 5).

Fig. 5 Correlation between Mayer, NPA-Wiberg and Natiello–Medrano atomic valences.

To summarize, all three methods show similar extent of conformational dependence of bond orders, whereas the Natiello–Medrano method produces slightly less conformationally dependent values. However, the latter method typically overestimates bond orders as compared to the idealized values anticipated from chemical intuition.

Conclusions

We analyzed the conformational dependence of NPA-Wiberg, Natiello–Medrano and Mayer bond orders, as well as NPA and Löwdin charges for 89 conformations of 2′-deoxycytidine molecule, a model DNA monomeric unit. Strong correlations between bond orders, as well as between the charges, obtained by the mentioned methods were observed. It has been found that depending on conformation of the molecule the bond orders of covalently bound pairs of atoms can change by up to 0.1. Natiello–Medrano method produces slightly less conformationally dependent bond orders whereas the absolute values of the bond orders produced by this method are typically 10% higher than the values anticipated from chemical intuition. It has been shown that NPA atomic charge is typically 2.6 times larger than the Löwdin charge by its absolute value, and that the mean absolute deviations (MADs) of NPA charges from their values averaged over all conformations of 2′-deoxycytidine are typically about 1.7 times higher than MADs of Löwdin charges thus indicating that NPA charges are slightly more conformation-dependent. Thus, the use of more localized atomic orbitals as basis functions for population analysis leads to more conformation-dependent values of electronic population descriptors.

Acknowledgements

The authors are grateful to Dr Roman Zhurakivsky for providing atomic coordinates of the 2′-deoxycytidine molecule conformers.

Notes and references

1. R. S. Mulliken, Electronic population analysis on LCAO–MO molecular wave functions. I., J. Chem. Phys., 1955, 23, 1833–1840.
2. D. Cortés-Arriagada and G. I. Cárdenas-Jirón, A through-space charge transfer mechanism for explaining the oxidation of 2-chlorophenol on a tetrasulphonated nickel(III) phthalocyanine, Comput. Theor. Chem., 2011, 963, 161–167.
3. I. Mayer, Charge, bond order and valence in the ab initio SCF theory, Chem. Phys. Lett., 1983, 97, 270–274.
4. E. R. Davidson, Properties and uses of natural orbitals, Rev. Mod. Phys., 1972, 44, 451–464.
5. E. R. Davidson, Reduced Density Matrices in Quantum Chemistry (Theoretical chemistry: a series of monographs, v. 6), Academic Press Inc., 1976.
6. M. S. Gopinathan and K. Jug, Valency. I. A quantum chemical definition and properties, Theor. Chim. Acta, 1983, 63, 497–509.
7. K. B. Wiberg, Application of the pople-santry-segal CNDO method to the cyclopropylcarbinyl and cyclobutyl cation and to bicyclobutane, Tetrahedron, 1968, 24, 1083–1096.
8. I. Mayer, Bond order and valence indices: A personal account, J. Comput. Chem., 2007, 28, 204–221.
9. I. Mayer, Löwdin population analysis is not rotationally invariant, Chem. Phys. Lett., 2004, 393, 209–212.
10. P. O. Löwdin, On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals, J. Chem. Phys., 1950, 18, 365–375.
11. B. C. Carlson and J. M. Keller, Orthogonalization procedures and the localization of Wannier functions, Phys. Rev., 1957, 105, 102.
12. F. Weinhold, Natural bond orbital methods, in Encyclopedia of Computational Chemistry, ed. P. v. R. Schleyer, N. L. Allinger, T. Clark, J. Gasteiger, P. A. Kollman, H. F. Schaefer III and P. R. Schreiner, Chichester, John Wiley & Sons, 1998, vol. 3, pp. 1792–1811.
13. E. D. Glendening, C. R. Landis and F. Weinhold, Natural bond orbital methods, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 1–42.
14. A. E. Reed, L. A. Curtiss and F. Weinhold, Intermolecular interactions from a natural bond orbital, donor–acceptor viewpoint, Chem. Rev., 1988, 88, 899–926.
15. A. E. Reed, R. B. Weinstock and F. Weinhold, Natural population analysis, J. Chem. Phys., 1985, 83, 735–746.
16. T. Yu. Nikolaienko, L. A. Bulavin and D. M. Hovorun, JANPA: An open source cross-platform implementation of the Natural Population Analysis on the Java platform, Comput. Theor. Chem., 2014, 1050, 15–22.
17. M. A. Natiello and J. A. Medrano, On the quantum theory of valence and bonding from the ab initio standpoint, Chem. Phys. Lett., 1984, 105, 180–182.
18. T. Kar, A. B. Sannigrahi and D. C. Mukherjee, Comparison of atomic charges, valencies and bond orders in some hydrogen-bonded complexes calculated from Mulliken and Löwdin SCF density matrices, J. Mol. Struct.: THEOCHEM, 1987, 153, 93–101.
19. I. Mayer, Bond orders and valences from ab initio wave functions, Int. J. Quantum Chem., 1986, 29, 477–483.
20. J. A. Lemkul, W. J. Allen and D. R. Bevan, Practical considerations for building GROMOS-compatible small-molecule topologies, J. Chem. Inf. Model., 2010, 50, 2221–2235.
21. B. A. Horta, P. F. Fuchs, W. F. van Gunsteren and P. H. Hünenberger, New interaction parameters for oxygen compounds in the GROMOS force field: Improved pure-liquid and solvation properties for alcohols, ethers, aldehydes, ketones, carboxylic acids, and esters, J. Chem. Theory Comput., 2011, 7, 1016–1031.
22. M. M. Reif, P. H. Hünenberger and C. Oostenbrink, New interaction parameters for charged amino acid side chains in the GROMOS force field, J. Chem. Theory Comput., 2012, 8, 3705–3723.
23. S. Nakagawa, Polarizable model potential function for nucleic acid bases, J. Comput. Chem., 2007, 28, 1538–1550.
24. J. Huang and A. D. MacKerell, Induction of Peptide Bond Dipoles Drives Cooperative Helix Formation in the (AAQAA)₃ Peptide, Biophys. J., 2014, 107, 991–997.
25. J. Song, C. Ji and J. Z. H. Zhang, The critical effect of polarization on the dynamical structure of guanine quadruplex DNA, Phys. Chem. Chem. Phys., 2013, 15, 3846–3854.
26. Y. Zhong and S. Patel Binding, Structures of Tri-N-Acetyl-β-Glucosamine in Hen Egg White Lysozyme Using Molecular Dynamics with a Polarizable Force Field, J. Comput. Chem., 2013, 34, 163–174.
27. J. A. Lemkul, J. Huang, B. Roux and A. D. MacKerell, An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications, Chem. Rev., 2016, 116, 4983–5013.
28. W. Wang and R. D. Skeel, Fast evaluation of polarizable forces, J. Chem. Phys., 2005, 123, 164107-1–164107-12.
29. S. P. Samijlenko, I. V. Alexeeva, L. H. Palchykivs'ka, I. V. Kondratyuk, A. V. Stepanyugin, A. S. Shalamay and D. M. Hovorun, ¹H NMR investigation on 6-azacytidine and its derivatives, Spectrochim. Acta, Part A, 1999, 55, 1133–1141.
30. Y. R. Mishchuk, A. L. Potyagaylo and D. M. Hovorun, Structure and dynamics of 6-azacytidine by MNDO/H quantum-chemical method, J. Mol. Struct., 2000, 552, 283–289.
31. Y. P. Yurenko, R. O. Zhurakivsky, M. Ghomi, S. P. Samijlenko and D. M. Hovorun, Comprehensive Conformational Analysis of the Nucleoside Analogue 2′-β-Deoxy-6-azacytidine by DFT and MP2 Calculations, J. Phys. Chem. B, 2007, 111, 6263–6271.
32. Y. P. Yurenko, R. O. Zhurakivsky, M. Ghomi, S. P. Samijlenko and D. M. Hovorun, Ab initio comprehensive conformational analysis of 2′-Deoxyuridine, the biologically significant DNA minor nucleoside, and reconstruction of its low-temperature matrix infrared spectrum, J. Phys. Chem. B, 2008, 112, 1240–1250.
33. Y. P. Yurenko, R. O. Zhurakivsky, S. P. Samijlenko and D. M. Hovorun, Intramolecular CH⋯O hydrogen bonds in the AI and BI DNA-like conformers of canonical nucleosides and their Watson–Crick pairs. Quantum chemical and AIM analysis, J. Biomol. Struct. Dyn., 2011, 29, 51–65.
34. R. O. Zhurakivsky and D. M. Hovorun, Comprehensive conformational analysis of canonical nucleoside 2′-deoxycytidine by density functional theory, Phys. Alive, 2006, 14(3), 33–46 , in Ukrainian, https://www.researchgate.net/publication/305496270_Comprehensive_conformational_analyses_of_canonical_nucleoside_2′-deoxycytidine_by_density_functional_theory.
35. N. Foloppe and A. D. MacKerell Jr, Intrinsic conformational properties of deoxyribonucleosides: implicated role for cytosine in the equilibrium among the A, B, and Z forms of DNA, Biophys. J., 1999, 76, 3206–3218.
36. T. Yu Nikolaienko, L. A. Bulavin and D. M. Hovorun, Structural flexibility of DNA-like conformers of canonical 2′-deoxyribonucleosides, Phys. Chem. Chem. Phys., 2012, 14, 15554–15561.
37. T. Y. Nikolaienko, L. A. Bulavin and D. M. Hovorun, How flexible are DNA constituents? The quantum-mechanical study, J. Biomol. Struct. Dyn., 2011, 29, 563–575.
38. A. D. MacKerell, Empirical force fields for biological macromolecules: overview and issues, J. Comput. Chem., 2004, 25, 1584–1604.
39. P. E. M. Lopes, O. Guvench and A. D. MacKerell, Current status of protein force fields for molecular dynamics simulations, in Molecular Modeling of Proteins, ed. A. Kukol, Springer, New York, 2015, ch. 3, pp. 47–71.
40. Z. Li, H. Li, B. Suo and W. Liu, Localization of Molecular Orbitals: From Fragments to Molecule, Acc. Chem. Res., 2014, 47, 2758–2767.
41. B. Meyer, B. Guillot, M. F. Ruiz-Lopez and A. Genoni, Libraries of Extremely Localized Molecular Orbitals. 1. Model Molecules Approximation and Molecular Orbitals Transferability, J. Chem. Theory Comput., 2016, 12, 1052–1067.
42. J. M. Turney, A. C. Simmonett, R. M. Parrish, E. G. Hohenstein, F. Evangelista, J. T. Fermann, B. J. Mintz, L. A. Burns, J. J. Wilke, M. L. Abrams, N. J. Russ, M. L. Leininger, C. L. Janssen, E. T. Seidl, W. D. Allen, H. F. Schaefer, R. A. King, E. F. Valeev, C. D. Sherrill and T. D. Crawford, Psi4: An open-source ab initio electronic structure program, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 556–565.
43. janpa: A cross-platform open-source implementation of NPA with Java; A part of JANPA package, http://janpa.sourceforge.net.
44. C. Altona and M. Sundaralingam, Conformational Analysis of the Sugar Ring in Nucleosides and Nucleotides. A New Description Using the Concept of Pseudorotation, J. Am. Chem. Soc., 1972, 94, 8205–8212.
45. O. P. Boryskina, M. Yu. Tkachenko and A. V. Shestopalova, Variability of DNA structure and protein-nucleic acid recognition, Biopolym. Cell, 2010, 26(5), 360–372.
46. D. Svozil, J. Kalina, M. Omelka and B. Schneider, DNA conformations and their sequence preferences, Nucleic Acids Res., 2008, 36, 3690–3706.
47. H. Berman, Crystal studies of B-DNA: The answers and the questions, Biopolymers, 1997, 44, 23–44.
48. B. Schneide, St. Neidle and H. M. Berman, Conformations of the Sugar–Phosphate Backbone in Helical DNA Crystal Structures, Biopolymers, 1997, 42, 113–124.
49. T. Yu. Nikolaienko and D. M. Hovorun, Quantum-mechanical conformational analysis of 2′-deoxycytidilic acid molecule – the DNA structural unit, Dopov. Nac. akad. nauk Ukr., 2010, 9, 173–184 , in Ukrainian, http://dopovidi.nas.gov.ua/2010-09/10-09-29.pdf.
50. J. Tomasi, B. Mennucci and R. Cammi, Quantum mechanical continuum solvation models, Chem. Rev., 2005, 105, 2999–3094.
51. PCMSolver, v1.1.3, an Application Programming Interface for the Polarizable Continuum Model electrostatic problem, written by R. Di Remigio, L. Frediani and K. Mozgawa with contributions from R. Bast, J. Juselius and V. Weijo, see http://pcmsolver.readthedocs.io/.
52. J. L. Pascual-ahuir, E. Silla and I. Tunon, GEPOL: An improved description of molecular surfaces. III. A new algorithm for the computation of a solvent-excluding surface, J. Comput. Chem., 1994, 15, 1127–1138.
53. These radii are actually the halves of the ‘atomic van der Waals distances’, ‘x_I’, introduced in Table I in a paper, A. K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard III and W. M. Skiff, UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations, J. Am. Chem. Soc., 1992, 114, 10024–10035.
54. R. O. Zhurakivsky and D. M. Hovorun, Intramolecular hydrogen bonds in conformers of 2′-deoxycytidine: results of quantum-chemical analysis of electron density topology, Ukr. Biochem. J., 2006, 78(6), 70–77.

---

Footnote

† Electronic supplementary information (ESI) available: A Python code for statistical and graphical analysis of bond orders properties and Cartesian coordinates of 89 conformers of 2′-deoxycytidine molecule. See DOI: 10.1039/c6ra17055b

---