Density-Functional Based Tight-Binding : an Approximate DFT Method

O método DFTB, bem como a sua extensão com carga corrigida auto-consistente SCC-DFTB, tem ampliado a faixa de aplicações das ferramentas teóricas com fundamentos bem estabelecidos. Como uma aproximação do método do funcional de densidade, o método DFTB mantém aproximadamente a mesma precisão, mas com custo computacional menor, permitindo a investigação da estrutura eletrônica de sistemas grandes que não podem ser explorados com métodos ab initio convencionais. No presente artigo, os fundamentos dos métodos DFTB, SCC-DFTB e da inclusão das forças de dispersão de London são revisados. Para mostrar um exemplo da aplicabilidade do método DFTB, o equilíbrio zwitteriônico de glicina em solução aquosa é investigado. Foram realizadas simulações de dinâmica molecular usando o hamiltoniano SCC-DFTB corrigido para incluir a dispersão e uma caixa periódica contendo 129 moléculas de água, a partir de uma abordagem puramente mecânico-quântica.


Introduction
Density functional theory (DFT) methods are the standard and the most used theoretical techniques for electronic structure calculations. 1-5The advent of the generalized gradient approximation (GGA) for the exchange-correlation functional enhanced the DFT accuracy 6 and the predicted molecular structures, relative energies and frequencies are nearly comparable to the second order Møller-Plesset perturbation theory (MP2) method, with remarkable success to treat transition metal complexes. 7Efficient algorithms to solve the Kohn-Sham equations have been implemented, scaling to N 3 with respect to the size of the basis sets and, hence, being much more efficient than the N 5 of the MP2 methods.DFT is the method chosen for a huge range of applications.The formalism of the DFT and its extension to the reactivity indexes are subject of intensive research and many empirical concepts such as electronegativity, chemical potential and hardness are now formally defined within the DFT framework. 3,5,8-10 With respect to the methodology, developments concerning improved exchange-correlation functionals and hybrid quantum-mechanics/molecular-mechanics (QM/MM) methods are still the main subjects of research of many theoreticians. 11Chemical property estimates based on DFT are now well established, and even optical properties are accessible through the generalization to time-dependent DFT, 7,12,13 a method which is nowadays implemented in many different computer codes.
Notwithstanding the marvelous ability of the DFT to treat systems of increasing complexity, many systems are still intractable at the actual stage of computer technology development.Biosystems, adsorption processes, nanostructures, molecular dynamics, clusters and aggregates with thousands of atoms, self-assembling systems, nanoreactors and supramolecular chemistry are some of the fields in which ab initio methods cannot be used with adequate chemical models.For this range of systems, semi-empirical methods seem to have their applicability.Semi-empirical methods such as AM1, 14 PM3 15-17 and, more recently, RM1 18 have many empirical parameters that are fitted to a set of molecular properties, estimated either theoretically or experimentally.Therefore, the applicability of such methods is restricted.
Density-functional tight-binding (DFTB) is an approximate method based on the density functional framework which does not require large amounts of empirical parameters.The virtues and weaknesses of the DFTB are a heritage from DFT.In fact the parameters are consistently obtained from DFT calculations of a few molecules per pair of atom types.On the other side, DFTB is closely connected to the tight-binding methods.In fact, it can be seen as a non-orthogonal tight-binding method parameterized from DFT.The self-consistent charge extension of DFTB (SCC-DFTB) improves very much the accuracy of the method.For improvement of physical approximations, all DFT extensions, such as treatment of relativistic effects and London dispersion, can be easily used in the DFTB method.Large number of applications has been reported showing its usefulness in the calculations of hyperfine coupling constants, magnetic properties, vibrational spectra of solids and molecules, nuclear magnetic shielding tensors, geometries, dynamic properties and many others. 19-24 Calculation of optical properties is also possible due to the time-dependent DFTB, 25-27 which is not covered in the present paper.
The goal of the present review is to call the attention of the chemistry community to the DFTB method, which can be a good complement of the set of semi-empirical methods available.Its advantages and weaknesses are highlighted.As an example of application, the zwitterionic and neutral forms of glycine in aqueous solution are discussed in terms of fully quantum mechanical molecular dynamics of this molecule in water.

Background Fundaments
Density functional theory has been extensively reviewed. 7,28In this section a very brief review of DFT is done in order to highlight its crucial aspects to the formulation of the DFTB method.
The Hohenberg-Kohn (HK) theorems 29 have rigorously made the electronic density acceptable as basic variable to electronic-structure calculations.However, development of practical DFT methods only became relevant after W. Kohn and L. J. Sham published their famous set of equations: the so-called Kohn-Sham (KS) equations. 30 The use of the electronic density within the KS scheme allows a significant reduction of the computational demand involved in quantum calculations.Furthermore, the KS method paved the way for studying systems that could not be investigated by conventional ab initio methods (which use the wave function as basic variable).
Even though DFT methods have been successfully applied for systems of increasing complexity, methods which can include approximations to reduce even more the computational demand, without compromising the reliability of results, are still required.
The application of tight-binding (TB) to the calculation of electronic structures starts with the paper by J. Slater and G. Koster. 31The main idea behind this method is to describe the Hamiltonian eigenstates with an atomic-like basis set and replace the Hamiltonian with a parameterized Hamiltonian matrix whose elements depend only on the internuclear distances (this requires the integrals of more than two centers to be neglected) and orbital symmetries.
Although the Slater-Koster method was conceived for the calculation of band structures in periodic systems, it was later generalized to an atomistic model, capable of treating finite systems as well.The transition to atomistic has three main requirements, as discussed by Goringe et al. 32  First, the elements of the Hamiltonian matrix must have a functional dependence on the interatomic distance.In the case of band structures one just has to know the matrix elements for discrete values of distance.This requirement was solved by Froyen and Harrison, 33 who proposed that the interatomic distance was related to the Hamiltonian elements by 1/r 2 .
The second requirement is to obtain an expression for the total energy and not only for the band energy.In 1979 Chadi 34 proposed that the total energy could be described as a sum of two contributions, where E bnd is the sum over the energies of all occupied orbitals obtained by diagonalization of the parameterized Hamiltonian matrix, and E rep is the repulsive contribution, obtained by the sum of the atomic-pair terms, (2) in which N is the number of atoms in the system.
The third and last requirement is the possibility to derive the atomic forces from the total energy.This is especially important for geometry optimization and molecular dynamics.By assuming differentiability of U ab in equation 2, the only problem is to derive E bnd , which depends on the parameterization method chosen for the Hamiltonian matrix.
The DFTB method attends these three requirements with the additional advantage of completely avoiding any empirical parameterization, since the Hamiltonian and overlap matrices are calculated using atom-like valenceorbitals which are derived from DFT.Therefore, the DFTB method can be considered as a simplification of the Kohn-Sham method.

The Kohn-Sham Method
Although the Hohenberg and Kohn theorems 29 have proven that the electronic energy of a system can be totally determined from its electronic density through the variational principle, they did not propose any procedure to perform this calculation.This was done about one year later, by Kohn and Sham, 30 with the publication of their equations known as Kohn-Sham equations.
The solution of Kohn and Sham starts from the idea of using monoelectronic orbitals to calculate the kinetic energy in a simple, yet reasonably precise, way leaving a residual correction that could be calculated apart.Thus, one starts with a reference system of M non-interacting electrons subjected to the external potential n S , with Hamiltonian (3) Where (4) in which there are no electron-electron repulsion terms and for which the electronic density is exactly the same as in the corresponding system of interacting electrons.By introducing the single particle orbitals ψ i all electronic densities physically acceptable for the system of non-interacting electrons can be written in the form (5) Therefore, the HK functional can be written as (6) where T S represents the kinetic-energy functional for the reference system of M non-interacting electrons, given by (7) J represents the classic Coulomb interaction functional (8) and the remaining interactions are grouped in E xc , the exchange-correlation functional, which contains the difference between the exact kinetic energy T and T S , besides the non-classic part of the electron-electron interactions V ee , i.e. (9) After combining equations 6, 7 and 8 within the second HK theorem, the chemical potential can be written as (10) with the KS effective potential \ (11) where n ext is the external potential, normally due to the atomic nuclei, and the exchange-correlation potential n xc is defined as (12) Equation 10, restricted by ∫r(r  )dr  = M, is exactly the same equation that would be obtained for a system of M noninteracting electrons submitted to the external potential n S = n KS .Thus, for a given n KS a suitable value of r can be calculated for equation 10 by solving the M monoelectronic equations (13) and by using the calculated ψ i in equation 5.
Equations 5 and 11-13 are the so-called Kohn-Sham equations.Since n KS depends on r through n xc the KS equations must be solved iteratively using a self-consistent procedure similar to the one depicted in Figure 1.An electronic density model r 0 is normally chosen to start the iterative procedure.In principle, any positive function normalized for the number of electrons would be applicable, but a good initial estimate of r can significantly accelerate convergence.
At the end of the iterative procedure, the total energy can be calculated, which is given in the KS method by the following expression: (14) The most difficult part of the KS scheme is to calculate n xc in equation 12.The existence of an exact density functional is assured by the first HK theorem, but the exact form of the E xc functional remains unknown.However, many approximations of this functional have been described in the scientific literature over the last 30 years.In practice, the approximation chosen for E xc and the way by which the KS orbitals are represented define the different DFT methods.

DFT as Basis for a Tight-Binding Method
Following Foulkes and Haydock 35 the electronic density is written as a reference density r 0 plus a small fluctuation dr, (15) This electronic density is then inserted in equation 14: (16) where r' 0 = r 0 (r  ') and dr' = dr(r  ') are defined as shorthand notations.The second term in equation 16 corrects the double counting in the Coulomb term; the third term corrects the new exchange-correlation contribution; and the fourth term results from splitting the Coulomb energy into one part related to r 0 and another related to dr.E nn is the nuclear repulsion.
Afterwards, E xc [r 0 + dr] is expanded in a Taylor series up to the second-order term: The sum in the first line of equation 18 is analogue to E bnd in equation 1.The terms in the second line of equation 18 define the repulsive contribution, (20) Finally, the last term in equation 18 includes the corrections related to the fluctuations in the electronic density.This term is defined as ( 21) Therefore, equation 18 can be rewritten as (22) In order to obtain a good estimate of the reference electronic density, r 0 is written as a superposition of atomlike densities centered on the nuclei a, (23) With this approximation it is assured that E rep does not depend on the electronic-density fluctuations.Furthermore, due to the neutrality of r a 0 the Coulomb contributions become negligible for long distances.Therefore, E rep can be expanded as (24) The contributions of 3 and more centers are rather small and can be neglected.These approximations can also be justified by Coulomb screening, i.e., since r a 0 is the electronic density of a neutral atom, the electron-electron interaction terms with more than two centers are canceled by the nucleus-nucleus interactions.
Due to the screening of terms of more than two centers, one can assume the two-center contributions to be short ranged.However, the repulsion energy does not decay to zero for long interatomic distances.Instead, it decays to a constant value given by the atomic contributions: (25) Thus, S N a E rep [r a 0 ] is assumed in order to make E rep dependent only on two-center contributions: (26) Although it would be possible to calculate E rep for known values of r a 0 , it is more convenient to adjust E rep to ab initio results.Thus, E rep is fitted to the difference between the DFT energy and E bnd as a function of the interatomic distance R ab using a suitable reference structure, i.e. (27) The value of E bnd can be obtained by diagonalization of the Hamiltonian matrix, which leads to (28) The value of E rep is usually fitted to a polynomial function or to a series of splines.Typical plots of E DFT , E bnd and E rep for a reference structure are shown in Figure 2.
Based on the considerations discussed so far, the DFTB model can be derived.

The Standard DFTB Model without Self-Consistency
In the standard DFTB scheme, the second-order correction term, E 2nd of equation 22, is neglected.Therefore, the calculation of the total energy does not depend on the electronic-density fluctuations dr and, accordingly, it does not have to be solved iteratively.
In DFTB the KS orbitals are represented with a linear combination of atomic orbitals (LCAO) centered on the nuclei.Denoting the basis functions by f n and the expansion coefficients by C in one can write the KS orbitals in the form (29) From this LCAO model, one obtains the secular problem (30) where the elements H 0 mn of the Hamiltonian matrix and S mn of the overlap matrix are defined as follows: (31) The second term of equation 22 can be transformed, with equations 29 and 11, into  (32) in which the elements of the density matrix P are defined as follows (33) In order to restrict the LCAO to valence orbitals only, it is necessary to assure the orthogonality of the basis functions with respect to the core basis-functions of the remaining atoms (by using atomic orbitals as basis functions the orthogonality between the core and valence functions within the same atoms is already assured).
Denoting |f) as a non-orthogonalized basis-function and |f b c ) as the core basis-functions of atom b, the corresponding orthogonalized basis-function of |f) is obtained by: (34) By using this orthogonalization procedure, equation 32 is transformed into (35) where e b c denotes the eigenvalue of the state c in atom b.The effective potential n KS and the core correction in equation 35 can be interpreted as a pseudo-potential (V pp ).Writing n KS as the sum of potentials V a centered on the atoms, (36) and using this definition in equation 35, the effective potential is transformed into a pseudo-potential for all atoms in the system, except for atoms to which f m and f n belong.Therefore, the pseudo-potential appears in the three-center terms and in the two-center terms whose valence orbitals belong to the same atom (so called crystal field terms).The pseudo-potential contributions are considerably smaller than the contributions of the full potentials and are neglected.Thus, the Hamiltonian matrix elements are defined as (37) where d ab is the Kronecker's delta.This approach, the potential superposition, has been used since the 1980's for the calculation of DFTB parameters.In 1998, Elstner et al. 36   presented an alternative approach to derive the DFTB equations through a second order expansion of the DFT total energy with respect to the electron density.As result the Hamiltonian matrix elements are calculated as density superpositions, which is identical to equation 37 except for the contribution of the exchange correlation potential.Indeed, due to the non-linear nature of n xc , the effective potential cannot be described as a simple sum of reference potentials within this approach, instead one obtains (38) Both approaches are physically motivated and their results are similar, which is not surprising if the potential difference between equations 37 and 38 is explicitly calculated.Both approaches have been used extensively in the past, the potential superposition being more popular for standard DFTB calculations, and the density superposition more widely used for SCC-DFTB.
The f n basis functions and the reference atom-like densities r a 0 are obtained by solving the Schrödinger equation (39) for the free atom within a self-consistent DFT method, as shown in Figure 1.The contraction potential (r/r 0 ) 2 in equation 39 constrains the wave functions, resulting in better basis sets for the study of condensed-phase systems and free molecules as well.The value for the parameter r 0 is normally chosen between 1.85r cov and 2r cov , with r cov being the atomic covalent radius. 37 In practice, the Hamiltonian matrix elements are calculated as follows: For the diagonal elements the energy level of the free atom is chosen, which ensures correct dissociation limits.Due to the orthogonality of the basis functions the offdiagonal elements of the intra-atomic blocks are exactly zero.The interatomic blocks are computed as given in equation 37 or 38, depending on the choice of potential generation.Within the density superposition approach the Hamiltonian matrix elements unfold as follows: (40) It should be noted that the Hamiltonian elements H 0 mn depend only on atoms a and b and, therefore, only the twocenter matrix elements are explicitly calculated, as well as two-center elements of the overlap matrix.According to Vol. 20, No. 7, 2009 equation 40 the free atom eigenvalues form the diagonal of the Hamiltonian matrix, which assures the correct limit for free atoms.
By using f n and r a 0 the Hamiltonian and overlap matrix elements can be calculated and tabulated as a function of the distance between atomic pairs.Thus, it is not necessary to recalculate any integrals during, e.g., a geometry optimization or molecular dynamics simulation.
At last, an analytic expression for atomic forces can be derived from the total energy with respect to the atomic space-coordinates, (41) By this approach, the DFTB method covers all three requirements for an atomistic tight-binding approach.

The Self-Consistent Charge Correction: SCC-DFTB
The non-self-consistent DFT scheme described so far is very suitable to study systems in which the polyatomic electronic density can be well represented by a sum of atomlike densities, i.e. homonuclear covalent systems or highly ionic systems.However, the uncertainties in the standard DFTB increase when the chemical bonds in the system are controlled by a more delicate charge balance between atoms, especially in the case of heteronuclear molecules and polar semiconductors.In order to have a better description of electronic systems and better transferability of DFTB in the cases where long-range Coulomb interactions are significant, the method has been improved, giving rise to the self-consistent charge correction DFTB (SCC-DFTB). 36In this new scheme, the electronic density is corrected through inclusion of the second-order contributions E 2nd in equation 22, which are neglected in standard DFTB.
In order to include the density fluctuations in a simple yet efficient way according to a tight-binding approach, dr is written as the superposition of atom-like contributions dr a , which fast decay along the distance from the corresponding atomic center, (42) where the atom-like contributions can be simplified with the monopole approximation: (43) Here Dq a is the Mulliken charge, difference between the atomic Mulliken population q a 38 and the number of valence electrons of the neutral free atom q 0 a (Dq a = q a -q 0 a ); F a 00 denotes the normalized radial dependence of the density fluctuation in atom a approximated to spherical by the angular function 00 .In other words, the effects of charge transfer are included, but changes in the shape of the electronic density are neglected.Equation 21then becomes (44) in which the notation g ab was introduced merely for convenience.
In order to solve equation 44, g ab must be analyzed.In the limit case where the interatomic separation is very large one finds, by GGA-DFT, that the exchange-correlation term goes to zero and g ab describes the interaction of two normalized spherical electronic densities, basically reducing to 1/|R In the opposite case, for which the interatomic distance tends to zero (|R , g ab describes the electron-electron interaction within the atom a and can be related with the chemical hardness h a , 39 or Hubbard parameter g aa = 2h a = U a .Typically, the atomic hardness can be calculated using the difference between ionization potential I a and electron affinity A a of atom a: 2h a = I a -A a .Due to practical problems, in particular related to the non-existence of various anions and accordingly missing experimental validation of the electron affinity of the corresponding elements, it is more convenient to exploit DFT to obtain these parameters.Application of Janak's theorem 40 relates the atomic hardness to the derivative of the HOMO energy with respect to the occupation number of the HOMO and hence the energy change with respect to electron change within the HOMO.This approach offers the possibility to treat the charge contribution shell-or even orbital-wise, which is important for the calculation of certain elements with sp and d bonding contributions, in particular for transition metals.
Orbital hardness values have been reported in the literature for elements from H to Xe. 41 In the following, we concentrate on the atomic SCC procedure, which implies that all sums over charges run over the atomic index a.For orbital-dependent SCC the summation index for the charge would run over the shell index x.Within the monopole approximation, U a can be calculated, using a DFT procedure, as the second derivative of the total atomic energy of atom a with respect to its atomic charge: In order to obtain a well-defined and useful expression for systems in all scales, and still keep consistence with the afore approximations, an analytical expression was developed 36 to approximate the density fluctuations with spherical electronic densities.In accordance with Slatertype orbitals (Gaussian-type orbitals can also be employed) used to solve the KS equations, 42,43 it is assumed an exponential decay of the normalized spherical electronic density: (47) Omitting the second-order contributions of E xc in equation 44 one obtains: (48) Integration over r  ' gives: where s is a short-range function with exponential decay, so that (51) Once it was assumed that the second-order contribution can be approximated by the Hubbard parameter when R = 0, according to equation 46, the exponents of equation 51 are obtained: (52) This result can be interpreted by noting that harder elements tend to have localized wave functions.The chemical hardness of a spin-depolarized atom is calculated by the energy derivative of the highest occupied atomic orbital with respect to its occupation number, equation 46, using a fully self-consistent ab initio method.Therefore, the influence of second-order contributions of the exchange-correlation energy is included in g ab for short distances, where it is important.The fact that, within GGA, the exchangecorrelation energy vanishes for large interatomic distances is taken into account.In the case of periodic systems, the long-range part can be calculated using the standard Ewald summation, whereas the short-range part s decays exponentially and can be summed over a small number of unit cells.Thus, equation 50 is a well-defined expression for extended and periodic systems.
Finally, the total energy within SCC-DFTB is written as (53) with Here the contribution due to the Hamiltonian H ^0 is exactly the same as in the standard DFTB scheme.Note that the first term in equation 53 does only simplify to the sum of MO energies, the convenient notation for DFTB, if all charges are zero.Like in the nonself-consistent method, the wave functions y i are expanded in a LCAO model, equation 29, and equation 53 gives: The charge fluctuations are calculated by Mulliken population analysis: 38 (55) and secular equations similar to those in equation 30 can be obtained, with modified elements in the Hamiltonian matrix: (56) The matrix elements H 0 mn and S mn are identical to those defined in the standard DFTB method, in equation 31.
Since the atomic charges depend on the monoatomic wave functions y i it is necessary to use a self-consistent procedure.Once the elements S mn extend to some neighboring atoms, multi-particle interactions are introduced.The second-order correction is achieved by introducing the elements H 1 mn , which depend on the Mulliken charges.
Identically to the standard DFTB, the repulsive potential is fitted according to equation 27 using a suitable reference system.
As the self-consistent charge correction allows for the explicit treatment of charge-transfer effects, the transferability of E rep is considerably better, in comparison with the non-self-consistent scheme.
As in the standard DFTB, a simple analytic expression for the atomic forces can be derived accordingly: (57) DFTB schemes have been successfully used in a wide range of applications, from molecular compounds 22,44 to systems in solid state. 19,45-47Indeed, a symposium dedicated to the DFTB methods was held during the 232 nd National Meeting of the American Chemical Society, from 10 th to 14 th of September, in 2006.A special section with contributions presented in this symposium was published in the Journal of Physical Chemistry A, issue 26 of 2007, 48 presenting the actual development state of DFTB with respect to its formalism, implementation and applications.

Weak Forces: Dispersion-Corrected (SCC-) DFTB
London interactions, also called dispersion forces, are defined as attractive forces between nonpolar molecules, due to their mutual polarizability. 49London dispersion forces are several orders of magnitude weaker than typical covalent or ionic interactions and also about 10 times weaker than hydrogen bridge interactions.Therefore, dispersion forces have negligible effect in short-range interactions and can be understood as the long-range component of van der Waals forces.
Despite their weak nature, London interactions affect many fundamental processes in chemistry, physics, and biology.They influence the formation of molecular crystals, the structure of biological molecules such as proteins and DNA, adsorption processes, π-π stacking interactions, among others.
However, as explained above, both the standard and self-consistent DFTB methods treat only short-range atomic potentials and terms with more than two centers are neglected.Therefore, the Hamiltonian matrix elements fall off quickly and become negligible at interatomic distances typically found in the region of the van der Waals minimum.Hence, DFTB completely disregards van der Waals interactions, especially dispersion forces.
Two treatments meant to include dispersion interactions a posteriori have been proposed. 50,51In both cases the dispersion energy E disp is calculated separately using empirical potentials and then added to the DFTB total energy expression.Since van der Waals forces are totally absent in DFTB, the addition of E disp does not introduce any double-counting errors to the energy.
Since both treatments are somewhat similar, we describe that used in the present work. 51This correction was implemented in an experimental version of the deMon code 52 and makes use of the UFF force field, 53 already available in deMon.The dispersion interaction U ab between atoms a and b at a distance R is given in Lennard-Jonestype form, which includes two parameters: van der Waals distance (R ab ) and well depth (d ab ): (58) The R ab and d ab parameters are reported in the original paper 53 and are available from H to Lw in the periodic table of elements.In UFF the van der Waals term is set to zero according to an adjacency criteria; however, this imposes an inflexible topology of the system, which is not desirable in a quantum-mechanical method.To overcome this problem, equation 58 is used only when U ab is attractive (London interactions are never repulsive), i.e.R < 2 -1/6 R ab .
In addition, a short-range potential is derived using the polynomial (59) where U 0 , U 1 , and U 2 are calculated to make the interaction energy and its first and second derivatives match equation 58 at R = 2 -1/6 R ab .The best value suggested for n is 5, which gives the following U 0 , U 1 , and U 2 parameters: 51   (60) (61) (62) Therefore, the dispersion potential for the DFTB method can be written as (63) and the dispersion energy is given by (64) This term is then added to the total DFTB energy calculated either using standard DFTB (section 5) or the SCC scheme (section 6).

Glycine in Aqueous Solution
Glycine (aminoethanoic acid) is the simplest among the essential aminoacids.In solution an intramolecular proton transfer from the carboxylic group to the amino group takes place, establishing the zwitterionic equilibrium shown in Figure 3.The charge separation in the zwitterionic form is stabilized by the solvent, which must have large dielectric constant, as it is the case in water.Thus, the neutral species is favored in nonpolar solvents.
In this work, Born-Oppenheimer molecular dynamics was carried out using the DC-SCC-DFTB method, as implemented in the deMon package. 52The glycine molecule was placed within a 16 Å periodic box containing 129 water molecules.
The data were collected during a 100 ps simulation time with a 0.5 fs time step.The simulation was carried out after a thermalization time of 50 ps.It is important to emphasize that both glycine and water molecules were calculated within a full quantum-mechanical approach.The radial distribution function (RDF) of water with respect to the glycine center of mass is shown in Figure 4.The first solvation shell integrates to 22 water molecules.
Table 1 shows the calculated geometrical properties of glycine in solution.The optimized geometric parameters are shown at the PBE/TZVP and DC-SCC-DFTB levels of theory.The estimated angles are in good agreement with the previously published results. 54The O-C=O angle presents the largest discrepancy for the zwitterionic form.The PBE/ TZVP estimated O-C=O angle is 13 degrees larger than the value estimated with DC-SCC-DFTB.Furthermore, the O-C=O angle is expected to increase from the neutral to the zwitterionic form due to the deprotonation of the carboxyl group.However, DC-SCC-DFTB seems to be insensitive to the large charge on the deprotonated carboxyl and the angle remains similar to that in the neutral form.The mean values of the angles and dihedrals from MD (last column in Table 1) are close to the optimized values with a standard deviation of about 4 degrees, except for the O-C-C-N dihedral.This dihedral involves rotation around a single C-C bond, therefore, large standard deviation is indeed expected, explaining the apparent disagreement with the gas phase PBE/TZVP results.
Wada et al. 55 estimated the Gibbs free energy variation between the two glycine forms in aqueous solution to be about -7.0 kcal mol -1 (DH = -10.3kcal mol -1 ).The change of the expected value between the two forms from the NVE molecular dynamics (DE NVE ) was estimated to be about -25.5 kcal mol -1 .We have also used the continuum model to estimate the DG of this reaction at the PBE/TZVP/ PCM level of theory and a value of -23.4 kcal mol -1 was obtained.

Final Remarks
DFTB is an approximate density-functional method which, in principle, does not employ any empirical parameter, in the sense that all quantities are calculated within DFT (Slater-Koster integrals) or they are calculated from reference structures by DFT calculations (E rep ).It has been implemented in many different codes. 56  Density Functional methods have along the time become a standard method for electronic structure calculations and substantially helped to unify organic chemistry, inorganic chemistry, surface chemistry, materials science and, more recently, biochemistry. 4With the advent of DFTB, the approximate DFT method, a plethora of challenging systems are now accessible for electronic structure calculations, enlarging the frontiers of the applicability of fundamentally well established theoretical tools.Nanostructured, self-assembled and nanoreactor systems are some of those for which DFTB can provide substantial help in the investigative work.

( 17 )
Substitution of equation 17 into 16 and use of the definition (d E xc /dr) r0 = n xc [r 0 ] results in (18) From equation 18 it is possible to define four important terms.The first is a reference Hamiltonian H ^0 depending only upon r 0 (19)

Figure 1 .
Figure 1.Flow-chart of a typical DFT calculation within the Kohn-Sham method.

Figure 2 .
Figure 2. Typical plots of E DFT , E bnd , E rep and total DFTB energy against the interatomic distance for a reference structure.

Figure 4 .
Figure 4. Radial distribution function of water around the glycine center of mass.