J Mol Model (2014) 20:2226 DOI 10.1007/s00894-014-2226-4

ORIGINAL PAPER

An innovative synergistic grid approach to the computational study of protein aggregation mechanisms Noelia Faginas-Lago & Margarita Albertí & Alessandro Costantini & Antonio Laganà & Andrea Lombardi & Leonardo Pacifici

Received: 25 November 2013 / Accepted: 26 March 2014 # Springer-Verlag Berlin Heidelberg 2014

Abstract Thanks to the advances in grid technologies, we are able to propose here an evolution of our molecular simulator that, when moving to larger systems, instead of reducing the granularity of the dynamical treatment (as is often done in molecular dynamics studies of such systems) exploits the extra power of the grid approach to the end of preserving the detailed nature of theatomistic formulation of the interaction. Key steps of such evolution are: (1) the assemblage of the interaction based on a composition of the ab initio intramolecular data and a portable parameterization of the intermolecular potential linking ab initio evaluation of intramolecular potentials and the partitioning of molecular polarizability; (2) the exploitation of an efficient coordinated porting and running of molecular dynamics codes on the European grid distributed computing infrastructure. As a prototype case study, the N-methylacetamide dimer in vacuo has been considered and the formation of possible conformers is analyzed.

Keywords Intermolecular interaction . Grid infrastructure . DL_POLY . Molecular dynamics

This paper belongs to Topical Collection 9th European Conference on Computational Chemistry (EuCo-CC9)

Introduction The progress of computer technologies is depriving computational sciences (among which are the molecular ones) of the role of being the driving force for computer technology evolution and converting it into that of being an advanced test field. This is indeed the case of computational investigations of systems of biological interest whose feasibility is increasingly dependent on the development of numerical procedures for accurately computing short- and long-range interactions, fast and slow molecular motions, and detailed and averaged observable properties specifically tailored for modern concurrent architectures. In order to deal with such molecular dynamics (MD) investigations and to have access to the necessary computing time, we have decided to further extend the so-called Grid-Empowered Molecular Simulator (GEMS). GEMS was born out of the original idea discussed in [1, 2] and was later implemented as a collaborative workflow of a set of codes based on different types of expertise [1–3]. It was designed with the purpose of offering users the possibility of carrying out ab initio simulations of the measured properties of realistic molecular systems (virtual experiment [4]). For this reason, GEMS is structured to offer as a service, in an a la carte approach programs, packages, and data of the following computational modules: &

N. Faginas-Lago (*) : A. Laganà : A. Lombardi : L. Pacifici Dipartimento di Chimica, Biologia e Biotecnologie, University of Perugia, Via Elce di Sotto, 8, 06123 Perugia, Italy e-mail: [email protected] M. Albertí Departament de Química Física, Universitat de Barcelona, Barcelona, Spain A. Costantini INFN-Perugia and Italian Grid Infrastructure, Perugia, Italy

&

Interaction to fully or partially evaluate, if not otherwise available, the potential energy values of the investigated systems. When needed, also a suitable functional form is fitted to the ab initio potential energy values. In order to facilitate interoperability among different quantum chemistry codes, the Q5cost [5] de facto standard format is adopted to store data-related files. Dynamics to integrate the equations determining the dynamics of the nuclei of the system. In a rigorous approach, the problem is dealt using full-dimensional quantum mechanics techniques that are the method of election for an

2226, Page 2 of 9

&

exact calculation of chemical reactivity. As an alternative, classical mechanics equations can be integrated for more approximate treatments. In this module, the D5cost (dynamics) [6] counterpart of Q5cost is adopted. Observables to perform the statistical treatments of the outcomes of the computed dynamical calculations needed to provide a priori estimates of the measured properties of the system. This implies a proper sampling of both initial and final conditions, a continuous evaluation of some intermediate quantities, and a final integration over some unobserved variables.

In particular, the synergism is obtained by combining expertise and software useful for: a) exploiting the advantage of using concurrently, on the distributed network of the grid, a large number of CPUs and cores; b) utilizing available ab initio information on the electronic structure of the considered molecular system (or generate them afresh using one of the packages available on the grid) originating from different sources; c) formulating the electronic energy values using suitable functional forms and producing a subroutine embodying analytical potential energy surfaces (PESs) d) incorporating the PES routine inside the dynamics code and performing extended runs of such code for evaluating related outcomes for the necessary range of initial conditions; e) statistically monitoring and averaging dynamical outcomes to evaluate experimental observables. To this end, the GEMS synergistic approach was adopted not only for exploiting the distributed features of the grid platform, but also for properly sampling the experimental conditions including, when needed, the geometrical characteristic of the experimental apparatus as is the case the cross beam virtual experiment of [4]. For more complex molecular systems relevant to biological studies, as is the case of the system discussed here, additional problems related to the above-mentioned b, c, and d items of the synergistic model had to be addressed. In particular, a more adequate representation of ab initio data and a more flexible formulation of the PES suited for properly dealing with the long-range forces acting on the system at low temperature had to be developed. The case study we consider here is the N-methylacetamide (NMA) dimer system (NMA–NMA) in vacuo. The long-term scientific target of such a choice is the use of GEMS for investigating the rearrangement of molecular systems made of a large number of proteins (in particular of their ordering in structures having specific functions). These effects involve typically electrostatic, hydrogen-bond and van der Waals (or dispersion) interactions which are, arguably, weaker than

J Mol Model (2014) 20:2226

conventional chemical bonds and more difficult to quantify both experimentally and theoretically. They also involve dynamical effects associated with mechanisms that go well beyond purely statistical and geometrical factors and are the manifestation of specific features of the long-range tails of the interaction. As a matter of fact, their study is important for the understanding of several fundamental and applicative issues like the formation of nanostructures in solution [7, 8], the alignment and orientation in molecular collisions, the stereodynamics of chiral molecules (see e.g., [9, 10]), the efficiency of reactive systems (see e.g., [8, 11, 12]), the sizeselectivity of transport properties of ion channels (see e.g., [13]). In this paper, we give a detailed account of the work carried out by extending the synergistic ab initio approach of GEMS to the understanding of the formation of NMA–NMA protein dimer building blocks. Accordingly, the paper is organized as follows: in section 2 the assemblage of the NMA–NMA PES is illustrated, in section 3 the outcomes of the MD study of the forming of NMA–NMA conformers are discussed. Some conclusions are drawn in section 4.

The synergistic approach to NMA–NMA electronic structure, data formats, and representations As already mentioned, the first step of our approach to highercomplexity molecular simulations is the Interaction module of GEMS that has been adapted to deal synergistically with an appropriate combination of ab initio energies at short range and empirical potentials at long range. For this purpose, the global interaction V is decomposed, as usual, into Vintra (intramolecular) and Vinter (intermolecular) components, which depend on the decomposition of the molecular polarizability [14] in effective potentials associated with atoms and bonds (or groups of atoms) of the molecular system [8, 11]. Such a model is grounded on past experiences concerning different neutral [15] and ionic [16, 17] systems, often involving weak interactions [15, 18], whose reliability has been assessed by comparison with different ab initio calculations [19]. Such adaptation of the Interaction module of GEMS, in addition to the already existing option of formulating the electronic structure of the stable fragments, Vintra, using Hartree–Fock (HF), post-Hartree–Fock, density functional theory (DFT) [20] packages offers the possibility of accessing databases hosting already computed ab initio values. The combining of ab initio data originating from different sources faces serious interoperability problems because of the different formats used for the electronic structure calculations and causes a disruption of the GEMS synergism. Compared with other molecular science applications, in fact, the quantum chemistry codes lag quite behind in terms of data format standardization. This is indeed the case, for example, of the

J Mol Model (2014) 20:2226

Page 3 of 9, 2226

widespread protein [21] and crystallographic structures [22] data banks. In the work reported here, this has been obviated by resorting into the adoption of the already mentioned de facto Q5cost standard and its development into the dynamicsoriented D5cost one with related libraries [6] when referring to the framework of the Born–Oppenheimer (BO) scheme [23] and to wave function-based electronic structure computing approaches. In GEMS, the usual utilization of ab initio data in dynamics calculations, when they are not used “on the fly”, is by best-fitting the parameters of empirical functional forms of the potential in its force field-like formulations. Among the popular force-field parameterizations usually considered in MD calculations, [AMBER [24] (Assisted Model Building and Energy Refinement), CHARMM [25] (Chemistry at HARvard Molecular Mechanics) and OPLS [26] (Optimized Potential for Liquid Simulations)], we adopted AMBER. For V inter, the GEMS synergism is further strengthened by the portability of the adopted force-field formulation (see the following discussion specifically tailored for the NMA dimer). This is vital when carrying out computationally intensive investigations of processes in which various proteins act as structural building blocks of biological systems and as catalyzers of large families of chemical reactions. In the particular case of the NMA–NMA dimer considered here (that is a minimal model of the peptide linkage forming the backbone of the proteins), the attention has been focused on the fact that the behavior of the NMA–NMA interaction in aqueous solution regulates the competition between peptide–peptide and peptide–water hydrogen bonds arising in protein hydration. The NMA hydrogens bind through the amide and carbonyl groups to form linear and branched chains in both pure liquid and aqueous solutions. The portability of such molecular-level description of protein– protein interaction is crucial to the understanding of aggregated protein functions and is central to our paper. More specifically, for Vintra we adopted, as already mentioned, the AMBER parameterization [27], while for Vinter the interaction was partitioned into a VvdW (van der Waals) and a Vel (pure electrostatic) component V inter ¼ V el þ V vdW

ð1Þ

In our case, reference was made to the trans-NMA conformer because of its higher abundance and stability [28, 29]. For this conformer, a specific FORTRAN routine (named coNMA.f) based on Eq. (1) was designed and implemented. In coNMA.f, VvdW explicitly considers a certain number of interaction centers (depending on the complexity of the considered molecules) placed on some atoms (or groups of atoms) and incorporates, as well, exchange (size) and dispersion contributions. Five interaction centers were used to describe the NMA interaction. In particular, two of such centers are associated with the two CH 3 groups and are placed on the related C (C N and C C in Fig. 1) atoms. Two other interaction

centers are associated with (and placed on) the C and O atoms of the carbonyl group (unlabeled C and O of Fig. 1). The last interaction center (associated with the NH group and the related lone pair) was placed on the N atom (N of Fig. 1). A polarizability contribution, whose sum reproduces the value of the molecular polarizability, was assigned to each center (see Table 1). For the NMA dimer, the portability of VvdW, which is expressed as a sum of 15 effective pairs of interactions formulated as an Improved Lennard Jones (VILJ) functional form, [31, 32] is the following:  V ILJ ðrÞ ¼ ε

m  r0 nðrÞ nðrÞ  r0 m − nðrÞ−m r nðrÞ−m r

 ð2Þ

In Eq. (2), ε, r0, and m are pair-specific parameters, while r is the distance between the two centers. The first term of the bracketed sum in Eq. (2) (the positive one) represents the sizerepulsion contribution arising from each pair, while the second term (the negative one) represents the effective dispersion attraction of the same pair. The n(r) exponent of the first term shapes its falloff as a function of r and is formulated as follows: nðrÞ ¼ β þ 4:0

 2 r r0

ð3Þ

with β being an adjustable parameter that introduces more ambient-like metrics (often named as “hardness” of the interacting partners [31, 32]) by modulating the repulsion and controlling the strength of the attraction. The introduction of this modulation (absent in the usual Lennard Jones (VLJ) potential) provides VILJ with the possibility of indirectly taking into account of induction, charge transfer, and atomclustering effects. This additional (with respect to VLJ) β parameter of VILJ corrects the dependence of the interaction on the internuclear distance and removes most of the VLJ inadequacies in the asymptotic region [33]. Such an approach has already shown in the recent past to be portable to other types of clusters [8, 19, 33–37]. The extra adaptability of VILJ is exploited to the end of increasing its portability by allowing the use of the same values of ε (the potential well depth) and r0 (the equilibrium distance) for different diatoms, if the same interaction centers are considered [38]. This assigns to ε and r0 a structural pair interaction meaning, while leaving to β the role of embodying the contributions of more collective and ambient effects. The values of ε and r0, adopted for NMA– NMA, are given in Table 2, using the same labeling of Fig. 1. The values were derived from the polarizabilities of the isolated atoms and from the effective group polarizabilities of CH 3 and NH (see Table 1) by applying the recipes of [14]. Such a choice guarantees the correct reproduction of molecular polarizabilities with the internal consistency of the parameters being enforced by the adopted procedure.

2226, Page 4 of 9

J Mol Model (2014) 20:2226

Fig. 1 (Left) the structure of NMA sketched as a peptide bond capped on both ends by a methyl group. Atoms belonging to the methyl group bonded to the carbonyl carbon are labeled by the subscript C, while those bonded to the nitrogen atom are labeled by the subscript N. The three

potential hydrogen-bonding sites are marked with dashed lines. (Right) a visualization of an NMA molecule taken from an actual simulation run. The color scheme is as follows: hydrogen (white); carbon (green); oxygen (red); nitrogen (blue)

Vel is formulated, as usual, as a sum of Coulombic potentials associated with a set of point-wise charges localized on the atoms of the corresponding molecule having a defined spatial distribution. For non-elementary systems, like the NMA aggregates, the choice of the spatial distribution of the charge bears a certain extent of arbitrariness [39] when different isomers have the same energy. The adopted point-wise charge values used for the calculation and quoted in Table 3 are those of AMBER because they lead to the most stable dimer (see [19]). The resulting value of dipole moment is 5.567 D. As pointed out in [40], it is worth stressing out here that the value of β is obtained starting from an initial guess set equal to the cubic root of the polarizability of the involved partners. Then, in order to account for the coming into play of charge transfer effects (as is the case of H bonding), the value of β is varied against the binding energy when different sets of pointwise charge values are adopted using the MD calculations iteratively.

equations of the DYNAMICS module. This module too has been structured for a synergistic use of some general purpose packages like GROMACS [41], NAMD [42], and DL_POLY [43]. The one used in this work is DL_POLY that can be used to simulate a wide variety of molecular systems. In order to better manage the complexity of its implementation into GEMS, a set of tools facilitating the use of different ensembles of data, models, and programs for an increasing number of applications have been developed. Of particular importance among them has been the use of the P-GRADE Grid Portal [44–47] and GriF [48, 49] framework. The P-GRADE Grid Portal is an open-source tool that provides intuitive graphical interfaces for the porting of computer codes without necessarily requiring the modification of their original structure for a distributed execution on grid platforms. GriF, instead, is a Service-Oriented Architecture Collaborative Framework designed to facilitate the use of

The synergistic NMA–NMA molecular dynamics calculations: porting, performances, and results Once completed, the formulation of the NMA–NMA potential GEMS undertakes the integration of the classical mechanics Table 1 Polarizability values assigned to the interaction centers [30]

Table 2 Values of the well depth (ε), equilibrium distance (r0), and m parameters for the NMA–NMA interactions NMA–NMA Pairs

ε/meV

r0/Å

m

C C -C C C C -C N C N -C N C-C O-O N-N C C -C C N -C C C -O C N -O C C -N C N -N C-O C-N O-N

12.64 12.64 12.64 6.52 5.16 9.12 8.81 8.81 7.28 7.28 10.64 10.64 5.64 7.66 7.51

3.952 3.952 3.952 3.628 3.398 3.773 3.805 3.805 3.721 3.721 3.867 3.867 3.521 3.704 3.670

6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0

J Mol Model (2014) 20:2226 Table 3 Point-wise charges Q (a.u.)

Page 5 of 9, 2226

CN

C

O

N

CC

HC

H

HN

−0.2078

0.5869

−0.5911

−0.4192

−0.0411

0.0173

0.2823

0.1127

distributed computing infrastructures by non-specialists aimed at optimizing the selection of different computing elements for running single and parameter study applications. To execute DL_POLY on the grid environment using PGRADE, we developed a script able to interact with any executable running on a grid node and adaptable to other applications. The whole grid execution process is transparent to the user with the advantage of requiring the evaluation of only the application related failures occurring during the execution. In Fig. 2, a screenshot generated by the visualization Jmol portlet implemented P-GRADE Grid Portal showing the time evolution of the NMA dimer, is illustrated. In order to check the new tools, the DL_POLY version of GEMS was first Intel compiled to the end of assuring binary compatibility when running on the grid environment [49]. The serial version of the DL_POLY package was compiled using the Intel Fortran Compiler, without the need of linking external libraries, due to the fact that the code makes use of its own mathematical routines. For the parallel version of the code, use was made of the same compiler of the MPICH2 library. For both implementations of the code, a static compilation was performed in order to obtain an executable suited to run on the grid. The serial and parallel binary versions were executed on the distributed infrastructure of EGI available to the Virtual Organization COMPCHEM. To the end of estimating the performance of the program on the grid environment, elapsed times were measured using up to eight cores on five EGI sites (though always belonging to the same machine in order to minimize the latency). Results obtained sequentially on one CPU and, in parallel, on 2, 4, and 8 CPUs are shown in Fig. 3, while technical details about the used platforms are given in Table 4. As can be argued from the figure, some of the EGI sites show parallel performances close to the ideal value, implying

Fig. 2 (Left) Output visualization of the time evolution of the NMA dimer produced by the DL_POLY package. (Right) Screenshot of the more stable conformer

a dedicated usage of the processors. Deviations from the ideal value were found to be mainly ascribable to the time-sharing regime as apparent from the data measured on the UNIPG cluster when using two processors. Obviously, also other features like the core-sharing regime adopted for those hyper threaded nodes or the distribution of the calculations among multiple nodes partially occupied by other jobs, and the time needed to transfer the input/output files to/from the grid site are to be taken into account. For this reason, the ability of GriF of adopting alternative distribution strategies turned out to be extremely useful for production runs together with the separate transfer of the output files (that in DL_POLY can be of the order of up to 1 GB). In particular, the use of GriF turned out to be particularly useful for reducing the turn around time caused by network latencies, for checking whether code failures depend either on the unsuitability of one or more starting conditions, or on some undeclared restrictions enforced by the site management. Related tests single out, then, that the different time-sharing regimes adopted by the various clusters lead to different execution times. Such differences, however, were also found to depend on the inadequate distribution of the calculations among nodes partially occupied by other jobs if an appropriate selection of the computing resources is not made. This problem was tackled using the already-mentioned GriF framework that is able to adopt alternative distribution strategies by exploiting the heterogeneity of the grid environment and selecting the machines most appropriate for execution. The advantage of using the GriF selection policy was also exploited, in general, both for reducing the turn around time of the jobs caused by network latencies and for checking whether code failures depend on the unsuitability of one or more starting conditions and for adopting undeclared restrictions enforced by the site management.

2226, Page 6 of 9

J Mol Model (2014) 20:2226

Fig. 3 Elapsed times measured on different EGI Computing Elements available to the VO COMPCHEM plotted as a function of the number of processors used

By leveraging on the use of both P-Grade and GriF we performed a large number of NMA–NMA MD calculations by implementing the subroutine coNMA.f in DL_POLY. The ensemble of two NMA molecules chosen for our calculations was set to be of the NVE type and no boundary conditions were imposed (both for equilibration and for dynamics). All simulations were initiated from the NMA–NMA geometry given in [50] and the system was equilibrated for 10 ps. After achieving equilibration (the equilibrium is reached when the temperature is converged by scaling the velocities of the atoms), trajectories were further integrated for 1.5 ns using a time step of 0.1 ps. Batches of calculations, starting from different initial configurations, were run at values of total energy (Etot) corresponding to temperatures of about 100 K for the purpose of investigating the formation of different isomers as the temperature varies. Then the geometries of the most stable isomers were selected as initial configurations for further runs at decreasing values of Etot (and consequently of the temperature T) down to the limit of T=2 K. By following the integration of the considered trajectories, the distributions of the relevant angles and distances have been analyzed to the end of assessing whether the proposed equilibrium geometries correspond to specific limiting values. In particular, the inspection of the evolution of V as a function

of the relevant distances has been directed to the goal of singling out the existence of different conformers having similar energies. It is also worth noticing here that under the considered conditions, V is mainly characterized by the behavior of VvdW and, therefore, the various effects can be traced back to the VvdW parameters illustrating the effect of the temperature on isomers formation. This is what is shown in Fig. 4, where the values of V calculated in different simulations well illustrate two opposite situations. Panel (a) shows the case of a lower temperature (5 K) corresponding to E tot = −0.31486 eV in which only one conformer is formed. On the contrary panel (b) shows a higher temperature (15 K) corresponding to E tot = −0.36484 eV case in which more than one conformer is formed. The higher temperature plot, in fact, shows separate blocks indicating that the same (or very similar) energy values can be associated with different geometries (i.e., different conformers) which are not necessarily related to the one formed at low temperature. Quantitative estimates of energetic and geometrical properties of the various NMA dimer conformers are given in Table 5 where corresponding ab initio values [28, 53] are given for comparison. The binding energies (De) were calculated by subtracting the energy of two fully optimized monomers from the energy of the dimer. We compare here only the geometry and the energy values estimated for one conformer because the geometries of the other (more stable) conformer

Table 4 Main features of the used EGI sites supporting the COMPCHEM VO Cluster

CPU vendor

CPU clock

RAM

Hellasgrid (GR) UniPg (IT) CESGA (ES) LAL (FR) UOI (GR)

Xeon Xeon Opteron Opteron Opteron

(MHz) 2393 2800 2800 2400 2800

(MB) 2048 4096 1524 2048 2048

Fig. 4 a Simulations carried out at T=5 K in which only one NMA– NMA conformer is formed. b Simulations carried out at T=15 K in which different NMA–NMA conformers are formed

J Mol Model (2014) 20:2226 Table 5 The values of binding energy, De, and selected geometry data for the NMA–NMA system

Page 7 of 9, 2226

De/eV

rOH Å

rNO/Å

−0.3123 −0.2550 −0.2746 −0.2888 −0.3100

1.844 – – 2.01 1.96 1.96

2.860 3.028 – 2.97 2.96 2.96

are not reported in the literature. For instance, Köddermann and Ludwig [28], from density functional calculations at the B3LYP/6-31+G* level of theory, calculated a binding energy (De) of −0.2885 eV (−0.2746 eV when introducing a counterpoise correction). At the same time, Vargas et al. [54], in their MP2/aug-cc-pVDZ level study aimed at optimizing the geometries via single-point calculations, found a lower binding energy value of −0.4193 eVat the best level of calculation that includes BSSE corrections. The results obtained from ab initio calculations were found to differ among them in some cases [51, 54] up to about 60 %. However, authors of [51 and 54] suggest −0.3 eV as a reference value for the NMA dimer binding energy De. The corresponding most stable geometry is a non-planar one, in which the NMA monomers are not parallel to each other.

Conclusions This paper describes in detail how the exploitation of the synergistic nature of the recently developed grid technologies has impacted not only the possibility of carrying out ab initio simulations of complex systems but also the feasibility of collaborative endeavors based on different competences and exploiting innovative tools. In particular, the consequent possibility of further developing GEMS to higher-level complexity applications has been exploited to the end of tackling the study of the aggregation of protein systems (the NMA dimer in our case) at an atomistic level. For this purpose, the functionality of selecting among different packages dealing with either ab initio or empirical (though ab initio validated) portable formulations of the interaction has been adopted. This has implied the tackling of various aspects associated with the synergistic nature of the approach (like standard data formats, accurate analytical formulations of the interaction, appropriate workflows, frameworks suited for platforms selection on the ground of quality evaluators, available scientific portals). In this context, we have developed a distributed implementation of the DL − POLY package on the grid and of the new coNMA routine generalizing the calculation of the potential

rO−HðCH3 Þ /Å

d /degrees NHO

d /degrees COH

2.701 – – – – –

171.3 – – 157.0 168.0 179.6

131.2 – – 122.0 118.0 133.6

Reference

Present [51] [28] [52] [52] [53]

for an aggregate of NMA molecule easily extendible to an arbitrary number of proteins. For the MD calculations reported in this paper, we have focused on the NMA dimer. To this end, the involved non-electrostatic contributions to the total intermolecular interaction has been described in terms of Improved Lennard Jones potentials and a proper set of coulombic terms assigned to various interaction centers and properly placed on each molecule of the aggregate. The scripts developed for distributing DL − POLY using P-GRADE have made it easy to run massive molecular dynamics calculations and to show that the VILJ function the necessary deal of flexibility to the formulation of the van der Waals interaction. In this way, it has also been possible to show that the VILJ function is able to cope with the different electrostatic environments of the NMA–NMA simulations without changing the structural parameters ε and r0 of the various pair potentials. The calculations eventually showed that the C=O-H(CH 3) hydrogen bond, of key importance for biological systems, plays an important stabilizing role when the conformer is formed (see right-hand side of Fig. 2). Hydrogens belonging to the methyl group bonded to the carbonyl carbon atom show directionality in their interactions with the oxygen atom suggestive of the formation of hydrogen bonding. This interaction (of the type of weak hydrogen bonding) appears to be associated with both inter- and intra-chain structural configurations.

Acknowledgments M. Albertí acknowledges financial support from the Ministerio de Educación y Ciencia (Spain, Projects CTQ201016709) and to the Comissionat per a Universitats i Recerca del DIUE (Generalitat de Catalunya, Project 2009-SGR 17). The Centre de Serveis Científics i Acadèmics de Catalunya CESCA-C4 and Fundació Catalana per a la Recerca are also acknowledged for the allocated supercomputing time. Noelia Faginas Lago acknowledges financial support from MIUR PRIN 2008 (contract 2008KJX4SN 003), Phys4entry FP72007-2013 (contract 242311) and EGI Inspire. A. Lombardi also acknowledges financial support to MIUR-PRIN 2010–2011 (contract 2010ERFKXL 002). Thanks are also due to INSTM, CINECA, IGI and the COMPCHEM virtual organization for the allocation of computing time. The research leading to the results presented in this paper has been made possible thanks to the grid resources and services provided by the European Grid Infrastructure (EGI) and the Italian Grid Infrastructure (IGI). For more information, please reference the EGI-InSPIRE paper (http://go.egi.eu/pdnon) and the IGI Web server (http://www.italiangrid.it).

2226, Page 8 of 9

References 1. Laganà A, Costantini A, Gervasi O, Faginas-Lago N, Manuali C, Rampino S (2010) J Grid Comput 8(4):571 2. Laganà A (2005) In: Laganà A, Lendvay G (eds) Theory of chemical reaction dynamics, NATO science series II: mathematics, physics and chemistry, vol. 145 (Springer Netherlands), pp. 363–380. DOI10. 1007/1-4020-2165-8_17 3. Gervasi O, Laganà A (2004) Future Gener. Comp. sy. 20(5), 703. DOI10.1016/j.future.2003.11.028. URL http://www.sciencedirect. com/science/article/pii/S0167739X03002498 4. Laganà A, Garcia E, Paladini A, Casavecchia P, Balucani N (2012) Faraday Discuss. 157, 415. DOI10.1039/C2FD20046E 5. Angeli C, Bendazzoli G, Borini S, Cimiraglia R, Emerson A, Evangelisti S, Maynau D, Monari A, Rossi E, Sanchez-Marin J, Szalay P, Tajti A (2007) Int J Quantum Chem 107:2082 6. Liuthi H, Evangelisti S, Rossi E, Laganà A, Morari A, Rampino S, Verdicchio M, Szalay P, Ruud K, Tajti A, Bendazzoli G, Borini S, Cimiraglia R, Angeli C, Kallay M, Sanchez-Marin J, Scemama A, Baldridge K (2013) J Comp Chem in press 7. Albertí M, Costantini A, Laganà A, Pirani F (2012) J Phys Chem B 116(14):4220. doi:10.1021/jp301124z 8. Lombardi A, Faginas-Lago N, Laganà A, Pirani F, Falcinelli S (2012) In: Murgante B, Gervasi O, Misra S, Nedjah N, Rocha A, Taniar D, Apduhan B (eds) Computational science and its applications ICCSA 2012, lecture notes in computer science, vol. 7333, (Springer Berlin Heidelberg) pp. 387–400. DOI 10.1007/978-3-642-31125-3_30 9. Barreto PRP, Albernaz AF, Caspobianco A, Palazzetti F, Lombardi A, Grossi G, Aquilanti V (2012) Comput Theor Chem 990:56 10. Lombardi A, Palazzetti F, Maciel GS, Aquilanti V, Sevryuk MB (2011) Int J Quantum Chem 111:1651 11. Bartolomei M, Pirani F, Laganà A, Lombardi A (2012) J Comput Chem 33:1806 12. Lombardi A, Faginas-Lago N, Pacifici L, Costantini A (2013) J Phys Chem A 117(45):11430. doi:10.1021/jp408522m 13. Arteconi L, Laganà A (2005) In: Gervasi O, Gavrilova M, Kumar V, Laganà A, Lee H, Mun Y, Taniar D, Tan C (eds) Computational science and its applications ICCSA 2005, lecture notes in computer science vol. 3480 (Springer Berlin Heidelberg) pp. 1093–1100. DOI10.1007/11424758_114 14. Pirani F, Cappelletti D, Liuti G (2001) Chem Phys Lett 350:286 15. Albertí M, Castro A, Laganà A, Pirani F, Porrini M, Cappelletti D (2004) Chem Phys Lett 392:514 16. Albertí M, Castro A, Laganà A, Moix M, Pirani F, Cappelletti D, Liuti G (2005) J Phys Chem A 109:2906 17. Albertí M, Aguilar A, Lucas JM, Pirani F (2009) Theor Chem Acc 123:21 18. Albertí M (2010) J Phys Chem A 114:2266 19. Albertí M, Faginas-Lago N, Laganà A, Pirani F (2011) Phys Chem Chem Phys 13(18):8422 20. Hohenberg P, Kohn W (1964) Phys Rev 136, B864. DOI 10.1103/ PhysRev.136.B864 21. Bernstein FC, Koetzle TF, Williams GJ, Jr EEM, Brice MD, Rodgers JR, Kennard O, Shimanouchi T, Tasumi M (1977) J Mol Biol 112, 535. URL http://www.rcsb.org/pdb/home/home.do 22. Gražulis S, Chateigner D, Downs RT, Yokochi AFT, Quirós M, Lutterotti L, Manakova E, Butkus J, Moeck P, Le Bail A (2009) J Appl Cryst 42(4):726. doi:10.1107/S0021889809016690 23. McQuarrie DA, Simon JD (1997) In: physical chemistry: a molecular approach. University Science, USA 24. Case D, Cheatham T III, Darden T, Gohlke H, Luo R, Merz K Jr, Onufriev A, Simmerling C, Wang B, Woods R (2005) J Comput Chem 26:1668. DOI10.1107/S0021889809016690. URL http://ambermd.org

J Mol Model (2014) 20:2226 25. Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M (1983) J Comput Chem 4(2):187. doi:10.1002/jcc. 540040211 26. Damm W, Frontera A, Tirado Rives J, Jorgensen WL (1997) J Comput Chem 1 8(16):1955. doi:10.1002/(SICI)109 6987X(199712)18:163.0.CO;2-L 27. Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM, Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA (1995) J Am Chem Soc 117(19):5179. doi:10.1021/ja00124a002 28. Koddermann T, Ludwig R (2004) Phys Chem Chem Phys 6:1867. doi:10.1039/B314702A 29. Ataka S, Takeuchi H, Tasumi M (1984) J Mol Struct 113(0):147. DOI10.1016/0022-2860(84)80140-4. URL http://www. sciencedirect.com/science/article/pii/0022286084801404 30. Pirani F, Cappelletti D, Liuti G (2001) Chem Phys Lett 350(34): 286. DOI10.1016/S0009-2614(01)01297-0. URL http://www. sciencedirect.com/science/article/pii/S0009261401012970 31. Pirani F, Albertí M, Castro A, Moix M, Cappelletti D (2004) Chem Phys Lett 394:37 32. Pirani P, Brizi S, Roncaratti L, Casavecchia P, Cappelletti D, Vecchiocattivi F (2008) Phys Chem Chem Phys 10:5489 33. Albertí M, Faginas-Lago N (2013) Eur Phys J D 67:73 34. Faginas-Lago N, Huarte Larrañaga F, Albertí M (2009) Eur Phys J D 55(1):75 35. Albertí M, Faginas-Lago N, Pirani F (2012) Chem Phys 399:232 36. Paolantoni M, Faginas-Lago N, Albertí M, Laganà A (2009) J Phys Chem A 113(52):15100 37. Albertí M, Faginas-Lago N, Pirani F (2011) J Phys Chem A 115(40): 10871 38. Albertí M, Aguilar A, Lucas JM, Pirani F (2010) J Phys Chem A 114(44):11964. doi:10.1021/jp105763h 39. Albertí M, Aguilar A, Lucas JM, Pirani F, Coletti C, Re N (2009) J Phys Chem A 113:14606 40. Capitelli M, Cappelletti D, Colonna G, Gorse C, Laricchiuta A, Liuti G, Longo S, Pirani F (2007) Chem Phys 338:62 41. Hess B, Kutzner C, van der Spoel D, Lindahl E (2008) J Chem Theory Comput 4:435 42. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, Chipot C, Skeel RD, Kalá L, Schulten K (2005) J Comput Chem 26(16):1781. doi:10.1002/jcc.20289 43. S.D. Laboratory, http://www.cse.clrc.ac.uk/ccg/software/ DL_POLY/index.shtml. URL http://cccbdb.nist.gov/. Available on line 44. Kacsuk P, Sipos G (2005) J Grid Comput 3(3–4):221. doi:10.1007/ s10723-005-9012-6 45. Skouteris D, Costantini A, Laganà A, Sipos G, Balaskà K, Kacsuk P (2008) In: Gervasi O, Murgante B, Laganà A, Taniar D, Mun Y, Gavrilova M (eds) Computational science and its applications ICCSA 2008, lecture notes in computer science vol. 5072 (Springer Berlin Heidelberg) pp. 1108–1120. DOI10.1007/978-3-540-698395_84 46. Costantini A, Gutierrez E, Cacheiro J, Rodriguez A, Gervasi O, Laganà A (2010) In: Taniar D, Gervasi O, Murgante B, Pardede E, Apduhan B (eds) Computational science and its applications ICCSA 2010, lecture notes in computer science vol. 6019 (Springer Berlin Heidelberg) pp. 41–52. DOI10. 1007/978-3-642-12189-0_4 47. Costantini A, Gutierrez E, Cacheiro JL, Rodriguez A, Osvaldo Gervasi AL (2010) Int J of Web and Grid Services 6(r27): 141. DOI10.1504/IJWGS.2010.033789. URL http://www.inderscience. com/info/inarticle.php?artid=33789 48. Manuali C, Laganà A, Rampino S (2010) Comp Phys Commun 181(7): 1179. DOI10.1016/j.cpc.2010.03.001. URL h t t p : / / w w w. s c i e n c e d i r e c t . c o m / s c i e n c e / a r t i c l e / p i i / S0010465510000810

J Mol Model (2014) 20:2226 49. ManualiC, Laganà A (2011) Future Gener Comp Sy. 27(3): 315. DOI10.1016/j.future.2010.08.006. URL http://www.sciencedirect. com/science/article/pii/S0167739X10001494 50. Zhang R, Li H, Lei Y, Han S (2004) J Mol Struct 693(13), 17. DOI10. 1016/j.molstruc.2004.01.035. URL http://www.sciencedirect.com/ science/article/pii/S0022286004000791 51. Trabelsi S, Bahri M, Nasr S (2005) J Chem Phys 122(2):024502. doi:10.1063/1.1824035

Page 9 of 9, 2226 52. Qian W, Mirkin NG, Krimm S (1999) Chem Phys Lett 315(12): 125. DOI10.1016/S0009-2614(99)01031-3. URL h t t p : / / w w w. s c i e n c e d i r e c t . c o m / s c i e n c e / a r t i c l e / p i i / S0009261499010313 53. Watson TM, Hirst JD (2002) J Phys Chem A 106(34):7858. doi:10. 1021/jp025551l 54. Vargas R, Garza J, Friesner RA, Stern H, Hay BP, Dixon DA (2001) J Phys Chem A 105(20):4963. doi:10.1021/jp003888m