Complete free energy landscape and statistical thermodynamics of single poly(ethylene glycol) molecules

In this paper, we revisit the statistical mechanics of poly(ethylene glycol) with the intention to check ad hoc approximations made in previous works, namely to keep all C–O bonds along the chain in the trans configuration. The theory developed here includes all the bonds on an equal footing. Starting with density functional theory calculations for the ethylene glycol monomer, dimer and trimer one can determine all the energetic and geometric parameters needed for a rotationally isomeric state model with interactions between nearest neighbour monomers included. The properties of short molecules are presented in detail, including the force-extension curves, end-to-end distributions, and also a discussion of the local helix properties. Finally, we use the transfer matrix method for long chains of 200 ethylene glycol monomers, i.e. 600 bonds. The agreement with the experimental force-extension curve is excellent.


Introduction
Poly(ethylene glycol) (PEG) is a simple polymer, yet complex enough to show a number of interesting features revealed, e.g. in atomic force microscopy (AFM) experiments. It consists of a sequence of ethylene glycol subunits (-O-CH 2 -CH 2 -, denoted by (EG)), with either methoxy or hydroxy terminal groups. Solvated in hexadecane, its force-extension curve shows the simple features of a freely joined chain (FJC) response, i.e. it is dominated solely by entropic effects, except at the highest extensions, where spring-like behaviour is observed [1,2]. On the other hand, when water is used as the solvent, a plateau-like hump appears that is reminiscent of two-state behaviour usually reserved for more complex molecules such as polysaccharides.
The first models concerning the statistics of such simple polymers in general were aimed at simply accounting for a one-dimensional connected structure as well as the volume exclusion effect, i.e. making sure that atoms do not overlap. To some general degree, the mean-field calculations of Edwards [3] or the renormalization group approach by de Gennes [4] are applicable to all macromolecules in solution. The phenomenological Flory-Huggings model [5] is based on calculating the free energy of a molecule in solution by accounting for the free energy of mixing as well as a molecule-solvent interaction, generally specified by a single parameter χ. An extension of the latter theory using a continuous mean-field theory [6,7] as well as a lattice model [8,9] made it possible to show that the PEG phase diagram can be explained by assuming clusters of molecules, held together by hydrogen bonds. This property has been corroborated by Smith and co-workers [10]- [14], who use exact calculations on dimethoxyethane ((EG) 1 ) for the development of molecular force fields, which in return are used in molecular dynamics calculations.
In this work, we outline a theory that starts from first principles calculations and systematically introduces approximations to approach larger scales in the modelling of single PEG molecules. Thus, an explanation of the macroscopic response and properties of a long polymer chain is ultimately given in terms of the structural and energetic properties of the (EG) monomer subunits [2]. The starting point for an appropriate statistical theory should be the energy spectrum and the density of states, which can be obtained from quantum mechanical electronic structure calculations. The latter can only be done for small systems because the computational effort increases at least as the fourth power of the number of electrons in the system. This paper 3 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT aims to develop a strategy that leads to an exact statistical theory of a single PEG molecule, potentially interacting with a solvent. To treat such a long polymer molecules the following approach has been developed.
In a first step one calculates the electronic properties of a short polymer chain from first principles by, for instance, employing density functional theory (DFT) or Hartree-Fock theory complemented by Møller-Plesset perturbation theory. The result of such calculations is the energy spectrum and the density of states (the number of independent wavefunctions or conformers at a given energy) for a short chain [15]. This provides all the input required for the statistical mechanics of short chain molecules.
Although experiments are generally done with long polymers, obtaining the force-extension curve for short chains is still a useful exercise. It provides a quantitative description of the late stages of the elastic response of long chains, where the molecule is more or less fully stretched already so that entropic contributions are less important. The total response of a long chain is more or less additive over the individual monomers, provided the first principles calculations are done as accurately as possible for all the conformers of the (small) molecule.
First principles calculations were reported for (EG) n with n = 2, 3, 4 and its interaction with water, using DFT with nonlocal gradient corrections [15]. At that time, a tractable number of conformers was selected, chosen to be representative of those states accessible at room temperature. We recall that the conformers of a macromolecule form local minima on its total energy surface. As such, a conformer is a static concept and its relevance for a system at nonzero temperatures hinges on timescales: if the activation barrier for escape from the local minimum is so high that the time spent by the molecule in this local minimum is much larger than the time to escape from it, we can make definitive measurements on this state and assign it the status of a conformer. It was argued in the previous study that most low energy conformers have the C-O bonds in the trans configuration and the C-C bonds in trans and clockwise and anticlockwise gauche configurations, resulting in three states of the EG subunit, (t|g + |t), (t|g − |t) and (t|t|t), to be placed on three positions along the chain, say (EG) 3 , for a total of 3 3 = 27 combinations or conformers.
Since the first principles approach just outlined can only be done for short chains of a few monomers, an extension of the theory to deal with long chains must invoke some approximations, by necessity. The latter are based on selecting from the large number of quantum states of a polymer, or conformers, those that are energetically accessible at the temperature of interest and dynamically relevant. This suggests multistate interacting chain models (ICMs) of which the three-state rotationally isomeric state (RIS) model is an early version. Such models have a number of parameters such as self-energies of the various states, interaction energies between nearest and next-nearest neighbours along the chain and geometric parameters like bond lengths, bond angles and dihedral angles. The result is an ICM with all parameters determined from first principles quantum mechanical calculations, valid for polymers of arbitrary length [16]. In recent years, we have developed a Green function or transfer matrix (TM) method which allows one to solve the complete and exact statistical mechanics of ICMs in both the Gibbs and Helmholtz ensembles [16]- [19].
The restriction to calculating only a few conformers for short chains was ultimately proved to be a good (and economical) strategy in that the multistate RIS model gave perfect agreement with experiment. Not only that, this theory also explains the new two-state feature when stretching is done in an aqueous environment [1,2,16]. Still, the question remains whether a more complete 4 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT theory will give further information. This paper is dedicated to this endeavour. Our strategy is as follows.
1. Using DFT, we calculate the geometry and energies of all 3 3 = 27 conformers of (EG) 1 , 3 6 = 729 conformers of (EG) 2 and also for a random selection of 493 (out of 3 9 = 19 683) conformers for (EG) 3 . Within the framework of quantum mechanics, a conformer is a local energy minimum in the configurational space spanned by the coordinates of the nuclear coordinates of all the atoms in these molecules. The results of such calculations are the respective densities of states, D N (E, L)dEdL, i.e. the number of conformers in intervals dEdL with energy E and length L. From this we can obtain the Helmholtz and Gibbs partition functions and free energies from which the mechanical equation of state, i.e. force-extension curve, and the thermal equation of state can be calculated by simple derivatives. The reason why we look concurrently at the Helmholtz and Gibbs ensembles is because the former lends itself to the calculation of force fluctuations and the latter to length fluctuations, and both are accessible in AFM experiments. 2. With the density of states we can define the free energy landscape which becomes apparent when we rewrite the partition functions as and similarly for the Helmholtz partition function. 3. Next we extract energetic and geometric parameters for the three-state RIS model including interactions between adjacent monomers. From this we recalculate the force-extension curve and find, not surprisingly, excellent agreement with the exact calculation. 4. Going beyond the shortest chains we use the RIS model to calculate the force-extension curve for (EG) 7 for which we also use the TM method (with the same geometric and energetic parameters), again finding excellent agreement.

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
5. The TM method can be used to calculate the partition function for long chains which we do for N = 200. The resulting force-extension curve we then compare with that obtained in previous work based on using only 10 non-redundant conformers of (EG) 3 as input. The agreement obtained is the ultimate justification of that simpler model as agreement with experiment is also guaranteed.
The paper is structured as follows: in the following section, we give the necessary details of the DFT calculations and results for short chains. We then construct the Hamiltonian for the RIS model, outline a method for the exact evaluation of its partition function and give details of the transfer TM. In section 3, we present force-extension curves for medium-length and long chains calculated in these frameworks. The last section contains a summary, conclusions and an outlook to future work. Some numerical details are given in appendices.

First principles calculations
Our calculations of all the conformers of (EG) n with n = 1, 2 and 3 employ DFT with gradient corrections and have been done with the Gaussian 03 suite of programs [20]. The exchangecorrelation potential used is BP86 and the basis set is 6-311++G**. A justification of this choice and a thorough comparison with earlier work has been given elsewhere [15]. We then adjust all these parameters and allow rotation around the C-C and/or C-O bonds to find a local minimum in the total electronic energy. That this is possible for a conformer, as defined above, is not surprising when one considers that in the 84-dimensional space of the coordinates of all 30 nuclei of (EG) 3 (3 × 30 − 6 for the centre of mass translation and rotation of the molecule) these minima are far from each other with activation barriers in between. Activation barriers do not enter the calculation of equilibrium properties but of course are important to estimate the time scale of the approach of the system to equilibrium as we indicate above. We also point out that, in the various conformers, bond lengths, bond angles, and dihedral angles are all slightly different from their standard values.
Symmetry considerations reduce the number of conformers to 10 for (EG) 1 , 196 conformers of (EG) 2 and 5002 for (EG) 3 . Convergence of the self-consistent DFT iterations can be reached for nine conformers of the monomer, only 135 for the dimer and 493 of 1040 randomly selected conformers of the trimer. Lack of convergence is interpreted as evidence for volume exclusion, because it occurs when a conformer contains the sequence g+g−g+ so that the monomer doubles up on itself. Moreover, all conformers that the RIS model predicts to be excluded did not converge either.
For the construction of the RIS model we need to extract average bond lengths, bond angles and dihedral angles from all conformers, which can be done by straightforward geometric averaging. However, this procedure would overestimate the contributions from the highly deformed, high-energy conformers. We therefore also use a weighted averaging with a Boltzmann factor involving the energy of a given conformer. Table 1. Geometry parameters obtained by averaging over (EG) 1 , (EG) 2 , and (EG) 3 : bond distances b, bond angles γ, and dihedral angles ϕ. The top number in each field is the Boltzmann-weighed average (used in analysis), while the bottom value in brackets denotes the arithmetic average. Table 1 presents the resulting bond lengths, bond angles and the energies relative to the all-trans conformer, for the monomer. For the dimer and trimer we give plots of all converged conformers in the energy verses length plane in figures 1 and 2. These plots represent the density of states D(E, L) which can be turned into a continuous function by appropriate binning of the discrete dataset. Important from the viewpoint of statistical mechanics is the band of conformers below about 50 meV or 2k B T at room temperature, well separated from the higher energy conformers.

Exact calculations in the RIS formalism
We use a variation of the RIS formalism to describe the thermodynamics of a long chain. A complete description of a molecule with a given chain length necessitates a summation over all possible conformers. Each conformer with N monomers is uniquely described by a set of vectors { n 1 , n 2 , . . . , n N } describing the rotational states s i ∈{g+, t, g−} of each atomic bond (of which there are three for the EG monomer) as A one-dimensional lattice gas Hamiltonian can be used to fully determine the energetics of a chain of N monomers, see [16] for the details with the (t| * |t)-approximation.
The nearest-neighbour potential V NN refers to the interaction of complete monomers and thus contains interactions including up to five neighbouring bonds along the backbone. This analysis allows for an exact representation of the energetics of (EG) 2 . As for the conformers of (EG) 1  196 of which are unique. Sixty-one of these interaction energies are also set to infinity due to volume exclusion effects. The calculation of these parameters is described in appendix A. Our Hamiltonian (5) corresponds to a first step from pure quantum mechanics toward a statistical model and still contains a maximally possible amount of microscopic information. The other ingredient for the RIS model is the geometry, which is based on the averages given in table 1. Given the sequence of dihedral states n, all bond vectors of a given conformer can easily be calculated with the successive application of rotation matrices. Defining the first bond to be along theẑ-axis of our coordinate system, the nth bond vector b n is in the direction where the operator R z (ϕ) rotates the coordinate system by an angle ϕ about theẑ-axis and R y (γ) rotates it by γ about theŷ-axis. From equation (6) one can calculate all required geometric properties of one conformer. As an example we give the end-to-end vector r N for a conformer 9 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT consisting of N b = 3N bonds (3 for each EG monomer), The location of the hydrogen atoms are obtained in a similar fashion. Since all of the atomic positions are known in this approach, it is a simple task to check whether or not a given conformer of (EG) n is allowed. A conformer is rejected due to excluded volume interactions if any two atoms in the molecule overlap, as defined by their respective van der Waals radii [21,22]. All exact RIS results presented in this work are shown twice, contrasting calculations that account for excluded volume and calculations that 'neglect' these effects.
Next comes the evaluation of the partition functions, equations (1), from which equilibrium thermodynamics follows straightforwardly. The computational effort for this calculation can be reduced by reducing the number of redundant conformers, which can be done with simple operations on the vector n defined in equation (4). For our purposes, a conformer is called redundant when its sequence vector n can be obtained from that of second conformer n either by inverting the sequence, by exchanging the gauche-states g+ ↔ g−, or by using a combination of both. The resulting possible degeneracies for unique conformers are g n = 1, 2, or 4.
According to equations (1), there is a choice between Helmholtz and Gibbs ensembles for a fixed number of bonds. Given that the overall number of conformers for (EG) n is still of the order of O(3 3n ), it is advisable to do all initial numerical analysis in the Helmholtz ensemble, where a complete description of the system can be obtained with each conformer taken into account exactly once. Due to the sheer number of available conformers, it is numerically sensible to bin the results of similar lengths over some small interval Using the same formalism, it is possible to write down the end-to-end distance probability distribution which describes the probability of finding the length of a molecule to be equal to L and is used to calculate the mean square length of the molecule, R 2 = dLL 2 P r (L).

TM method
The all-inclusive RIS calculations described above, are only possible for short chains with a limited number of conformers; our longest chain was (EG) 7 . An additional approximation enables one to use the TM method for which the limit on the number of (EG)-units is at about 200 [16]. This method was developed primarily to calculate force-extension curves and various other properties of single molecules [16,18] , but it has been shown that the TM method is a powerful tool to use for modelling single molecules in external potentials or even in a self-assembled monolayer [19,23,24].
In the TM method, the state of the nth bond in a molecule is described via a Green function G n ( ) [3,5]. The Green function for the following monomer can be obtained by defining a transfer operator T ( , ) and integrating the Chapman-Kolmogoroff equation, G n+1 ( ) = d T ( , )G n ( ). A given bond in a molecule is completely described by specifying a set of coordinates , whose choice depends on the problem at hand and the underlying polymer model used.
In the particular case of the RIS model, must include the spatial location and orientation of the nth bond as well as its rotational state (g−, t, or g+). The spatial symmetry of the problem allows the use of a single spatial coordinate, which is chosen to be the z-coordinate in our coordinate system for computational reasons. Starting with the orientation of the nth bond, specified by θ n and z n in the usual cylindrical polar coordinates, the orientation of bond b n+1 is calculated with the same rotation matrices used in the RIS geometry, equation (6) b n+1 (s n+1 ) = (b n+1 /b n )R z (ϕ n+1 )R y (γ n+1 )(b n cos θ n , b n sin θ, z n ) T .
The dependence on the rotational state s n+1 enters through the choice of the dihedral angle φ n+1 as per table 1. In the numerical implementation it proves advantageous to do the actual TM calculation by using z-projections of bonds n and n + 1 as independent variables, rather than z n and the angle θ n . This effectively saves one (out of four) integrations in the calculation of the Green function for each bond. The independent variables for the transfer operator are i = {z i , z i , z i−1 , s i } and the operator itself has the form The specification of proper initial conditions for this procedure involves an artificial bond with index zero that precedes the first real bond. This bond has to be taken into account in such a way that the Green function of the first real bond contains all possible orientations and rotational states with equal probability. Finally we should note that an appropriate discretization of the variables is necessary to actually calculate properties of a large molecule. The Green function then becomes a vector containing some 10 8 elements and the transfer operator (11) takes on the form of a square matrix of this dimension. Each integration of the Chapman-Kolmogoroff equation then corresponds to a matrix multiplication.
The evaluation of statistical equilibrium properties, such as the force-extension curve from the formalism outlined in this section, follows from the fact that the Green function for the last bond in effect corresponds to the Helmholtz partition function when it is integrated over all orientation variables [16], N (z N , z N , z N−1 , s N ). (12) The TM also requires microscopic geometry and energy information for the proper evaluation of equation (11). The bond lengths, bond angles, and dihedral angles can be taken directly from the DFT calculation as was done in the RIS model. However, the Hamiltonian 11 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT underlying this model is somewhat different from the RIS Hamiltonian (5). In the TM case, everything is based on single bonds and their first neighbour, rather than on monomers. Thus the Hamiltonian H TM has the form The interactions V NN only take into account the immediately neighbouring bonds, in contrast to the whole monomer that was part of the nearest neighbour interaction V NN in the RIS case. The nature of the energies E self and V NN does not allow for a direct calculation from DFT data. However, it can be obtained as a degeneracy-weighed best fit to all 493 available (EG) 3 conformers using a Hamiltonian (5) with interactions up to nearest-neighbour dihedrals. This procedure still accounts for the exact geometry up to four bonds along the chain. The results are presented in appendix A, along with the RIS energy parameters.
Although the TM in its present implementation accounts fully for interactions up to four neighbouring bonds along the chain, it does not include volume exclusion between more distant atoms as is done in the exact calculations for short chains. This neglect, however, should be of little consequence when the TM method is used to calculate the stretching response for very long molecule: when a force is applied, the relevant extended configurations are less susceptible to the excluded volume interactions. This will be shown explicitly with examples below.

Results
In this section, we will present the results for the thermodynamic properties (such as forceextension curves) of short and long chain molecules based on the various analytic tools that we developed in earlier sections. We begin with short chains. In the bottom panel of figure 1 we show the force extension curves for (EG) 2 calculated at various levels of approximation. The solid line results from the inclusion of all conformers obtained from the converged DFT calculations. For the dash-dotted line we only included those conformers with energies less than 50 meV or 2k B T (at room temperature). As this excludes some of the high-energy short conformers the curve goes into the compressive mode a little earlier. Lastly, the dashed line results when we only account for the conformers where the CO bonds are kept in the trans configuration, i.e. conformers of the form (t * t − t * t) where * stands for trans or the two gauche configurations of the C-C bonds. This excludes almost all the short conformers and the compressional regime starts very early. The latter approximation was used in our previous work [16]. This clearly shows that for (EG) 2 this is not a good approximation. Note however, that in the regime of high stretching all curves agree because there only the longer trans configurations are important. The same conclusions are reached for (EG) 3 as shown in the bottom panel of figure 2.
Moving on to longer chains we look in detail at (EG) 7 . This is the longest chain for which we can, in a reasonable time, do exact calculations with excluded volume interactions. We start in figure 3 with a plot of the density of states obtained from the RIS model showing the results without and with volume exclusion. Volume exclusion reduces the number of viable conformers, in particular for higher energies. This is simply due to the fact that with excluded volume conformers with multiple gauche defects are suppressed. Of interest here are the bands with a high density of states at fairly low energies, i.e. around 50 meV. Those are, of course, the conformers that will dominate the thermodynamics of a single PEG molecule. Clearly the probability distributions for the end-to-end length with and without volume exclusion must differ with the former peaking at higher lengths, as seen in figure 4. The effect of volume exclusion is almost negligible in the force-extension curve. This is due to the fact that in the low-force regime the curve is dominated by entropic effects for which the exclusion of a few conformers is not crucial. On the other hand, in the high-force regime the chain is more or less fully stretched so that again volume exclusion does not play a role. We have also employed the TM method to calculate the force-extension curve and find that with all conformers included the result is almost identical to the curve calculated with the exact calculation apart from small numerical discrepancies. On the other hand, including only the trans configurations of the C-O bonds, i.e. only EG monomers of the form (t * t − t * t), the polymer is softer. This is due to two facts: with fewer conformers entropy is lower so that less force is needed to unfold the polymer; and these conformers are among the longest so that the unfolding immediately results in longer conformers.
Next, we discuss the experimental claim that local helical structures are an important contribution to the statistics of single PEG molecules. A local helix of n monomers can be defined as a sequence (tgt) n of consecutive gauche bonds (all of the same helicity) anywhere in the backbone.As an example, we choose the 2/7 helix that was used as one of the two fundamental conformers to model the force response in previous studies [1], based on the argument that this type of helix is observed in crystalline PEG. Within the RIS formalism, it is an easy task to sort out these local helical contributions to the radial and other distribution functions. Figure 5 shows the details of the radial distribution functions with at least n = {2, 3, 4} consecutive helical monomers. These curves demonstrate that even sequences with n = 4 are highly improbable and have a negligible contribution to P(R) and hence for all thermodynamic properties.  More data is available to support this claim. The relative probability of finding two neighbouring bonds among all possible rotational states has been calculated directly from averaging over all conformers in (EG) 7 . For the O-C-C and C-C-O bond pairs we find P gt = 0.531, P gg = 0.149, P tt = 0.232, P tg = 0.040, and P gg = 0.048, while the results for the C-O-C bond pairs are P tt = 0.585, P tg = P gt = 0.376, P gg = 4 × 10 −4 and P gg = 0.038. From the above probabilities it can be estimated that the probabilities of finding one, two, or three consecutive (EG) units in the same helical state are 7 × 10 −2 , 3 × 10 −3 and 2 × 10 −4 respectively. Smith et al [11] have measured this relative occurrence of bond states in PEG melts using smallangle neutron diffraction and their results agree well with our estimates. There is a simple physical (entropic) reason why one should not find large local helical sequences in nature. As indicated in the introduction, a piece of (EG) 3 has several thousand possible configurations, only two of which are of the type (tgt) 3 .
All results presented in this paper are for PEG in vacuo. There have been speculative arguments in the literature [1] that solvation in water stabilizes local helical structures. We briefly want to address this problem. In a detailed ab initio study [2] of the energetics and geometry of short PEG molecules closely related to the 2/7 helix, it was shown that in the absence of water the helical structure is the most stable one. This result has been confirmed here, but we also show that this conformer is among many of similar energetics and length. In the presence of water it was found [2] that gauche defects (rather than continuous helices) are stabilized greatly by water cross-linking between three to four neighbouring monomers. Thus it appears that solvation in water would destroy rather than enhance local helical structures. Although this argument is sound it needs further elaboration which we hope to do with further ab initio calculations of EG dimers and trimers in the presence of water.
In our last calculation we compute the force-extension curve for a long chain stretched in an AFM, (EG) 200 , of 200 monomers, see figure 6. Again the curve obtained with all conformers included lies above the one where only the (t * t − t * t) conformers are accounted for. The reason is once more that the restriction to a few conformers, i.e. 3 N instead of 3 3N significantly decreases the entropy. In our earlier work we had advanced two scenarios to correct for this flaw, namely to use (EG) 3 as the base unit in a FJC model but modified in such a way that the base unit itself is stretched appropriately with the applied force [2]. The FJC model then obviously supplies the missing entropy. Another scenario was to enlarge the three-state RIS model to a multistate ICM. The motivation was the following: in the RIS model one takes into account only the configuration at the minimum of the energy for a particular conformer. However, small changes in its length will result in a stretched or compressed conformer of the same overall configuration, albeit with somewhat higher energies. The multistate ICM takes a few additional points (five or seven, typically) in this potential trough thus enhancing the total number of contributing configurations and results in a larger entropy.

Conclusion
In this paper we have revisited the statistical mechanics of single PEG with the intention to check ad hoc approximations made in previous works, namely that of keeping all C-O bonds along the chain in the trans configurations. The theory developed here includes all the bonds on an equal footing. Starting with DFT calculations for the (EG) dimer and trimer, one can determine all the energetic and geometric parameters needed for a RIS model with interactions between nearest neighbour monomers included. The results are force-extension curves and end-to-end distributions for short chains. We use the TM method for long chains of 200 EG monomers, i.e. 600 bonds.
A comparison of our theoretical force-extension curve to experiment, figure 6, shows that the agreement of our best theory, namely accounting for all degrees of freedom, is excellent. Recall that this theory is parameter-free because all energetic and geometric parameters were obtained from the DFT calculations. We have also shown that some previously employed approximation schemes are justified. We have also included a discussion of the issue of helical dominance.
The exact analysis approach employed in this work does not, at present, take into account the interaction between PEG and the surrounding solvent as well as intermolecular interactions (such as hydrogen bonds) between different PEG molecules. Based on the results in section 2.2 and section 3 it will be possible to do exact calculations of such a system in the future. An adaptation of an all-conformer sum in the RIS model will facilitate exact calculations on small PEG clusters, say up to three short molecules. Once the interactions between neighbouring molecules is settled, a mean-field theory along the lines developed by Szleifer and co-workers [25] can be used to model the proper thermodynamics at larger scales. Using ab initio calculations such as presented in [15], it might even be possible to account for intermolecular hydrogen bonding in such a theory. This work will have to involve similar methods as employed in the current paper and rigorous checks of all approximations can and must be made along the way.
least-square fitting to all 493 equations of type (A.2), with the results presented graphically in figure A.1(a). The self-energies of all the bonds can be written in vector form as All the values are in meV. Similarly, the nearest neighbour interactions, can be written in matrix form Symmetry requires that the energetics of all g+ and g− states are equal. Moreover, the order of the sequence should not matter. These conditions lead to the following equalities: Note that some combinations were never observed in any of our DFT calculations due to the excluded volume effects discussed previously. Thus, the matrix elements V NN,C−C have been set to infinity.