Density functional theory calculations of the geometry and vibrational frequencies. Equilibrium geometries of the investigated G·C base pairs and transition states (TSs) of their tautomerizations, as well as their harmonic vibrational frequencies have been calculated, using Gaussian’09 program package  at the B3LYP/6-311 + + G(d,p) level of theory [26–30], which approved itself successfully for the calculations of the similar systems and processes and shown acceptable level of accuracy and adequacy of the obtained results [30, 31]. A scaling factor that is equal to 0.9668 has been applied in the present work for the correction of the harmonic frequencies for all complexes [20, 21, 32].
We have confirmed local minima and transition states, localized by Synchronous Transit-guided Quasi-Newton method , on the potential energy landscape by the absence or presence, respectively, of one imaginary frequency in the vibrational spectra of the complexes.
All reaction pathways have been reliably confirmed by providing intrinsic reaction coordinate (IRC) calculations  from each TS in the forward and reverse directions at the B3LYP/6-311 + + G(d,p) level of theory.
All calculations have been performed in the continuum with ε = 1, that adequately reflects the processes occurring in real biological systems without deprivation of the structurally-functional properties of the bases in the composition of the DNA or RNA molecules and satisfactorily models the substantially hydrophobic recognition pocket of the DNA-polymerase machinery as a part of the replisome [34, 35].
Single point energy calculations. We continued geometry optimizations with electronic energy calculations as the single point calculations at the MP2/6-311 + + G(2df,pd) level of theory [36, 37].
The Gibbs free energy G for all structures was obtained in the following way:
G = Eel+Ecorr, (1)
where Eel – electronic energy and Ecorr – thermal correction.
QTAIM analysis. Bader's quantum theory of Atoms in Molecules (QTAIM)  was applied to analyse the electron density distribution, using program package AIMAll .
The presence of the bond critical point (BCP), namely the so-called (3,-1) BCP, and a bond path between the donor and acceptor of the H-bond, as well as the positive value of the Laplacian at this BCP (Δρ > 0), were considered as criteria for the formation of the H-bond [40–43]. Wave functions were obtained at the B3LYP/6-311 + + G(d,p) level of theory, used for geometry optimisation.
The atomic numbering scheme for the bases is conventional and rare tautomeric forms of the G and C bases are marked by an asterisk (*) .
Obtained results and their discussion.
In this work investigated tautomerization pathways of the G·C base pairs are presented on Fig. 1 and their discussion is outlined below.
It is interesting to note that G·C*tO2(WC) base pair (Fig. 1) tautomerizes through the double proton transfer along the N1H…N3 and O2H…N2 H-bonds and via the TSG·C*tO2(WC)↔G*NH3·C*t(WC), which are stabilized by the participation of the two intermolecular (С)N4H…O6(G) and (C)O2H…N2(G) H-bonds and (G)N1-H-N3(C) covalent bridge. Eventually, this reaction leads to the formation of the G*NH3·C*t(WC) base pair, stabilized by three intermolecular N4H…O6, N3H…N1 and N2H…O2 H-bonds. Exactly the proton transfer along the lower O2H…N2 H-bond leads to the formation of the NH3 group at the G base.
Formed G*NH3·C*t(WC) base pair can transform via the mutual rotation of the bases around the middle N3H…N1 H-bond, leading to the new reverse G*NH3·C*t(rWC) base pair, stabilized by the N3H…O6, N4H…N1 and N2H…N4 H-bonds. Transition state of this interconversion TSG*NH3·C*t(WC)↔G*NH3·C*t(rWC) is joined by three N4H…N1, N3H…N1, N2H…N4 H-bonds and single N2…N3 van der Waals contact.
Another reaction – G·C*O2(WC)↔G*NH3·C*(WC) – occurs via the transfer of the proton, localized at the N1 nitrogen atom of the G base, to the N3 nitrogen atom of the C*O2 base and of the proton, localized at the O2 oxygen atom of the C*O2 base to the N2 atom of the NH2 amino group of the G base and finally leads to the G*NH3·C*(WC) base pair by the participation of the G*NH3 base with NH3 group. Transition state TSG·C*O2(WC)↔G*NH3·C*(WC) of this reaction is characterized by the (G)O6…N4(C) van der Waals contact and two (G)N1-H-N3(C) and (G)N2-H-O2(C) covalent bridges.
Formed G*NH3·C*(WC) base pair can transform by the mutual rotations of the bases around the middle N3H…N1 H-bond into the reverse G·C(rwWC) base pair: G*NH3·C*(WC)↔G·C(rwWC). Transition state TSG*N3·C*(WC)↔G·C(rwWC) of this reaction is joined by three intermolecular N3H…N1, N2H…N4, N2H…O2 H-bonds and N2…N3 van der Waals contact.
The most interesting case represents the last transformation – G*·C*O2(WC)↔G*NH3·C(wWC)↓, since proton transfer within the G*·C*O2(WC) base pair leads not only to the changing of its tautomeric status, but also to its geometrical rearrangement. This reaction occurs through the proton transfer along the upper O6H…N4 and lower O2H…N2 H-bonds from the O2 atom of the G* base to the N4 atom of the C*O2 base and from the O2 atom of the C*O2 base to the N2 atom of the G* base, respectively. Finally, C base shifts down accordingly the G*NH3 base, forming the G*NH3·C(wWC)↓ base pair by the participation of the G base with NH3 group.
Altogether it was revealed four G·C base pairs, involving G*NH3 base with NH3 group – G*NH3·C*t(WC), G*NH3·C*t(rWC), G*NH3·C*(WC) and G*NH3·C(wWC)↓ (Fig. 1).
In general, considered here G·C base pairs form the following order in terms of their relative Gibbs free ΔG and electronic ΔE energies (in kcal·mol− 1): G·C(rwWC) (0.00 and 0.00) < G*NH3·C(wWC)↓ (16.08 and 14.54) < G·C*tO2(WC) (19.43 and 17.53) < G*·C*O2(WC) (19.44 and 17.56) < G*NH3·C*t(WC) (24.60 and 22.32) < G*NH3·C*(WC) (31.67 and 30.17) < G*NH3·C*t(rWC) (32.97 and 31.80).
Notably, that difference in Gibbs free and electronic energies between the classical G·C(WC) and reverse G·C(rwWC) base pairs consists 11.53 and 13.09 kcal·mol− 1, respectively, while the G·C*tO2(WC) and G·C*O2(WC) base pairs are iso-energetical (19.43 and 19.44 kcal·mol− 1, respectively).