SciELO - Scientific Electronic Library Online

Home Pagealphabetic serial listing  

Services on Demand




Related links


Boletín de la Sociedad Chilena de Química

Print version ISSN 0366-1644

Bol. Soc. Chil. Quím. vol.46 n.1 Concepción Mar. 2001 


F. D. González-Nilo†, N. E. Robertson†† and M. L. Contreras†*

†University of Santiago de Chile, Faculty of Chemistry and Biology, Chemical Sciences
Department, Casilla 40, Correo 33, Santiago, Chile.
†† University of Tarapacá, Faculty of Sciences, Casilla 7-D, Arica, Chile
(Received: May 31, 2000 - Accepted: December 7, 2000)


A suitable model for molecular dynamics simulation studies of an acrylic acid (AA) polyelectrolyte gel considering explicit water molecules is presented. It consists of a molecule of crosslinker bound to four AA strands completely ionized. Radial distribution functions and nonbonding energy results calculated with both molecular and quantum mechanics methods show that a differentiated solvation occurs for monomer units located at different distances from the crosslinker molecule. The used calculation methodology and the found differentiated solvation could be very interesting for polyelectrolyte gel studies and their applications.

KEYWORDS: poly(acrylic acid), molecular simulation, polyelectrolite gel, radial distribution function (RDF), combined HF/MD calculations.


Se presenta un modelo apto para estudios de simulación de dinámica molecular de un gel polielectrolito de ácido acrílico (AA) considerando moléculas de agua explícitas, el cual consiste en una molécula de entrecruzante unida a cuatro cadenas de AA completamente ionizadas. Los resultados de funciones de distribución radial y energías no-enlazantes obtenidas tanto por métodos de mecánica molecular como de mecánica cuántica muestran que hay una solvatación diferenciada para las unidades monoméricas ubicadas a distintas distancias del entrecruzante. La metodología de cálculo usada y la solvatación diferenciada encontrada podrían ser un aporte interesante para posteriores estudios sobre geles polielectrolitos y sus aplicaciones.

PALABRAS CLAVES: poli(ácido acrílico), simulación molecular, gel polielectrolito, función de distribución radial (FDR), cálculos HF/DM combinados.


Polyelectrolyte gels are infinite networks of ionizable monomers and offer many interesting applications in several areas such as drug delivery and mechanochemical systems [1-5]. These gels show a significant volume change (gel swelling) with external solution composition, temperature, osmotic pressure, or electrical field [6,7]. A study of these amorphous systems whose molecular knowledge is experimentally limited is necessary for understanding gel molecular interactions and from this knowledge to envisage practical applications. Molecular modeling could be an interesting approach for visualising gel behaviour, for instance, in an aqueous media.

Free energy involved in gel swelling could be related to the number of water molecules associated with monomer molecules in gel. This quantity depends on structural monomer properties and crosslinker content. One of the principal aims of this work is to determine the interaction energy of monomer units with associated water. Theoretical calculation of solvation free energy for these systems directly with solvent continuum using quantum mechanics is not possible yet for most laboratories. However combined Molecular Dynamics (MD) and ab initio methods (Hartree-Fock, HF) considering explicit solvent could give a good approximation to solvation free energies because of the combination of an efficient conformational search algorithm with an accurate energy estimation method. Some interesting studies [8-10] have validated the linear response approach for solvation free energy, where solvation free energy correlates linearly with electrostatic molecule-solvent interaction energy (correlation coefficients higher than 0.95). This relation can be expressed [9] as:

(DGsol) @ 1/2 (<DVm – s>)el. When using this approach relative solvation free energies can be calculated from electrostatic nonbonding energies. These energies consider implicitly the entropy contribution and therefore it is not necessary its specific calculation. In this way a quite good estimation of relative solvation free energies for charged systems could be obtained from total nonbonding energies because the electrostatic term predominates in these systems.

MD simulation studies on gels are complicated due to the need of defining an appropriate molecular model for the desired system. For solving this challenge a few studies [11,12] have used a counterion condensation model for linear poly(acrylic acid), PAA. These works [11,12] did not consider the crosslinker in the PAA neither explicit solvent.

In the present work a molecular model considering both explicit water molecules and crosslinker is presented. MD simulation is done using this model. Trajectory analysis is carried out both through radial distribution functions (RDF) and through nonbonding interaction energy calculations between repeating monomer units and associated water. From these values a differenciated solvation was determined. Results obtained in this work are given below.

2. Model and Simulation Details

2.1. Program and Potential FunctionsSimulation programs [13] Insight II and Discover_3, (MSI) run on a SGI Challenge machine with two R4400 co-processors (256 Mb RAM) using CVFF forcefield (Consistent-Valence Forcefield [13]). HF calculations were done on a Pentium III 450 MHz (164 Mb RAM) using Gaussian 98 [14] for Windows NT.

2.2. Molecular Simulation
Molecular simulation of an acrylic acid (AA) gel model was done at 300 K using a 40 Å cubic periodic box with 1,989 water molecules, NVT ensemble (i.e., the simulation is performed at constant number of molecules, volume and temperature), group-based van der Waals and electrostatic cutoffs of 12 Å and 15 Å respectively, and all monomer units completely ionized. A cutoff has to be defined for limiting the number of both van der Waals and coulombic interaction calculations just to a predetermined distance beyond which interactions are not considered. In this way it is possible to decrease CPU time without generating structural perturbations. The system was considered diluted enough to ignore counterions. Strand elasticity was not accounted here because it is not significant for these chain lengths.

The model consists of a macroion molecular structure having one molecule of crosslinker (N,N´-methylenebis(acrylamide)) bound to four AA strands constituted by 7 AA monomer units bonded head-to-tail as shown in Figure 1. First, geometry of this wholly ionized macroion was minimized in vacuum and then it was located in a water periodic box where just the solvent was minimized. Next, a total minimization was done using the steepest descent and conjugate gradient algorithms. Finally, a 150 ps MD simulation (1 fs timestep) was done until solvent was stabilized.

RDF for every monomer residue was calculated using MD trajectory data corresponding to the equilibrium zone (50 - 150 ps). Nonbonding energy calculations (HF/6-31G*), correcting for basis set superposition error (BSSE) were done on isolated systems including second solvation spheres for each monomer. These second solvent shells (i.e. the whole configuration of the system) were defined by the RDF studies. HF calculations for one of these geometries were performed and the necessary hydrogen atoms were added to head and tail carbon atoms in order to complete a precise definition of the system.

3. Results and Discussion

3.1. Model and MD simulations
Figure 1 shows the model used in this work, including (for clarity) only its first water coordination sphere. It considers a similar global behavior for the four strands and allows for focussing energy calculations just on two of the strands and their environment. This helps to save CPU time. Monomer units were numbered according to its relative position from 1 to 7, starting from the crosslinker. MD trajectory obtained for this system is depicted on Figure 2. Systems having counterions require trajectories large enough to assure their stabilization. Present system contains no counterion neither any isolated punctual charge and essentially needs only stabilization of water molecules around the macroion. That should be achieved with trajectories as the one shown on Figure 2.

3.2. RDF analysis
RDF was calculated for carboxylic carbon atom in each monomer and water molecules, from MD trajectory data. A bidimensional representation of radial distribution functions for monomers in first position is depicted in Figure 3a. Figure 3b shows a tridimensional representation of these results as an isocontour map (X = carboxylic carbon atom - water molecules distance; Y = monomer positions and Z = g(r). This quantity represents the probability of finding a water molecule at a distance r). RDF were calculated for all monomers and averaged according to their relative positions.

Figure 3a shows two solvation spheres: one located between 1.4 Å and approximately 2.2-2.4 Å. The second ends between 3.3-3.5 Å. The higher g(r) value observed is located inside the first solvation sphere close to monomers in third position (Figure 3b), at 1.7 Å approx. Another important g(r) value, of lower intensity, is observed around the same distance near monomers 2, 3 and 5. Within the second solvation sphere two high g(r) value zones of equivalent intensities are found at 2.7 Å approximately near monomers 2, 3 and 4. These results show that monomers in third positions are the most solvated ones.

3.3. Nonbonding interaction energies
As a first approach molecular mechanics (MM) was used for calculating nonbonding energies between each monomer unit and the solvent in a periodic box without cutoff and using CVFF force field. This force field has been validated [15] for calculating solvation DG values using LIE (Linear Interaction Energy) approximation [16] in this way demonstrating that the use of LIE approximation with CVFF force field allows for a good estimation of these values. Validation of the found results was done using HF methods. One configuration was taken every 2 ps of MD trajectory within the range 50 - 150 ps for MM calculations and every 10 ps for the HF calculations. MM nonbonding interaction energy calculations were done without cutoff (InsightII package [13]) between two subsets: one of them is the isolated monomer unit and the other consists of all the periodic system water molecules. Analysis of these energies is done by averaging the calculated energies for all equidistant monomer units (same relative distance to the crosslinker; see Figure 4).

On the other side the nonbonding interaction energy was calculated using HF theory which is a more accurated calculation method. For each monomer unit in the model the supermolecule approach was used and a set of three single point HF calculations was done: the monomer in vacuum, the second sphere water molecules, and finally, the whole system (monomer unit and associated water). In total 42 calculations (two strands of seven monomer units each) were done for each configuration taken from MD trajectory data. However, this procedure presents a disadvantage. The calculated interaction energies include a spurious attractive contribution (the basis set superposition error, BSSE) which arises because the components are stabilized in the supermolecule, not for any physical reason, but simply because the basis functions on one component add flexibility to the basis set on the other component. In the enlarged supermolecule basis set the components have lower energies. For weak interaction energies and small basis set the BSSE can dominate the calculated interaction potential and, it can distort our picture of the nature of the forces between the components. It is common to correct for BSSE by the counterpoise method of Boys and Bernardi [17] in which component energies are calculated in the full supermolecule basis set. In this work a BSSE correction was made through the massage command provided by Gaussian 98. Results are presented in Table 1 where BSSE error is given for every monomer. As it was expected the higher superposition errors are produced when a higher number of molecules are interacting at shorter distances like monomers at position 3 which are the most solvated. The BSSE decrease with distance to the crosslinker.

On Figure 4 final results are shown for the average nonbonding energies versus monomer residue relative positions. MM results show a general increase of nonbonding energies with monomer-crosslinker distance. Monomers at positions 2 and 3 present the lowest energies. HF results show that monomer residues at position 3 have the lowest nonbonding energies while monomer residues at position 6 have the highest energies. The average energy difference between these two monomer positions is near 21 kcal/mol. This difference is almost twice for MM values due probably to forcefield parameterization.

Intrachain repulsion in the ionized polymer model, where all carboxylic groups are fully charged, contribute unfavorably to the electrostatic potential energy and the flexibility of the chain segment decreases [11]. In this work 150 ps were enough to stabilize the system as it becomes apparent from Figure 2. At crosslinker neighborhood there is a highly charged zone, which concentrates a great number of water molecules. So this effect could favor solvation of monomers nearby crosslinker. The fact that monomers at position 1 are less solvated than monomers at positions 2 and 3 can be explained by the steric hindrance produced at the neighborhood of the crosslinker. Residues at position 2 and 3 have more available space for interacting with environmental water molecules. As a result, lower nonbonding energies are found for these monomers, probably due to the formation of a greater number of hydrogen bonds. The effect produced by the crosslinker decreases with distance. In this way monomers at positions 6 and 7 have the highest energy. Monomers at positions 7 are the last in the strands and have lower HF energies because they interact with a larger number of water molecules (frontier effect). Both methods show similar results: monomers nearest crosslinker are the best solvated. MM values does not show the frontier effect found at HF for unit 7 due to the use of periodic boundary conditions. Extension of these calculations, using a higher number of monomers for each strand should show that solvation energy should not vary significantly beyond a given distance, obtaining a plateau.


The simulation model and the methodology presented in this work allowed for the calculation of nonbonding interaction energies between each monomer unit in the polyelectrolyte and the solvent. These energies allow for a good approximation of solvation free energies implicitly considering entropy contribution as validated with experimental data [15]. RDF studies allowed for a clear definition of two solvation spheres for each monomer unit. This definition let a quantum mechanics evaluation of the nonbonding energies to be carried out. Also nonbonding energy interactions were evaluated by MM methods. Both methods showed the same trend. So, this methodology constitutes a study tool and a contribution to molecular knowledge analysis of polyelectrolyte gels and could be easily extended to other systems.

MD simulation RDF analysis of a crosslinked acrylate macroion and both HF and MM evaluations of nonbonding interaction energies showed a differentiated solvation for each monomer unit at different crosslinker distances (i.e. crosslinker effect favours solvation of nearest monomer units). This fact could help to increase the spectra of possible applications for these gels.


Thanks are given to Dr. Raúl E. Cachau at the National Cancer Institute, Frederick, USA; Dr. Gloria Cárdenas-Jirón at University of Santiago de Chile and Dr. Alejandro Toro-Labbé at P. Catholic University of Chile, for collaboration and fruitful discussions. Also to CONICYT and the University of Santiago de Chile for F. D. González-Nilo Fellowships and FONDECYT Project # 2970090.

Presented in part at the Centennial Meeting of the American Physical Society, APS, March 1999, Atlanta, Georgia, U.S.A.


1. T. Miyata, N. Asami, T. Uragami, Nature 399 (1999) 766.        [ Links ]

2. C. Bell, N. Peppas, Biopolymer II serie: Adv. in Polym. Sci. 122 (1995) 125.        [ Links ]

3. Scranton, B. Rangarajan, J. Klier, Biopolymer II, serie: Adv. in Polym. Sci. 122 (1995), 1.        [ Links ]

4. P. Stayton, T. Shimoboji, C. Long, A. Chilkoti, G. Chen, J. Harris, A. Hofman, Nature, 378:6556 (1995) 472.        [ Links ]

5. H. Okuzaki, Y. Osada, Electrochim. Acta, 40 (1995) 13.        [ Links ]

6. Y. Osada, J. P. Gong, Adv. Mater. 10, (1998) 827.        [ Links ]

7. S. Beltran, H. H. Hooper, H. W. Blanch, J. M. Prausnitz, J. Chem. Phys. 92 (1990) 3, 2061.        [ Links ]

8. Roux, H-A. Yu, M. Karplus, J. Phys. Chem. 94 (1990) 4683.        [ Links ]

9. J. Åqvist, C. Medina, J. E. Samuelsson, Protein Eng. 7 (1994) 385.        [ Links ]

10. J. Åqvist, J. Phys. Chem. 94 (1990) 8021.        [ Links ]

11. M. Ravi, A. Hopfinger, J. Comput. Polym. Sci. 4 (1994) 19.         [ Links ]

12. M. Ravi, A. Hopfinger, J. Comput. Polym. Sci. 3 (1993) 137.        [ Links ]

13. Insight II and Discover_3 (V97.0), San Diego, Molecular Simulations Inc. (1998).        [ Links ]

14. Gaussian 98, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, V. G. Zakrzewski, J. A. Montgomery, R. E. Stratmann, J. C. Burant, S. Dapprich, J. M. Millam, A. D. Daniels, K. N. Kudin, M. C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G. A. Petersson, P.Y. Ayala, K. Cui, K. Morokuma, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. Cioslowski, J. V. Ortiz, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, C. González, M. Challacombe, P. M. W. Gill, B. G. Johnson, W. Chen, M. W. Wong, J. L. Andres, M. Head-Gordon, E. S. Replogle, and J. A. Pople, Gaussian Inc., Pittsburgh, PA, (1998).        [ Links ]

14. F. D. González-Nilo, Ph. D. Thesis, University of Santiago de Chile, 2000;         [ Links ] F. D. González-Nilo and M. L. Contreras, TMMeC, Vol 2, 1998-1999; ISSN 0797-9274.        [ Links ]

16. J. Åqvist, J. Comput. Chem. 17 (1996) 1587.        [ Links ]

17. S. F. Boys, F. Bernardi, Mol. Phys. 19 (1970) 553.        [ Links ]

Creative Commons License All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License