The four classical, biomolecular force fields designed to study hexopyranose-based carbohydrates (GROMOS 56a6CARBO/56a6CARBO_R, GROMOS 53a6GLYC, CHARMM and GLYCAM06) have been tested in the context of ring-inversion properties. These properties were evaluated for both unfunctionalized monomers of all hexopyranoses of the d series and for residues in a chain composed of uniform units connected by α(1→4) and β(1→4) glycosidic linkages. The results indicate that the tested force fields differ in their predictions of the ring-inversion properties of both monomers and residues in a chain. The comparison with the available experimental data and with the semi-empirical Angyal scheme reveals that, at the level of monomers, GROMOS 56a6CARBO, GROMOS 53a6GLYC and CHARMM correctly reproduce the ring-inversion free energies. However, due to the lack of analogous reference data we cannot state which force field is more or less accurate in the context of ring distortion of residues in a chain. Therefore, the use of ab initio potentials is recommended in the prospective, quantitative studies on the related subject.
The flexibility of pyranose ring is a degree of freedom that often is neglected when studying the three-dimensional structure of carbohydrates. This is because in many cases the ring-inversion free energies are high enough to assume that the ring shape is ‘conformationally locked’ in the energetically favorable state (which, for hexopyranoses of the d series, is the 4C1 conformation, as a rule). However, there is a growing number of reports pointing out that the ring shape and the related conformational rearrangements may play an important role in various, biologically-relevant processes , , , . Furthermore, there exists a significant collection of various hexopyranoses for which a large degree of ring flexibility has been proven in terms of the NMR-based studies, performed for molecules in solution , , , , . In particular, such studies usually concern idose , altrose  and their derivatives ,  as well as some functionalized pyranoses such as mannuronic  or iduronic  acids. Even if the population of non-standard ring geometries is not high, it may significantly increase under influence of external factors. The example are atomic force microscopy (AFM) experiments ,  during which the carbohydrate chain is pulled by the mechanical force. Contrary to the case of monomers, there are no experiental studies that would quantitatively describe the ring flexibility for residues in carbohydrate chain and the available data are either of qualitative nature , ,  or concern the static crystal structures (see discussion in Ref. ). Currently, there is no quantitative interpretation of the AFM ‘stretching curves’ in terms of population of particular ring conformers in the pulled chain. Thus, the results of such experiments (although giving an important assessment of the favorable ring conformation), provide only indirect evidences which can hardly be translated into the quantitative ring-inversion properties , .
Because of scarce number of experimental reports that provide the quantitative description of hexopyranose ring flexibility, the theoretical approaches based on the molecular dynamics (MD) technique seem to be an attractive alternative to investigate this subject. The results of recent theoretical studies suggest that the ring flexibility may be enhanced when inserting a given residue in a chain of specific linkage type or upon functionalization (β-O1-methylation; see Fig. 1a for atom numbering) . Furthermore, the ring shape may influence the favorable orientation of glycosidic linkages while the significance of this effect depends not only on the ring-inversion properties but also on the linkage type , .
The accuracy of the MD-based methods relies primarily on the quality of potentials that approximate various interactions present in the molecular system of interest. In spite of several cases of successful application of the ab initio potentials into MD calculations aimed at studying the carbohydrate conformation , , the high computational cost makes such attempts limited to molecular systems of small size. Larger number of studies focused on puckering of the hexopyranose ring were conducted by using classical, biomolecular force fields dedicated to study the conformational properties of carbohydrates , , , , , , , , . So far, the accuracy of such force fields has been comparatively studied only in the context of small sets of molecular systems (e.g. β-d-glucopyranose)  or only for unfunctionalized hexopyranose monomers . The collection of all hexopyranoses of the d series either in monomeric forms or connected by α(1→4) and β(1→4) linkages has never been considered.
The main aims of our investigation are:
Describing and comparing the ring inversion properties of unfunctionalized monomers of all hexopyranoses of the d series calculated within four different force fields (GROMOS 56a6CARBO/56a6CARBO_R , , GROMOS 53a6GLYC , CHARMM ,  and GLYCAM06 ). Contrary to already existing studies (see e.g. table 4 in Ref. ) these properties will be expressed here as the two separate quantities: (i) the free energy change associated with the full chair-chair ring inversion from 4C1 to 1C4 (ΔFi); (ii) the free energy change associated with the ‘half’ ring inversion from 4C1 to boat/skew-boat conformers (B/S) (ΔFh).
Describing and comparing the ring inversion properties of residues in a chain composed of uniform monomers connected by glycosidic linkages of the α(1→4) and β(1→4) types calculated within three different force fields (GROMOS 56a6CARBO/56a6CARBO_R, GROMOS 53a6GLYC and CHARMM). In this stage of the study the ring inversion properties are also investigated in the context of both the ΔFi and ΔFh values.
All MD simulations were carried out with the GROMACS 4.5.5 program , using the following four force fields: GROMOS 56a6CARBO/56a6CARBO_R , , GROMOS 53a6GLYC , CHARMM ,  and GLYCAM06 . Note that one of these force fields exists in two versions, namely: GROMOS 56a6CARBO and GROMOS 56a6CARBO_R; the latter version contains parameters applicable for simulating the alkylated and glycosylated carbohydrates whereas parameters for unfunctionalized monomers are identical in both variants. For the sake of simplicity, in the present article we will be using only the 56a6CARBO term. The GROMACS-compatible CHARMM parameters were downloaded from the official CHARMM website (http://mackerell.umaryland.edu/charmm_ff.shtml); the GROMOS 56a6CARBO_R parameters were prepared during the procedure of the force field parameterization (see Ref. ); the GROMOS 53a6GLYC parameters were downloaded from the website of Biomaterial Modeling Group of Federal University of Pernambuco (http://dqfnet.ufpe.br/biomat); the GLYCAM parameters for monomers were provided by Dr V.H. Rusu. In the case of trisaccharides, the pdb2gmx program was used to produce the molecular topologies from respective monomeric building blocks.
For each of these force field an appropriate water model was applied to describe the explicit solvent present in the system, i.e. SPC  (GROMOS force fields) or TIP3P  (GLYCAM06 and CHARMM). The systems considered consisted of one solute molecule solvated by water molecules within a cubic computational box simulated under periodic boundary conditions at constant volume and at a temperature of 298 K. Given the selected numbers of water molecules (~900 to ~1500 for the different systems), the box edges (3.1–3.6 nm for the different systems) were all preoptimized by a 1 ns constant-pressure MD equilibration at 1 bar and 298 K, ensuring an effective solvent density appropriate for these conditions in the subsequent isochoric production simulations.
The solutes considered were mono- or trisaccharides. The monosaccharides included the eight hexopyranoses of the d-series, namely allose (All), altrose (Alt), glucose (Glc), mannose (Man), gulose (Gul), idose (Ido), galactose (Gal) and talose (Tal), considering both the α- and β-anomers. The trisaccharides considered systematically were α(1→4)- and β(1→4)-linked homooligomers of the eight monosaccharides without any functionalization, the anomeric configuration at the reducing end being the same as that of the linkage.
The calculations involved the determination of one-dimensional free-energy profiles (FEPs) along the Cremer-Pople θ parameter  describing the conformation of specific rings. For the FEP analyses along θ, the intervals [0–60°], [60–120°] and [120–180°] were ascribed to regular chair (4C1), boat/skew-boat (B/S) and inverted chair (1C4) conformations, respectively, as illustrated in Fig. 1. The free-energy changes (ΔFi and ΔFh) associated with the ring inversions (4C1→1C4 and 4C1→B/S, respectively) were defined as the difference between the values at the FEP minima in the two corresponding regions. These analyses were carried out systematically for the monomers and for the central residue in the trimers.
The calculation of FEPs relied on an enhanced-sampling scheme combining  parallel tempering  and well-tempered metadynamics  as implemented in the PLUMED 1.3 plug-in  of the GROMACS program . The well-tempered metadynamics  relied on Gaussian local functions of widths 2.86°, an initial deposition rate of 0.01 kJ mol−1 ps−1 and a temperature parameter ΔT (see eq. (2) in Ref. ) of 7152 K (simulations within the GROMOS 53a6GLYC force field) or 1788 K (remaining force fields). The parallel-tempering  relied on 16 metadynamics simulations carried out in parallel at different temperatures ranging from 298.0 to 363.2 K in steps of about 4.3 K, along with replica-exchange attempts performed at 2 ps intervals.
After equilibration, all simulations were carried out until the convergence of ΔFi and ΔFh values were detected which resulted in total simulation time varying from 10 ns (GROMOS 56a6CARBO, GLYCAM06) to 110 ns (GROMOS 53a6GLYC). The convergence of the ΔFi and ΔFh values was checked during the course of simulations by using the PLUMED built-in sum_hills tool and hand-written scripts.
In all simulations, the equations of motion were integrated using the leap-frog scheme  with a timestep of 2 fs. The full rigidity of the water molecules was enforced by application of the SETTLE procedure . The temperature was maintained close to its reference value by application of the V-rescale thermostat  with a relaxation time of 0.1 ps. The translational center-of-mass motion was removed every timestep. The solute bond lengths were constrained by application of the LINCS procedure ,  with a relative geometric tolerance of 10−4 to all covalent bonds (the GROMOS force fields) or to hydrogen-containing bonds (GLYCAM06 and CHARMM). In the case of GROMOS 56a6CARBO, the non-bonded interactions were calculated using a twin-range scheme , with short- and long-range cutoff distances set to 0.8 and 1.4 nm, respectively, and an update frequency of five timesteps for the short-range pairlist and intermediate-range interactions. A reaction-field correction ,  was applied to account for the mean effect of the electrostatic interactions beyond the long-range cutoff distance, using a relative dielectric permittivity of 61 . In the case of GROMOS 53a6GLYC, long-range electrostatic interactions were treated by the reaction-field method with a relative dielectric permittivity of 66. A 1.2 nm cutoff was used for the short range electrostatics and Lennard-Jones interactions. Within GLYCAM06-based simulations long-range electrostatic interactions were treated by using the particle mesh Ewald summation  (a 0.9 nm distance cutoff was used for direct space nonbonded calculation). In the case of the CHARMM force field, a force-switched smoothing function  was applied to Lennard-Jones interactions in the range of 1–1.2 nm. Long-range Coulomb interactions were handled using particle mesh Ewald with a real-space cutoff of 1.2 nm.
The results of simulations are gathered in Tables 1 and 2 as the free energy values ΔFi and ΔFh. The correlation between selected sets of data is presented graphically in Figs. 2–5. All possible monomers of d-hexopyranoses as well as homotrimers composed of corresponding residues were considered, except of the case of GLYCAM06 where, due to the lack of available parameters, investigations were limited to non-Ido monomers.
The calculated values rely on Cremer-Pople coordinates  for the state assignment. The correlation between the selected values is illustrated graphically in Figs. 2–5.
Let us mention that both force fields belonging to the GROMOS family (56a6CARBO and 53a6GLYC) have been parametrized with attention paid to possible deformations of pyranose rings. The respective papers (Refs. , , ) contain a detailed description of the hexopyranose ring-inversion properties which were studied in terms of all possible ring conformers  or just the chair-chair inversion , . Initially, these analyses have been performed only at the level at which the ring flexibility-related features were adjusted, i.e. for unfunctionalized monomers; therefore, in the case of GROMOS 56a6CARBO, the additional patch was needed to account for the more accurate description of the ring distortion in O1-functionalized pyranoses . Conversely, GROMOS 53a6GLYC, GLYCAM06 and CHARMM have never been systematically tested with respect to the ring-puckering properties in the context of varying local environment of the distorted ring (i.e. unfunctionalized monomer vs. functionalized monomer vs. residue in a chain).
As the experimental data providing the desired quantities (i.e. the relative populations of at least the two most abundant ring conformers) are limited to idose  and altrose , for the rest of hexopyranoses the semi-empirical Angyal scheme  was used for comparative purposes (Table 1). Note that this scheme concerns only unfunctionalized hexopyranose monomers and the chair-chair inversion energy. Thus, for the second type of considered data, i.e. free energies of chair-boat/skew-boat transition, there is neither experimental data nor semi-empirical schemes that could be used to verify the results obtained by MD simulations within different force fields. The same is true for any type of ring inversions occurring for residues in a chain.
Both the GROMOS 56a6CARBO and CHARMM force fields show an excellent agreement with both the Angyal scheme prediction and (where available) experimental data. The largest divergences are observed for compounds that are assumed to exhibit very rigid rings, independently of the method of calculations, e.g. β-Gal, β-Tal or α-Man. In the case of experimental data the largest divergence corresponds to the 56a6CARBO force field and β-anomers of Ido and Alt and equals to ~6 kJ/mol. Note that for each compound the sign of the ring-inversion free energy remains uniform with that of the Angyal scheme and/or experimental data.
A similar level of qualitative agreement is offered by the GROMOS 53a6GLYC force field. Also in this case the sign of ΔFi values remains in agreement with the Angyal scheme prediction and the NMR data. However, now a much broader dispersion of the ring-inversion free energies is suggested. This results in both more extreme values of ΔFi observed for most of compounds (|ΔFi|>20 kJ/mol for 10 out of 16 compounds) and larger divergences between 53a6GLYC and Angyal scheme which, in the extreme case of β-Glc, can reach ~50 kJ/mol. Moreover, in spite of the correct sign of ΔFi for α-Ido, the absolute difference between the calculated and NMR-inferred values equals to ~25 kJ/mol. The related difference for β-Alt is also large (~12 kJ/mol).
The least degree of agreement between MD simulations and the available experimental data and/or the Angyal scheme was exhibited by the GLYCAM06 force field. GLYCAM06 is clearly at odd with the experimental data concerning α-Alt and suggests strongly dominating population of 1C4 conformers (ΔFi<−25 kJ/mol). Furthermore, GLYCAM06 predicts the negative values of ΔFi for α-All and α-Gul which is rather surprising in the context of both the semi-empirical scheme-inferred data  and the ab initio calculations . The lack of parameters for Ido additionally hinders the comparison with the second set of existing experimental data.
The reasons for relatively worse GLYCAM06 performance in comparison with other force fields are most likely associated with the method of parametrization accepted during its development. Contrary to both GROMOS force fields and CHARMM, each molecule of unfunctionalized hexopyranose modeled by GLYCAM is associated with its own unique set of parameters (e.g. partial atomic charges) developed on the basis of ensemble-averaged conformation. As the dominating conformer is usually 4C1, the derived charges are suited very well to describe the ensemble of states associated with the regular chair conformation but may fail to represent alternative ring shapes as well as their relative populations. GROMOS and CHARMM are designed to keep the set of parameters identical or nearly identical with respect to all hexopyranoses, thus, they suffer less from this effect.
Let us emphasize that the divergence between predictions of the given force field and the other types of data may influence the force field applicability only in some very special cases. This is due to the following reasons: (i) the behavior of the most frequently studied hexopyranose monomers (Glc, Man, Gal) is roughly independent of the applied force field and corresponds to strongly dominating population of the 4C1 conformers; (ii) the free energy barriers separating particular ring conformers are high enough to prevent spontaneous ring ‘flips’ at the timescale of tens or even hundreds of ns of unbiased MD simulation.
Finally, let us refer to the ring-inversion free energies reported elsewhere . The differences between our calculations and the data shown in table 4 of Ref.  are of quantitative nature and can be ascribed to different schemes of metadynamics-based calculations. Namely, calculations reported in Ref.  employed ‘plain’ metadynamics, without bias scaling, whereas in the present case, the well-tempered metadynamics scheme was used in combination with the replica-exchange method. Appreciating the existence of these minor differences, let us note that also in this case, the signs of ΔFi values remain uniform.
The tested force fields, compared in the context of correlation between ring-inversion properties of monomers and those of residues in a chain, reveal a series of differences in their predictions. As discussed in details in Ref. , in the case of GROMOS 56A6CARBO, the values of ΔFi are tendentially reduced in magnitude upon insertion of the free residue into a chain (for 12 out of 18 trisachcarides). Furthermore, the correlation between the ΔFi values for the monomers and for the corresponding residues in a chain is very poor along the d-hexopyranose series. In contrast, there exists a much higher correlation between the same quantities calculated within both CHARMM and GROMOS 53a6GLYC force fields. While the monomer- and trimer-related ΔFi values for β-anomers are very close to each other (except of decrease of ΔFi by ~20 kJ/mol for β-Alt and β-Tal simulated within GROMOS 53a6GLYC), the tendential decrease of the trimer-related ΔFi values is observed for 4–6 out of 8 α-anomers (i.e. for Glc, Gal, Man and Tal in the case of CHARMM and Glc, Gal, Man, Tal, Gul and Alt in the case GROMOS 53a6GLYC). Interestingly, Glc, Gal, Man and Tal are the same compounds for which also the GROMOS 56A6CARBO force field predicts the largest decrease of ΔFi upon inserting residue into a α(1→4)-linked chain. All of these compounds contain equatorial hydroxyl group at C3 in the 4C1 chair (see Fig. 1 for notation). Note that both the CHARMM and GROMOS 56A6CARBO force fields agree in their predictions that the dominating ring conformation at the level of (non-Ido) monomer is the same as that of residue in a chain. On the contrary, in the case of GROMOS 53a6GLYC, there are three compounds for which the monomer- and trimer-related ΔFi values differ in their sings, i.e. α-Alt, α-Gul and α-Tal. The reason for that is a tendency to systematic decrease of ΔFi exhibited by GROMOS 53a6GLYC in the case of α(1→4)-linked residues. Such decrease in combination with the residues of relatively small, monomer-related, ΔFi values results in the negative values of ring-inversion free energies for residues in a chain.
Analogous investigations focused on the possible correlation between the ΔFh values calculated for the monomer and for the central residue in trimer give similar results, i.e. a poor correlation for GROMOS 56A6CARBO and much better ones for CHARMM and GROMOS 53a6GLYC. Interestingly, ΔFh values characteristic of residues in a chain are again decreased in comparison to those of monomers for the same group of four α anomers containing equatorial hydroxyl group at C3 in the 4C1 chair conformer (the related decrease is additionally exhibited by α-Gul in the case of GROMOS 53a6GLYC).
There is no quantitative experimental data that would allow to assign the ring-inversion free energy value to the residue being a part of a larger chain. The existing qualitative data determining the preferred conformation of such residues concern only selected systems and agree with the predictions of all tested force fields (e.g. Refs. , ,  suggesting the 4C1 chair conformation as dominating in residues of the amylose chain). Similarly, there is no semi-empirical schemes generalized for the case of glycosylated hexopyranoses and no ab initio methods have been used so far to systematically study this issue. For these reasons, currently there is no possibility to compare the quality of the force fields on the basis of the data collected in Table 2. However, the most striking difference between force fields is the presence or the absence of the correlation between ring-inversion properties of monomers and those of residues in a chain. This is especially important in the context of the lack of quantitative data concerning the ring conformation of residues in a chain and the possibility of predicting of such properties on the basis of relatively better recognized properties of respective monomers.
The mentioned difference can be explained when considering the monomer-inherent flexibility of the rings predicted by each of the tested force fields. The chair-chair ring-inversion free energies suggested by GROMOS 56A6CARBO for monomers vary from ~−2 to ~28 kJ/mol (i.e. in the range of ~30 kJ/mol). The corresponding ranges for the CHARMM and GROMOS 53a6GLYC force fields are larger and equal to ~40 kJ/mol and ~100 kJ/mol, respectively. Let us assume that the non-systematic alterations of the ring flexibility upon inserting residue into a chain, as suggested by GROMOS 56A6CARBO, originate from the residue-residue interactions and that the strength of such interactions is of comparable magnitude for each of the force fields. Then, the ‘decorrelation effect’ will be more significant in the case of force field that predicts larger flexibility of rings, i.e. GROMOS 56A6CARBO. The residue-residue interactions of the same strength are less important for rings of greater degree of rigidity, i.e. those described in terms of CHARMM and, especially, GROMOS 53a6GLYC. Such hypothesis would link the possible decorrelation of ring-inversion free energies with the capability of the tested force fields to reflect very subtle differences between ring puckering at the level of monomers. As most of the currently available data allow only to assign a particular ring conformation as the dominating one and are not capable to characterize this domination by providing the ring-inversion free energy (or equivalent quantity), the possible alternative methods of choice would include MD-related techniques employing ab initio potentials of increased accuracy (see e.g. Ref. , ). However, due to the high significance of intramolecular hydrogen bonding  that strongly influences the conformation of hexopyranoses in the absence of solvent, such simulations would require the presence of explicit solvent, possibly modeled at the molecular-mechanics level of theory.
Interestingly, there exists a relatively good correlation of the chair-chair inversion free energies (ΔFi) values with those of chair-boat/skew-boat ones (ΔFh). The correlation is present for all tested force fields, however, it becomes very poor in the case of residue in a chain simulated within GROMOS 56A6CARBO and is of (better) comparable quality for the remaining cases (see Figs. 4 and 5). In the case of monomers, with the exception of very rigid rings (ΔFi>20 kJ/mol), the relation ΔFi≤ΔFh is always satisfied. For residues in a chain the same inequality holds except of α- and β-Ido simulated within the CHARMM force field. This means that the moderately flexible rings should be described in terms of the 4C1 → 1C4 transition rather than the 4C1 → B/S one. For more rigid rings the population of B/S conformers may be larger in comparison to 1C4, however both of them will be very scarce in relation to 4C1. Finally, the most flexible rings (All and Ido, according to predictions of all force fields) may require a more complex description involving both types of transitions or, possibly, introducing additional coordinate(s) accounting for the ring shape within the B/S set of conformers (e.g. Cremer-Pople φ).
Four different classical force fields have been tested in the context of their capabilities to simulate the ring distortion in hexopyranose molecules (in monomers and in residues in a chain exploiting either α(1→4) or β(1→4) glycosidic linkages). Both the GROMOS 56a6CARBO and CHARMM force fields offer very similar predictions for unfunctionalized monomers which, at the same time, remain in agreement with the available experimental data and with the semi-empirical Angyal scheme. GROMOS 56a3GLYC also reflects this agreement at the qualitative level but differs quantitatively from GROMOS 56a6CARBO and CHARMM by predicting a much broader dispersion of ring-inversion free energies varying significantly from one compound to another. GLYCAM06 exhibits the worst accuracy at the level of unfunctionalized monomers, being at odd with the NMR data for α-Alt and predicting the dominating population of inverted chairs for α-Gul and α-All. The most frequently studied hexopyranoses (Glc, Man, Gal) exhibit similar ring-inversion properties, independently of the force field used. All force fields predict alterations of the chair-chair inversion free energies upon inserting residue into a chain. However, the scale and direction of such alterations are dependent on the force field: GROMOS 56a3GLYC and CHARMM suggest a strong correlation between monomer- and trimer-related properties and increased flexibility of some of the α(1→4)-linked residues, relative to the properties of corresponding monomers. According to GROMOS 56a6CARBO, the ring-inversion free energies for residues in a chain are typically smaller in magnitude compared to those of the monomers (for both the α(1→4)- and β(1→4)- linkages) and correlate poorly with the latter. Except of very rigid rings, chair-boat/skew-boat inversion free energies are higher than those corresponding to the chair-chair inversion. For all tested force fields, there exists a correlation between chair-chair and chair-boat/skew-boat inversion free energy for both monomers and residues in a chain. However, it becomes rather poor for trimers simulated within GROMOS 56a6CARBO.
Finally, as the availability of the reliable quantitative data related to the ring distortion of residues in a chain is very limited, the use of the force field-independent, ab initio potentials in combination with the presence of explicit solvent may be suggested as an alternative method to study this issue.
A collection of invited papers based on presentations at the XXVIII International Carbohydrate Symposium (ICS-28), New Orleans, July 17–21 2016
The authors acknowledge the financial support of the Polish National Science Centre (contract financed in 2012–2015 under Project No. 2011/03/D/ST4/01230 SONATA). We thank Dr V.H. Rusu for providing the topology files of some of the molecular structures used in this work.
 O. Guvench, S. N. Greene, G. Kamath, J. W. Brady, R. M. Venable, R. W. Pastor, A. D. MacKerell, Jr. J. Comput. Chem.29, 2543 (2008).10.1002/jcc.21004Search in Google Scholar PubMed PubMed Central
 K. N. Kirschner, A. B. Yongye, S. M. Tschampel, J. Gonzalez-Outeirino, C. R. Daniels, B. L. Foley, R. J. Woods. J. Comput. Chem.29, 622 (2008).10.1002/jcc.20820Search in Google Scholar PubMed PubMed Central
 S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. Van Der Spoel, B. Hess, E. Lindahl. Bioinformatics29, 845 (2013).10.1093/bioinformatics/btt055Search in Google Scholar PubMed PubMed Central
 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, J. Hermans. “Interaction models for water in relation to protein hydration”, in Handbook of Intermolecular Forces, B. Pullman (Ed.), Springer-Science+Business Media, B. V. Israel (1981).10.1007/978-94-015-7658-1_21Search in Google Scholar
 M. Bonomi, D. Branduardi, G. Bussi, C. Camilloni, D. Provasi, P. Raiteri, D. Donadio, F. Marinelli, F. Pietrucci, R. Broglia, M. Parrinello. Comput. Phys. Comm. 180, 1961 (2009).10.1016/j.cpc.2009.05.011Search in Google Scholar
 R. W. Hockney. Methods Comput. Phys.9, 136 (1970).Search in Google Scholar
 B. Hess, H. Bekker, H. J. C. Berendsen, J. G. E. M. Fraaije. J. Comput. Chem.18, 1463 (1997).10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0.CO;2-HSearch in Google Scholar
©2017 IUPAC & De Gruyter. This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License. For more information, please visit: http://creativecommons.org/licenses/by-nc-nd/4.0/