A Coarse-Grained Force Field Parameterized for ... - ACS Publications

Subscriber access provided by NEW YORK UNIV

Article

A Coarse-Grained Force Field Parameterized for MgCl and CaCl Aqueous Solutions 2

2

Zheng Gong, and Huai Sun J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00206 • Publication Date (Web): 28 Jun 2017 Downloaded from http://pubs.acs.org on June 29, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

A Coarse-Grained Force Field Parameterized for MgCl2 and CaCl2 Aqueous Solutions Zheng Gong and Huai Sun*

School of Chemistry and Chemical Engineering, and Ministry of Education Key Laboratory of Scientific and Engineering Computing, Shanghai Jiao Tong University, Shanghai, 200240, China

Keywords: calcium chloride, magnesium chloride, coarse-grained force field, hydration shell, osmotic coefficient

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract

Calcium and magnesium ions play important roles in many physicochemical processes. To facilitate the investigation of phenomena related to these ions that occur over large length- and time-scales, a coarse-grained force field (CGFF) is developed for MgCl2 and CaCl2 aqueous solutions. The ions are modeled by CG beads with characteristics of hydration shells. To accurately describe the non-ideal behavior of the solutions, osmotic coefficients in a wide range of concentration were used as guidance for parameterization. The osmotic coefficients were obtained from the chemical potential increments of water calculated using Bennett acceptance ratio (BAR) method. The result CGFF accurately reproduces experimental osmotic coefficients, densities, surface tensions and cation-anion separations of calcium chloride and magnesium chloride solutions at molalities up to 3.0 mol/kg. As a preliminary application, the force field is applied to simulate aggregations of sodium dodecyl sulfate (SDS) in calcium chloride solution, and the simulation results are consistent with experimental observations.


Page 2 of 50

Page 3 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.

Introduction

Magnesium and calcium ions are the third and fifth most abundant dissolved ions in seawater. Calcium ions plays important roles in environmental science, providing an important link between tectonics, climate and the carbon cycle. Calcium and magnesium are essential elements for living organisms. Movement of calcium ions through the cytoplasm signals many cellular processes. Magnesium ions drive the folding of RNA into stable states and participate in metabolism by binding to adenosine triphosphate. On the other side, calcium and magnesium ions have undesirable effects in consumer and industry application such as forming limescale in hot-water tubes and precipitating with detergents. Because of their omnipresence and importance, calcium and magnesium ions have been the subject of many experiments. Vibrational sum frequency generation spectroscopy1, X-ray absorption fine structure spectroscopy2, neutron diffraction3, X-ray diffraction4-6, Raman spectroscopy7 and nuclear magnetic resonance8 have been widely applied to investigate the interactions between calcium or magnesium ions and water molecules or other solutes in aqueous solutions. Molecular dynamics (MD) simulation is a powerful tool to study the behavior of ions interacting with other molecules in aqueous solution. Both ab initio and all-atom force field (AAFF) MD simulations have been used to study divalent electrolyte solutions. The hydration structure9, transport properties10 of divalent ions, and the binding between



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

divalent ions and surfactant11 or protein12 have been studied using MD simulations. These simulations revealed the mechanisms of ionic interactions. However, due to limitations in time and length scales, previous studies employing quantum and classical simulations have not focused on slow physiochemical processes such as self-aggregation or phase separation. In order to extend the length- and time-scales of simulations on phenomenon involving divalent ions, a coarse-grained force field (CGFF) for calcium and magnesium chlorides is desired. To the best of our knowledge, three CGFFs that cover calcium ion are UNRES (united-residue)13, 3SPN (three sites per nucleotide)14 and MARTINI force field15,16. However, UNRES and 3SPN does not use explicit water and MARTINI treats the calcium ion using the same Lennar-Jones parameters as the sodium ion without validation. In this work we develop a CGFF for calcium and magnesium chlorides with explicit water models. We chose the coarse-grained model of Shinoda-DeVane-Klein (SDK)17,18 as the base for development. The SDK force field has been parameterized and validated primarily for surfactant solutions including ionic species Cl- and Na+. The three-molecule water bead in SDK force field has been parameterized to reproduce the density, compressibility and surface tension of pure water. Thus, the focus of this study is on divalent cations (Ca2+ and Mg2+) and their interactions with water and anions (Cl-). In order to describe the non-ideal behavior of high-concentration solutions of calcium chloride and magnesium chloride, we included osmotic coefficients in the


Page 4 of 50

Page 5 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


training set for the parameterization, as simply fitting the force field parameters to the hydration free energy of ions in infinitely dilute solutions19-22 fails represent non-ideal behavior23. The osmotic coefficient provides information on the combined effect of solute−solute, solute−solvent, and solvent−solvent interactions in a system. Because experimental osmotic coefficients are available for many molecular system in the literature, using these data provides a means of parametrizing the strengths of solute−solute interactions for a wide range of molecule types. An approach to calculate osmotic pressure directly using MD simulation with virtual semi-permeable wall model has been proposed by Murad and Powles24. Since Luo and Roux25 optimized ion−ion interactions in NaCl and KCl aqueous solutions, this method has been used to reparametrize other ions26-30 and carbohydrates31,32. Very recently, Miller et al.33 re-optimized several popular protein force fields using the same approach. In this work, we took a different approach by calculating the osmotic coefficients from chemical potential of solvent (water). The CG beads for ions are assumed hydrated with associated water molecules. We choose to use a modified Lennard-Jones (LJ) function and a distance-dependent dielectric constant34 to present the hydration shell structure. The initial parameters were derived from all-atom force field, and then optimized by simultaneously reproducing osmotic coefficients, liquid densities and surface tensions in a broad range of concentrations. The result CGFF is demonstrated to be capable of predicting the non-ideal properties of CaCl2



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 50

and MgCl2 Solutions. This paper is organized as following: In section 2, we present the hydrated CG model and parameterization protocol. The parameterization and validation results and discussions are presented in section 3. As a preliminary application, this force field has been used to simulate aggregations of SDS in calcium chloride solution, the results are presented in section 4. The conclusion and remarks are given in section 5.

2.

Hydrated Ion Model, Parameterization and Simulation

In the SDK force field17, the water is represented by an uncharged bead (W) each represents three water molecules. The ions are treated as charged and hydrated particles. A chloride ion bound with two water molecules is represented by a single coarse-grained bead denoted as “CLA”18. In this study we add two beads: “CA”, which represents a calcium ion associated with three water molecules, and “MG”, a magnesium ion bound to three water molecules. The four CG beads considered in this work, and the ions/molecules they represent, are illustrated in Figure 1. The general Lennard-Jones function with variable exponents (LJ-m-n) is used to represent the van der Waals (VDW) interactions between the CG beads.

= ()( ) [( ) − ( ) ]

(1)

The LJ-12-4 function is used for interactions between water beads (W-W), as well as ion-water (CA-W, MG-W, CLA-W), cation-cation (CA-CA, MG-MG) and anion-anion


Page 7 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(CLA-CLA) interactions. A modified LJ-6-4+g function (m-LJ) is used to represent cation-anion van der Waals interactions (CA-CLA, MG-CLA) in which the “+g” represents a Gaussian function34 added to the LJ-6-4 function. = + exp −

( )

(2)

The Gaussian represents the shell structure, and prevents the approach between oppositely charged ions with a tunable degree of penetration. The constant H in Eq. 2 represents the resistance to this penetration, and and ! represent the location and width, respectively, of the hydration shells corresponding to each ion-pair. The electrostatic energy is represented by the standard Coulomb function with fixed (formal) charges: "#$ = +2.0 e, "%& = +2.0 e, "#$ = -1.0 e. *

'(') = +,

-

./ .

(3)

,

For interactions between the same ion type (CA-CA, MG-MG, CLA-CLA), the dielectric constant of bulk water 0 = 78 is used in Eq. 3. However, this constant does not apply for interactions between oppositely charged ions separated by short distances where not enough water molecules to shield the electrostatic interactions between them. A distance-dependent function34 for the dielectric constant is applied for interactions between anions and cations: 0=

,3 4,5

+

,3 ,5

tanh

: :

(4)

where 0; is the dielectric constant of bulk water (78), 0< is a constant representing the limiting value of the dielectric constant as the ions are close. We set 0< to 5.2 as



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 50

suggested by Lenart et. al.34 The variables , and !, are the location of the midpoint and the range of transition of dielectric constant curve, respectively, which are subject to optimization. A correction term to the electrostatic energy is evaluated as '('),)> =

*

./ . ,? @,

+,- ,?

− 1B

Since this correction term approaches nearly zero soon ( : :

(5) ,? ,

− 1 = 0.001 when

= 3.4), it is truncated at the same cutoff length as VDW interaction. Then the total

energy of the system is = ∑(H* + ∑ + ∑'(') + ∑'('),)>

(6)

The parameters for the pair-interactions between CLA-CLA and CLA-W beads are taken from the SDK force field. Six new sets of parameters are optimized in this work for interactions between CA-CLA, MG-CLA, CA-W, MG-W, CA-CA, and MG-MG pairs. Among them the most critical ones are those belonging to the cation-anion interactions, CA-CLA and MG-CLA. An initial guess for each of the parameters was obtained by fitting them against the potential of mean force (PMF) generated from all-atom simulations. Two sets of all-atom force fields were used in PMF calculations: AAFF-1 and AAFF-2. The AAFF-1 used the cation parameters suggested by Mamatkulov et al.35 and the anion parameters by Smith and Dang36. AAFF-2 used the cation parameters by Aqvist19 and the anion parameters of Chandrasekhar et. al.37 The SPC/E parameters were adopted for water38. These two all-atom force fields generated different PMF, which will be discussed in the following section. The PMF calculated using AAFF-1 were used to


Page 9 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


estimate initial parameters, prior to more rigorous optimization, between anions and cations in our CG force field. The parameters were then optimized in a trial-and-error fashion to match experimental data of osmotic coefficients, densities and surface tensions at different ion concentrations. The osmotic coefficient I is calculated from the chemical potential ∗ difference of water between pure solvent J; and solution J; by: * N N ∗

I = − LM %? ∑O? ?

5

(7)

where P; is the molar mass of water, Q< is the molality of each solute. The density and surface tensions were calculated using molecular dynamics simulations as reported previously39. The optimized parameters are given in Table 1. There are a total of 22 adjustable parameters. Forty-two (42) experimental data points were used to optimize the adjustable parameters of the force field. Constant pressure and temperature (NPT) simulations at 298 K and 1 bar were performed in most cases. Constant volume and temperature (NVT) simulations at 298 K were performed for the surface tension calculations. MD simulations using AAFF were performed using GROMACS 4.6.740. The cutoff of van der Waals and Coulomb forces was set to 12 Å, the continuum correction41 and the particle mesh Ewald (PME) methods42 were used for calculating long-range interactions. The time step was set to 2 fs. The SETTLE algorithm43 was used to constrain the bond lengths and angle in water molecules. The temperature and pressure were controlled using the velocity-rescaling



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

thermostat44 and Berendsen barostat45, respectively. MD simulation using CGFF were performed using LAMMPS46. The cutoff for van der Waals and Coulomb forces was set to 15 Å. Long-range electrostatics were treated using the particle–particle particle–mesh (PPPM) solver47. The time step was set to 10 fs. Smaller steps (2 fs, 5 fs) were tested, but resulted in nearly identical results. The temperature and pressure were controlled using the Nose-Hoover thermostat48,49 and barostat50,51 respectively. Simulation boxes with 3-D periodic boundary conditions and different simulation durations were used for different purposes. 0.5 ns was used for equilibration of all systems. For liquid density calculations, a cubic box with edge length of 35 Å was used in conjunction with a duration of 1.0 ns for data collection. For surface tension calculations, a box with X- and Y- dimensions of 50 Å, and a Z-dimension of 150 Å was used to construct a liquid-vapor slab model. For these simulations, production time of 20 ns was used. For osmotic coefficient calculation, a cubic box with an edge-length of 45 Å was used, in which the production duration was 10 ns, respectively. The excess free energy required to determine the osmotic coefficient was calculated using the BAR method52 implemented in pyMBAR53. The soft-core interaction54 is used to avoid singularities when dismissing non-bonded interactions. For PMF calculations, a cubic box with an edge length of 40 Å was used. The distance between two ionic species was sampled using the umbrella harmonic potentials, evenly distributed with spacing of 0.3 Å and a force constant of 30 kcal/mol/Å2 between 2.5 to 7.0 Å, and a spacing of 0.6 Å with


Page 10 of 50

Page 11 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


force constant of 10 kcal/mol/Å2 between 7.6 to 14.8 Å. The production duration was 10 ns for each window. The software Plumed 2.1.055 was used to perform the umbrella sampling. The PMF was extracted using weighted histogram analysis method (WHAM)56 implemented by Grossfield57. The software PACKMOL58 was used to build the initial coordinates of the simulation boxes. For each simulation, block averages or multiple independent runs from different initial configurations were carried out in order to estimate statistical uncertainties.

3.

Results and Discussions

The potential energy profile between calcium (Ca2+) and chloride (Cl-) ions, resulting from the optimized force field, is presented in Figure 2. The van der Waals (m-LJ) and electrostatic (Elec.) contributions to the total potential energy (PE) are depicted separately. The subtle balance between hydration repulsion and strong electrostatic attraction at short distances promotes the formation of two energy minima on the PE curve, separated by a barrier at 3.4 Å, indicating the existence of contact ion pair (CIP) near 2.9 Å and solvent separated ion pair (SSIP) near 5.2 Å in calcium chloride solution. The existence of two classes of ion pairs is pronounced in the PMF curve which exhibits a double-well feature. The PMF curve has two minima located at 2.8 A and 5.0 A, separated by a barrier of 1.7 kcal/mol. The PMF is lower than the PE at short distance, higher than PE at medium distance and converges to the PE at long distance. The lower PMF at short distance is due



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

to the entropy increase from released water beads from the shells as the two solutes approach each other. By analyzing the trajectory of a molecular dynamics simulation, we found that 24 water beads are in the first shells of the two ions when they are separated at a distance of 10.0 Å, while the number is reduced to 16 when the ions are separated at a distance of 2.8 Å. A water bead is defined to be in the first shell of an ion when the distance between the ion and a water is less than 6.4 Å, which corresponds to the first minima of the RDF between the ion and water beads. The potential energy between magnesium (Mg2+) and chloride (Cl-) ions is presented in Figure 3. The shape of the vdW and electrostatic potential energy contributions resemble those between Ca2+ and Cl-. However, the total potential energy curve only shows one minima located at 4.9 Å, corresponding to the SSIP. There is no energy minimum corresponding to the CIP. This reflects the stronger hydration for smaller cation Mg2+. The PMF curve shows a shoulder at short distance at approximately 3 Å and a minimum at 4.5 Å. Again, the PMF is lower than the PE at short distance, higher than the PE at medium distance, and converges to the PE at long distance. For comparison, the PMFs between Ca2+ and Cl-, Mg2+ and Cl- respectively are calculated using two composite all-atom force fields (AAFF-1 and AAFF-2) and the results are shown in Figure 4. The parameters between counter ions in the AAFF-1 are specifically optimized targeting the activity coefficient35, while in AAFF-2 the parameters between counter ions were generated using Lorentz-Berthelot combination rules. In the


Page 12 of 50

Page 13 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


case of CaCl2, these two force fields produced similar PMF curves but the first free energy minima differ by 1.2 kcal/mol. For MgCl2, the PMFs calculated using the two force fields are dramatically different. The PMF of AAFF-1 shows no sign of CIP but AAFF-2 shows a CIP minima at 2.6 Å. This shows the ion pair structure is highly sensitive to the force field parameters of the ions. Figure 5 shows the radial distribution functions between the cations and chloride anions calculated at a concentration of 1.0 mol/kg using the optimized CGFF. For calcium chloride, a small peak at 2.8 Å shows the presence of CIP and a larger and broader peak located at 5.0 Å represents the SSIP. This results support two X-ray diffraction experiments4,5 that reveal the existence of CIP at distance of about 2.75 Å and SSIP at distance of 4.86-5.01 Å. For magnesium chloride, CIP has not been observed in our CG simulations. There is a single peak located at 4.7 Å, representing the SSIP. Due to the small radius of Mg2+, water molecules bind to Mg2+ strongly and the complex rarely breaks away from the first hydration shell. This phenomena has been revealed by both experiments8 and ab initio simulation59. Raman spectroscopy7 and an X-ray diffraction6 experiments establish the absence of CIP and the presence of SSIP at around 4.5 Å. We also calculated RDFs between cations which show the first peak at around 9 Å for either Mg2+-Mg2+ or Ca2+-Ca2+ pairs. This confirms that short-range corrections in the LJ function and distance-dependent dielectric constant are not necessary to describe these ion pairs.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 50

Figure 6 shows the simulated and experimental60,61 densities of the electrolyte solutions of CaCl2 and MgCl2 versus the molality from 0.5 to 3.5 mol/kg. The data are obtained by 1 ns data collected after the equilibrium runs in the NPT simulations. The uncertainties evaluated using five blocks are less than 0.002 g/mL. Excellent agreement is obtained over the wide range of concentrations simulated for both electrolytes. RS and The chemical potential of water J; can be split into two parts, the ideal J; 'T the excess J; contributions: RS 'T J; = J; + J;

(8)

RS The ideal contribution J; depends on the number density of water ρ in solution and

can be calculated analytically: RS J; = VWln(ρYZ; /"; )

(9)

Here, Y; is the thermal de Broglie wavelength, and "; is the partition function of 'T water molecule. The excess contribution J; can be calculated by MD simulations. Since

a CG water bead “W” consists of three water molecules, a thermodynamic cycle as 'T from shown in Figure 7 is used to calculate the excess chemical potential of water J;

the excess chemical potential of CG water-bead “W”: 'T J; =

]&^_`a,bac` 4]&de,f ]&de,g

(10)

Z

Here, the hydration free energy of CG water-bead hiO'jS,

kSj

was calculated using the

BAR method and the value obtained in pure water is -8.44 kcal/mol. The coarse-graining free energy change hi#&,l or hi#&,( can be divided into two parts. One can be


Page 15 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


calculated as the association free energy of trimer from the monomer and trimer molecular partition functions. Another includes all factors that are not in the association free energy, such as the configuration entropy of the trimer (the number of states of the monomers in a trimer). Assuming the trimers are the same in gas and liquid states, the second part can be canceled out so that we replace the difference between hi#&,l and hi#&,( by the difference in association free energies: hi#&,l − hi#&,( ≈ hij

A Coarse-Grained Force Field Parameterized for ... - ACS Publications

Recommend Documents